C#實現FFT（遞迴法）
1. C#實現複數類
我們在進行訊號分析的時候，難免會使用到複數。但是遺憾的是，C#沒有自帶的複數類，以下提供了一種複數類的構建方法。
複數相比於實數，可以理解為一個二維數，構建複數類，我們需要實現以下這些內容：
- 複數實部與虛部的屬性
- 複數與複數的加減乘除運算
- 複數與實數的加減乘除運算
- 複數取模
- 複數取相位角
- 尤拉公式（即\(e^{ix+y}\)）
C#實現的程式碼如下：
public class Complex
{
double real;
double imag;
public Complex(double x, double y) //建構函式
{
this.real = x;
this.imag = y;
}
//通過屬性實現對複數實部與虛部的單獨檢視和設定
public double Real
{
set { this.real = value; }
get { return this.real; }
}
public double Imag
{
set { this.imag = value; }
get { return this.imag; }
}
//過載加法
public static Complex operator +(Complex c1, Complex c2)
{
return new Complex(c1.real + c2.real, c1.imag + c2.imag);
}
public static Complex operator +(double c1, Complex c2)
{
return new Complex(c1 + c2.real, c2.imag);
}
public static Complex operator +(Complex c1, double c2)
{
return new Complex(c1.Real + c2, c1.imag);
}
//過載減法
public static Complex operator -(Complex c1, Complex c2)
{
return new Complex(c1.real - c2.real, c1.imag - c2.imag);
}
public static Complex operator -(double c1, Complex c2)
{
return new Complex(c1 - c2.real, -c2.imag);
}
public static Complex operator -(Complex c1, double c2)
{
return new Complex(c1.real - c2, c1.imag);
}
//過載乘法
public static Complex operator *(Complex c1, Complex c2)
{
double cr = c1.real * c2.real - c1.imag * c2.imag;
double ci = c1.imag * c2.real + c2.imag * c1.real;
return new Complex(Math.Round(cr, 4), Math.Round(ci, 4));
}
public static Complex operator *(double c1, Complex c2)
{
double cr = c1 * c2.real;
double ci = c1 * c2.imag;
return new Complex(Math.Round(cr, 4), Math.Round(ci, 4));
}
public static Complex operator *(Complex c1, double c2)
{
double cr = c1.Real * c2;
double ci = c1.Imag * c2;
return new Complex(Math.Round(cr, 4), Math.Round(ci, 4));
}
//過載除法
public static Complex operator /(Complex c1, Complex c2)
{
if (c2.real == 0 && c2.imag == 0)
{
return new Complex(double.NaN, double.NaN);
}
else
{
double cr = (c1.imag * c2.imag + c2.real * c1.real) / (c2.imag * c2.imag + c2.real * c2.real);
double ci = (c1.imag * c2.real - c2.imag * c1.real) / (c2.imag * c2.imag + c2.real * c2.real);
return new Complex(Math.Round(cr, 4), Math.Round(ci, 4)); //保留四位小數後輸出
}
}
public static Complex operator /(double c1, Complex c2)
{
if (c2.real == 0 && c2.imag == 0)
{
return new Complex(double.NaN, double.NaN);
}
else
{
double cr = c1 * c2.Real / (c2.imag * c2.imag + c2.real * c2.real);
double ci = -c1 * c2.imag / (c2.imag * c2.imag + c2.real * c2.real);
return new Complex(Math.Round(cr, 4), Math.Round(ci, 4)); //保留四位小數後輸出
}
}
public static Complex operator /(Complex c1, double c2)
{
if (c2 == 0)
{
return new Complex(double.NaN, double.NaN);
}
else
{
double cr = c1.Real / c2;
double ci = c1.imag / c2;
return new Complex(Math.Round(cr, 4), Math.Round(ci, 4)); //保留四位小數後輸出
}
}
//建立一個取模的方法
public static double Abs(Complex c)
{
return Math.Sqrt(c.imag * c.imag + c.real * c.real);
}
//建立一個取相位角的方法
public static double Angle(Complex c)
{
return Math.Round(Math.Atan2(c.real, c.imag), 6);//保留6位小數輸出
}
//過載字串轉換方法,便於顯示覆數
public override string ToString()
{
if (imag >= 0)
return string.Format("{0}+i{1}", real, imag);
else
return string.Format("{0}-i{1}", real, -imag);
}
//尤拉公式
public static Complex Exp(Complex c)
{
double amplitude = Math.Exp(c.real);
double cr = amplitude * Math.Cos(c.imag);
double ci = amplitude * Math.Sin(c.imag);
return new Complex(Math.Round(cr, 4), Math.Round(ci, 4));//保留四位小數輸出
}
}
2. 遞迴法實現FFT
以下的遞迴法是基於奇偶分解實現的。
奇偶分解的原理推導如下：
\[\begin{split}
X(k)=DFT[x(n)]&=\sum_{n=0}^{N-1}x(n)W_N^{nk}\\
&=\sum_{r=0}^{N/2-1}x(2r)W_N^{2rk}+\sum_{r=0}^{N/2-1}x(2r+1)W_N^{(2r+1)k}，將x(n)按奇偶分解\\
&=\sum_{r=0}^{N/2-1}x(2r)(W_N^{2})^{rk}+W_N^k\sum_{r=0}^{N/2-1}x(2r+1)(W_N^2)^{rk}
\end{split}
\]
\(x(2r)\)和\(x(2r+1)\)都是長度為\(N/2-1\)的資料序列，不妨令
\[x_1(n)=x(2r)\\
x_2(n)=x(2r+1)
\]
則原來的DFT就變成了：
\[\begin{split}
X(k)&=\sum_{n=0}^{N/2-1}x_1(n)(W_N^{2})^{nk}+W_N^k\sum_{n=0}^{N/2-1}x_2(n)(W_N^2)^{nk}\\
&=F(x_1(n))+W_N^kF(x_2(n))\\
&=X_1(k)+W_N^kX_2(k)
\end{split}
\]
於是，將原來的N點傅立葉變換變成了兩個N/2點傅立葉變換的線性組合。
但是，N/2點傅立葉變換隻能確定N/2個頻域資料，另外N/2個資料怎麼確定呢？
因為\(X_1(k)\)和\(X_2(k)\)週期都是\(N/2\)，所以有
\[X_1(k+N/2)=X_1(k),X_2(k+N/2)=X_2(k)\\
W_N^{k+N/2}=-W_N^k\\
\]
從而得到：
\[\begin{split}
X(k+N/2)&=X_1(k+N/2)+W_N^{k+N/2}X_2(k+N/2)\\
&=X_1(k)-W_n^kX_2(k)
\end{split}
\]
綜上，我們就可以得到遞迴法實現FFT的流程：
-
對於每組資料，按奇偶分解成兩組資料
-
兩組資料分別進行傅立葉變換，得到\(X_1(k)\)和\(X_2(k)\)
-
總體資料的\(X(k)\)由下式確定：
\[X(k)==X_1(k)+W_N^kX_2(k)\\ X(k+N/2)=X_1(k)-W_n^kX_2(k)\\ 0\le k \le N/2 -1 \] -
對上述過程進行遞迴
具體程式碼實現如下：
public Complex[] FFTre(Complex[] c)
{
int n = c.Length;
Complex[] cout = new Complex[n];
if (n == 1)
{
cout[0] = c[0];
return cout;
}
else
{
double n_2_f = n / 2;
int n_2 = (int)Math.Floor(n_2_f);
Complex[] c1 = new Complex[n / 2];
Complex[] c2 = new Complex[n / 2];
for (int i = 0; i < n_2; i++)
{
c1[i] = c[2 * i];
c2[i] = c[2 * i + 1];
}
Complex[] c1out = FFTre(c1);
Complex[] c2out = FFTre(c2);
Complex[] c3 = new Complex[n / 2];
for (int i = 0; i < n / 2; i++)
{
c3[i] = new Complex(0, -2 * Math.PI * i / n);
}
for (int i = 0; i < n / 2; i++)
{
c2out[i] = c2out[i] * Complex.Exp(c3[i]);
}
for (int i = 0; i < n / 2; i++)
{
cout[i] = c1out[i] + c2out[i];
cout[i + n / 2] = c1out[i] - c2out[i];
}
return cout;
}
}
3. 補充：窗函式
順便提供幾個常用的窗函式：
- Rectangle
- Bartlett
- Hamming
- Hanning
- Blackman
public class WDSLib
{
//以下窗函式均為periodic
public double[] Rectangle(int len)
{
double[] win = new double[len];
for (int i = 0; i < len; i++)
{
win[i] = 1;
}
return win;
}
public double[] Bartlett(int len)
{
double length = (double)len - 1;
double[] win = new double[len];
for (int i = 0; i < len; i++)
{
if (i < len / 2) { win[i] = 2 * i / length; }
else { win[i] = 2 - 2 * i / length; }
}
return win;
}
public double[] Hamming(int len)
{
double[] win = new double[len];
for (int i = 0; i < len; i++)
{
win[i] = 0.54 - 0.46 * Math.Cos(Math.PI * 2 * i / len);
}
return win;
}
public double[] Hanning(int len)
{
double[] win = new double[len];
for (int i = 0; i < len; i++)
{
win[i] = 0.5 * (1 - Math.Cos(2 * Math.PI * i / len));
}
return win;
}
public double[] Blackman(int len)
{
double[] win = new double[len];
for (int i = 0; i < len; i++)
{
win[i] = 0.42 - 0.5 * Math.Cos(Math.PI * 2 * (double)i / len) + 0.08 * Math.Cos(Math.PI * 4 * (double)i / len);
}
return win;
}
}