﻿//------------------------------------------------------------------------------
//	係数計算用メソッド
//
//		作成者：三上直樹，2023/03/30 作成，©三上直樹 2023
//------------------------------------------------------------------------------

using System;

using static System.Math;

namespace FIR_Design_Windowing
{
	public partial class Form1
	{		
		private void Design()
		{
			hk_ = new double[nOrder_+1];

			Func<double, int, double> win = (order, n) => 1.0;	// rectangle 窓
			switch (windowType_.Checked)
			{
				case "Hanning":
					win = (order, n) => 0.5 - 0.5*Cos(2.0*PI*n/order); break;
				case "Hamming":
					win = (order, n) => 0.54 - 0.46*Cos(2.0*PI*n/order); break;
				case "Blackman":
					win = (order, n) => 0.42 - 0.5*Cos(2.0*PI*n/order)
											 + 0.08*Cos(4.0*PI*n/order); break;
				case "Kaiser":
					win = (order, n) => Kaiser(order, n); break;
			}
			
			// LPF の係数を求める
			double fc = fcLow_;
			switch (passBand_.Checked)
			{
				case "HPF": fc = fSample_/2 - fc; break;
				case "BPF":
				case "BRF": fc = (fcHigh_ - fcLow_)/2; break;
			}
			for (int n=0; n<=nOrder_; n++)
			{
				double nSh = n - nOrder_/2.0;
				if (n != nOrder_/2.0)
					hk_[n] = (1.0/(PI*(nSh)))*Sin(2.0*PI*nSh*fc/fSample_)*win(nOrder_, n);
				else
					hk_[n] = 2.0*fc/fSample_;
			}

			if (passBand_.Checked == "HPF")		// HPF の場合
			{
				double pm1 = ((nOrder_ % 4 == 0) || ((nOrder_ - 1) % 4 == 0)) ? 1 : -1;
				for (int n=0; n<=nOrder_; n++) hk_[n] = pm1*Pow(-1, n)*hk_[n];
			}

			if ((passBand_.Checked == "BPF") ||	(passBand_.Checked == "BRF"))	// BPF, BRF の場合
			{
				double w0 = PI*(fcLow_ + fcHigh_)/fSample_;
				for (int n=0; n<=nOrder_; n++)	// BPF の係数
					hk_[n] = 2.0*Cos((n - nOrder_/(double)2)*w0)*hk_[n];

				if (passBand_.Checked == "BRF")	// BRF の場合
					for (int n=0; n<=nOrder_; n++)
						if (n != nOrder_/2.0)	hk_[n] = -hk_[n];
						else					hk_[n] = 1.0 - hk_[n];
			}
		}

		// Kaiser 窓
		private double Kaiser(double order, int n)
		{
			double alpha = 0.0;
			if (att_ >= 50.0) alpha = 0.1102*(att_ - 8.7);
			else if (att_ >= 21) alpha = 0.5842*Pow(att_ - 21.0, 0.4) + 0.07886*(att_ - 21.0);
			double fm = 2.0*(n - order/2.0)/order;
			return I0(alpha*Sqrt(1.0 - fm*fm))/I0(alpha);
		}

		// 0 次の第１種変形ベッセル関数
		private double I0(double x)
		{
			double sum = 1.0;
			double xj  = 1.0;
			for (int j=1; j<20; j++)
			{
				xj = 0.5*xj*x/j;
				sum += xj*xj;
			}
			return sum;
		}
	}
}
