傅里葉變換在信息處理應用中具有很實用的價值,而快速傅里葉變換,即FFT,是實用的計算算法。
本文介紹FFT和2維FFT的C#實現方法。
1. FFT編程依據
FFT是按照如圖結構(也稱蝶形結構)進行運算(圖片來源於網絡)。
圖中,箭頭代表數據流向,箭頭與箭頭的合並點代表相加,箭頭下面的常數代表相乘,WN(注:N為下標) = exp(-j * 2 * PI / N),j為虛數單位,DIT的輸入(或DIF的輸出)序列重排點y[m]與初始輸入點x[n]的關系為:m與n互為倒位序(例:001 -- 100,011 -- 110)。
該圖給出的是8點FFT例程,其余2^N點FFT可按照這個圖進行類推,類推后可以發現,2^N點FFT只需用三個循環套嵌即可實現。
對於idft的算法,根據公式X[k] = 對n求和{x[n] * (WN) ^ kn}以及x[n] = 對k求和{X[k] * (WN) ^ (-kn)} /N可知,只需將FFT算法中的WN換為(1 / WN)並將最后結果除以N即可得到idft。
2. FFT算法實現
由於作者並不擅長面向對象編程,僅僅是在功能上實現了本文需求,代碼在閱讀、效率上有很多缺陷。
下面是FFT算法的實現,完整的實現在文章后面。
// 快速傅里葉變換 // 作者:死貓 // ComplexList src:數據輸入 // ref int lenght:欲變換數據的長度,函數調用后此值更改為實際變換長度 // int flag:區分fft或dtft,為1為fft,為-1為idft private static ComplexList fft_core(ComplexList src, ref int lenght, int flag) { //按時間抽取FFT方法(DIT) //補零后長度 int relog2N = ReLog2N(lenght); int bitlenghth = relog2N; int N = 0x01 << bitlenghth; //重新復制數據,同時進行 // 逆序排放,並補零 int index; ComplexList Register = new ComplexList(N); for (int i = 0; i < N; i++) { index = ReArrange(i, bitlenghth); Register[i] = (index < src.Lenght) ? src[index] : new Complex(0); } //生成WN表,以免運行時進行重復計算 ComplexList WN = new ComplexList(N / 2); for (int i = 0; i < N / 2; i++) { WN[i] = new Complex(Math.Cos(2 * Math.PI / N * i), -flag * Math.Sin(2 * Math.PI / N * i)); } //蝶形運算 int Index0, Index1; Complex temp; for (int steplenght = 2; steplenght <= N; steplenght *= 2) { for (int step = 0; step < N / steplenght; step++) { for (int i = 0; i < steplenght / 2; i++) { Index0 = steplenght * step + i; Index1 = steplenght * step + i + steplenght / 2; temp = Register[Index1] * WN[N / steplenght * i]; Register[Index1] = Register[Index0] - temp; Register[Index0] = Register[Index0] + temp; } } } //若為idft if (-1 == flag) { for (int i = 0; i < Register.Lenght; i++) { Register[i] /= N; } } //賦值 lenght = N; /* //清理內存 WN = null; temp = null; GC.Collect(); */ //返回 return Register; }
3. 二維FFT變換
2維傅里葉變換是在1維傅里葉變換基礎上實現的,實現方法為:
1. 對2維輸入數據的每一行進行FFT變換,變換結果仍然按行存入一個二維數組中;
2. 對按行FFT變換后的結果,對每一列進行FFT變換,變換結果仍然按列存入一個二維數組中,該數組即為2維FFT變換結果。
代碼實現如下:
// 2維快速傅里葉變換 // 作者:死貓 // int width, int height:欲變換數據的長度和寬度,函數調用后此值更改為實際變換長度 // int flag:區分fft或dtft,為1為fft,為-1為idft private static Complex2D fft_2D_core(Complex2D src, ref int width, ref int height, int flag) { //補零后長度 int width_Log2N = ReLog2N(width); int height_Log2N = ReLog2N(height); int Relog2N = Math.Max(width_Log2N, height_Log2N); int ReWidth = 0x01 << Relog2N; int ReHeight = 0x01 << Relog2N; //重新復制數據,清零 Complex2D ReList2D = new Complex2D(ReWidth, ReHeight); ReList2D.Clear(); int width_temp = Math.Min(src.Width, width); int height_temp = Math.Min(src.Height, height); for (int i = 0; i < height_temp; i++) { for (int j = 0; j < width_temp; j++) { ReList2D[i, j] = new Complex(src[i, j].Re, src[i, j].Im); } } ComplexList Xn; ComplexList Xk; int Lenght_temp; //第1遍fft for (int i = 0; i < ReHeight; i++) { Xn = ReList2D.GetRow(i); Lenght_temp = Xn.Lenght; Xk = fft_core(Xn, ref Lenght_temp, flag); ReList2D.SetRow(i, Xk); } //第2遍fft for (int i = 0; i < ReWidth; i++) { Xn = ReList2D.GetColumn(i); Lenght_temp = Xn.Lenght; Xk = fft_core(Xn, ref Lenght_temp, flag); ReList2D.SetColumn(i, Xk); } //賦值 width = ReWidth; height = ReHeight; //清理內存 Xn = null; Xk = null; GC.Collect(); //返回 return ReList2D; }
3. 二維FFT測試結果
我對一張圖片進行低通、高通濾波,其傳輸函數分別為:
低通:Hlp(z) = (1 - a) / 2 * (1 + 1 / z) / (1 - a * 1 / z),a = 0.8;
高通:Hhp(z) = (1 + a) / 2 * (1 - 1 / z) / (1 - a * 1 / z),a = 0.7。
低通、高通濾波測試結果分別如下:
4. 其余相關代碼實現
由於作者並不擅長面向對象編程,僅僅是在功能上實現了本文需求,代碼在閱讀、效率上有很多缺陷。
代碼1:
// 快速傅里葉變換 // ComplexList src:數據輸入,若長度非2的n次冪, // 則函數對末尾補零后的數據進行運算 public static ComplexList fft(ComplexList src) { int lenght = src.Lenght; return fft_core(src, ref lenght, 1); } // 快速傅里葉變換 // ComplexList src:數據輸入 // ref int lenght:欲變換的長度,函數調用后此值更改為實際變換長度 public static ComplexList fft(ComplexList src, ref int lenght) { return fft_core(src, ref lenght, 1); } // 快速傅里葉變換的逆變換 // ComplexList src:數據輸入 public static ComplexList idft(ComplexList src) { int lenght = src.Lenght; return idft(src, ref lenght); } // 快速傅里葉變換的逆變換 // ComplexList src:數據輸入 // ref int lenght:欲變換數據的長度,函數調用后此值更改為實際變換長度 public static ComplexList idft(ComplexList src, ref int lenght) { return fft_core(src, ref lenght, -1); } public static Complex2D fft_2D(Complex2D src) { int width = src.Width; int height = src.Height; return fft_2D_core(src, ref width, ref height, 1); } public static Complex2D idft_2D(Complex2D src) { int width = src.Width; int height = src.Height; return fft_2D_core(src, ref width, ref height, -1); }
代碼2:
// 獲取按位逆序,bitlenght為數據長度 // fft函數內使用 private static int ReArrange(int dat, int bitlenght) { int ret = 0; for (int i = 0; i < bitlenght; i++) { if (0 != (dat & (0x01 << i))) ret |= ((0x01 << (bitlenght - 1)) >> i); } return ret; } // 獲取擴展長度后的冪次 // 由於fft要求長度為2^n,所以用此函數來獲取所需長度 public static int ReLog2N(int count) { int log2N = 0; uint mask = 0x80000000; for (int i = 0; i < 32; i++) { if (0 != ((mask >> i) & count)) { if ((mask >> i) == count) log2N = 31 - i; else log2N = 31 - i + 1; break; } } return log2N; }
代碼3:
public class Complex2D { double[] _Complex2D_Re; double[] _Complex2D_Im; public int Width { get; private set; } public int Height { get; private set; } public Complex this[int Row, int Column] { get { return new Complex(_Complex2D_Re[Row * Width + Column], _Complex2D_Im[Row * Width + Column]); } set { _Complex2D_Re[Row * Width + Column] = ((Complex)value).Re; _Complex2D_Im[Row * Width + Column] = ((Complex)value).Im; } } public Complex2D(int width, int height) { Width = width; Height = height; int lenght = Width * Height; _Complex2D_Re = new double[lenght]; _Complex2D_Im = new double[lenght]; } public void Clear() { Array.Clear(_Complex2D_Re, 0, _Complex2D_Re.Length); Array.Clear(_Complex2D_Im, 0, _Complex2D_Im.Length); } public ComplexList GetColumn(int index) { ComplexList ret = new ComplexList(Height); for(int i = 0; i < Height; i++) { ret[i] = this[i, index]; } return ret; } public int SetColumn(int index, ComplexList src) { for (int i = 0; i < Height; i++) { this[i, index] = (i < src.Lenght) ? src[i] : new Complex(0); } return 0; } public ComplexList GetRow(int index) { ComplexList ret = new ComplexList(Width); for(int i = 0; i < Width; i++) { ret[i] = this[index, i]; } return ret; } public int SetRow(int index, ComplexList src) { for (int i = 0; i < Width; i++) { this[index, i] = (i < src.Lenght) ? src[i] : new Complex(0); } return 0; } public Complex2D GetAmplitude() { Complex2D ret = new Complex2D(Width, Height); for (int i = 0; i < Height; i++) { for (int j = 0; j < Width; j++) { ret[i, j] = new Complex(this[i, j].Modulus()); } } return ret; } } public class ComplexList { double[] _ComplexList_Re; double[] _ComplexList_Im; public int Lenght { get; private set; } public Complex this[int Index] { get { return new Complex(_ComplexList_Re[Index], _ComplexList_Im[Index]); } set { _ComplexList_Re[Index] = ((Complex)value).Re; _ComplexList_Im[Index] = ((Complex)value).Im; } } public ComplexList(int lenght) { Lenght = lenght; _ComplexList_Re = new double[Lenght]; _ComplexList_Im = new double[Lenght]; } public ComplexList(double[] re) { Lenght = re.Length; _ComplexList_Re = re; _ComplexList_Im = new double[Lenght]; } public ComplexList(double[] re, double[] im) { Lenght = Math.Max(re.Length, im.Length); if (re.Length == im.Length) { _ComplexList_Re = re; _ComplexList_Im = im; } else { _ComplexList_Re = new double[Lenght]; _ComplexList_Im = new double[Lenght]; for (int i = 0; i < re.Length; i++) _ComplexList_Re[i] = re[i]; for (int i = 0; i < im.Length; i++) _ComplexList_Im[i] = im[i]; } } public void Clear() { Array.Clear(_ComplexList_Re, 0, _ComplexList_Re.Length); Array.Clear(_ComplexList_Im, 0, _ComplexList_Im.Length); } public double[] GetRePtr() { return _ComplexList_Re; } public double[] GetImPtr() { return _ComplexList_Im; } public ComplexList Clone() { return new ComplexList((double[])(_ComplexList_Re.Clone()), (double[])(_ComplexList_Im.Clone())); } public ComplexList GetAmplitude() { double[] amp = new double[Lenght]; for (int i = 0; i < Lenght; i++) { amp[i] = this[i].Modulus(); } return new ComplexList(amp); } } public class Complex { public double Re; public double Im; public Complex() { Re = 0; Im = 0; } public Complex(double re) { Re = re; Im = 0; } public Complex(double re, double im) { Re = re; Im = im; } public double Modulus() { return Math.Sqrt(Re * Re + Im * Im); } public override string ToString() { string retStr; if (Math.Abs(Im) < 0.0001) retStr = Re.ToString("f4"); else if (Math.Abs(Re) < 0.0001) { if (Im > 0) retStr = "j" + Im.ToString("f4"); else retStr = "- j" + (0 - Im).ToString("f4"); } else { if (Im > 0) retStr = Re.ToString("f4") + "+ j" + Im.ToString("f4"); else retStr = Re.ToString("f4") + "- j" + (0 - Im).ToString("f4"); } retStr += " "; return retStr; } //操作符重載 public static Complex operator +(Complex c1, Complex c2) { return new Complex(c1.Re + c2.Re, c1.Im + c2.Im); } public static Complex operator +(double d, Complex c) { return new Complex(d + c.Re, c.Im); } public static Complex operator -(Complex c1, Complex c2) { return new Complex(c1.Re - c2.Re, c1.Im - c2.Im); } public static Complex operator -(double d, Complex c) { return new Complex(d - c.Re, -c.Im); } public static Complex operator *(Complex c1, Complex c2) { return new Complex(c1.Re * c2.Re - c1.Im * c2.Im, c1.Re * c2.Im + c2.Re * c1.Im); } public static Complex operator *(Complex c, double d) { return new Complex(c.Re * d, c.Im * d); } public static Complex operator *(double d, Complex c) { return new Complex(c.Re * d, c.Im * d); } public static Complex operator /(Complex c, double d) { return new Complex(c.Re / d, c.Im / d); } public static Complex operator /(double d, Complex c) { double temp = d / (c.Re * c.Re + c.Im * c.Im); return new Complex(c.Re * temp, -c.Im * temp); } public static Complex operator /(Complex c1, Complex c2) { double temp = 1 / (c2.Re * c2.Re + c2.Im * c2.Im); return new Complex((c1.Re * c2.Re + c1.Im * c2.Im) * temp, (-c1.Re * c2.Im + c2.Re * c1.Im) * temp); } }
==================
20200208補充源工程,下載地址:
https://files.cnblogs.com/files/Si-Mao/FFT.rar