快速傅里叶算法(FFT)示例

作者：追风剑情 发布于：2018-2-23 18:25 分类：Algorithms

连续非周期信号的傅里叶变换公式

111.png(公式中的i是虚数单位)

现在假设在x(t)的某一段连续区间上以周期T进行采样，得到N个采样点，则每个采样点的离散傅里叶变换公式就是：

2222.png(n=0,1,...,N-1)

令33333.png,则上式可简单记为:

44444.png(n=0,1,...,N-1)

上式可用线性方程组表示为:

111111.png 

for (int n=0; n<N; n++)
	for (int k=0; k<N; k++)
		int p = n * k;//旋转因子指数


逆变换公式

111111.png

11111.png被称为旋转因子。

77.png的周期性可以表示为5555.png

77.png的对称性可以表示为6666.png

77.png的可约性可以表示为22222.png

原始信号与两个分组信号之间的关系:

每次分解，都将原始信号x(n)按照时间顺序(也就是n的序号)分成奇、偶两个组x1(r)和x2(r)，其中r与n的关系是:

当n为偶数时，令n=2r;

当n为奇数时，令n=2r+1;

也就是说，原始信号与两个分组的信号存在以下关系:

x(2r)=x1(r), x(2r+1)=x2(r)  其中r=0,1,...,N/2-1

将一个N点DFT分解为两个N/2点DFT：

88888.png

因为9999.png，代入上式，得到：

00000.png

由周期性可知:111.png,因此：

2222.png

同理可得：3333.png

由以上关系可知：

一个N点DFT变换的前N/2个周期数据的转换关系是:

4444.png (n=0,1,...,N/2-1)

一个N点DFT变换的后N/2个周期数据的转换关系是:

5555.png(n=0,1,...,N/2-1)

2222.png

欧拉公式

3333.png

从上面的公式可知，FFT是递归求解过程。不断对数据进行奇偶分组，如下图:

11111.png

蝶形算法流程图(借用网上一张截图)

输出倒位序的FFT流程图.png

11111.png

旋转因子WNP规律
P = i*(1<<L); (i为第几次蝶形运算, L为当前级数(从0开始算))

111111.png

FFT算法的复量运算量为

111111.png

图像处理：先进行行变换，再对变换结果进行列变换，或者先进行列变换，再对结果进行行变换。

二维离散傅里叶变换
二维离散傅里叶变换(C#版)

示例 


  1. using System;
  2. using System.Collections.Generic;
  3.  
  4. namespace ConsoleApp34
  5. {
  6.     class Program
  7.     {
  8.         static void Main(string[] args)
  9.         {
  10.             TestFFT();
  11.             TestFFT2();
  12.             Console.Read();
  13.         }
  14.  
  15.         private static void TestFFT()
  16.         {
  17.             const int N = 8;
  18.             //初始化时域数据
  19.             Complex[] TD2FD = new Complex[N];
  20.             for (int i = 0; i < N; i++)
  21.             {
  22.                 Complex cpx = new Complex();
  23.                 cpx.re = i;
  24.                 cpx.im = 0;
  25.                 TD2FD[i] = cpx;
  26.             }
  27.             Console.WriteLine("------ 一维原始信号 -----");
  28.             Print(TD2FD);
  29.  
  30.             Console.WriteLine("----- FFT -----");
  31.             FFT(TD2FD);
  32.             Print(TD2FD);
  33.  
  34.             Console.WriteLine("----- IFFT -----");
  35.             IFFT(TD2FD);
  36.             Print(TD2FD);
  37.         }
  38.  
  39.         private static void TestFFT2()
  40.         {
  41.             const int W = 8, H = 8;
  42.             Complex2D complex2D = new Complex2D(W, H);
  43.             for (int i=0; i<H; i++)
  44.             {
  45.                 for (int j=0; j<W; j++)
  46.                 {
  47.                     Complex cpx = new Complex();
  48.                     cpx.re = i * W + j;
  49.                     cpx.im = 0;
  50.                     complex2D.SetComplex(i, j, cpx);
  51.                 }
  52.             }
  53.  
  54.             Console.WriteLine("------ 二维原始信号 -----");
  55.             complex2D.Print();
  56.  
  57.             Console.WriteLine("----- FFT2 -----");
  58.             FFT2(complex2D);
  59.             complex2D.Print();
  60.  
  61.             Console.WriteLine("----- IFFT2 -----");
  62.             IFFT2(complex2D);
  63.             complex2D.Print();
  64.         }
  65.  
  66.         //快速傅里叶变换
  67.         public static void FFT(Complex[] TD2FD)
  68.         {
  69.             FFT_Core(TD2FD, WT_LUT(TD2FD.Length, 1));
  70.         }
  71.  
  72.         //快速傅里叶变换 (二维)
  73.         public static void FFT2(Complex2D TD2FD)
  74.         {
  75.             //对每一行做FFT
  76.             for (int i=0; i<TD2FD.Height; i++)
  77.             {
  78.                 Complex[] row = TD2FD.GetRow(i);
  79.                 FFT(row);
  80.             }
  81.  
  82.             //对每一列做FFT
  83.             for (int i = 0; i < TD2FD.Width; i++)
  84.             {
  85.                 Complex[] column = TD2FD.GetColumn(i);
  86.                 FFT(column);
  87.             }
  88.         }
  89.  
  90.         //快速傅里叶逆变换
  91.         public static void IFFT(Complex[] FD2TD)
  92.         {
  93.             //做FFT变换
  94.             Complex[] WT = WT_LUT(FD2TD.Length, -1);
  95.             FFT_Core(FD2TD, WT);
  96.             //除以N
  97.             for (int i = 0; i < FD2TD.Length; i++) {
  98.                 FD2TD[i].re /= FD2TD.Length;
  99.                 FD2TD[i].im /= FD2TD.Length;
  100.             }
  101.         }
  102.  
  103.         //快速傅里叶逆变换 (二维)
  104.         public static void IFFT2(Complex2D FD2TD)
  105.         {
  106.             //对每一行做IFFT
  107.             for (int i = 0; i < FD2TD.Height; i++)
  108.             {
  109.                 Complex[] row = FD2TD.GetRow(i);
  110.                 IFFT(row);
  111.             }
  112.  
  113.             //对每一列做IFFT
  114.             for (int i = 0; i < FD2TD.Width; i++)
  115.             {
  116.                 Complex[] column = FD2TD.GetColumn(i);
  117.                 IFFT(column);
  118.             }
  119.         }
  120.  
  121.         // 返回旋转因子查询表(twiddle factor lookup table)
  122.         private static Complex[] WT_LUT(int N, int flag = 1)
  123.         {
  124.             Complex[] WT = new Complex[N];
  125.             for (int i = 0; i < N; i++)
  126.             {
  127.                 Complex cpx_wt = new Complex();
  128.                 float angle = (float)(-i * Math.PI * 2 / N);
  129.                 cpx_wt.re = (float)Math.Cos(angle);
  130.                 //IFFT flag=-1, FFT flag=1
  131.                 cpx_wt.im = flag * (float)Math.Sin(angle);
  132.                 WT[i] = cpx_wt;
  133.             }
  134.             return WT;
  135.         }
  136.  
  137.         private static void FFT_Core(Complex[] TD2FD, Complex[] WT)
  138.         {
  139.             int power = (int)Math.Log(TD2FD.Length, 2);
  140.             int butterfly;
  141.             int p, s;
  142.             Complex x1, x2, wt;
  143.  
  144.             //倒位排序
  145.             BitReverse(TD2FD);
  146.  
  147.             //蝶形运算
  148.             for (int k = 0; k < power; k++) //级数
  149.             {
  150.                 for (int j = 0; j < 1 << (power - k - 1); j++) //组数
  151.                 {
  152.                     //每组有几个元素
  153.                     butterfly = 1 << k + 1;
  154.                     //计算参与蝶形运算的两个元素的索引
  155.                     p = j * butterfly;
  156.                     s = p + butterfly / 2;
  157.                     for (int i = 0; i < butterfly / 2; i++) //蝶形运算次数
  158.                     {
  159.                         x1 = TD2FD[i + p];
  160.                         x2 = TD2FD[i + s];
  161.                         wt = WT[i * TD2FD.Length / butterfly];
  162.                         TD2FD[i + p] = x1 + x2 * wt;
  163.                         TD2FD[i + s] = x1 - x2 * wt;
  164.                     }
  165.                 }
  166.             }
  167.         }
  168.  
  169.         private static int BitReverse(int x)
  170.         {
  171.             //倒位排序
  172.             //0   1   2   3   4   5   6   7   十进制
  173.             //000 001 010 011 100 101 110 111 十进制对应的二进制
  174.             //000 100 010 110 001 101 011 111 码位反转
  175.             //0   4   2   6   1   5   3   7   码位反转后对应的十进制
  176.             int[] table = new int[8] { 0, 4, 2, 6, 1, 5, 3, 7 };
  177.             return table[x];
  178.         }
  179.  
  180.         // 倒位排序——雷德算法
  181.         private static void BitReverse(Complex[] array)
  182.         {
  183.             int i, j, k;
  184.             int N = array.Length;
  185.             Complex temp;
  186.             j = 0;
  187.  
  188.             for (i = 0; i < N - 1; i++)
  189.             {
  190.                 if (i < j)
  191.                 {
  192.                     temp = array[i];
  193.                     array[i] = array[j];
  194.                     array[j] = temp;
  195.                 }
  196.                 // 求j的下一个倒序位
  197.                 // N/2的二进制最高位为1，其他位都为0
  198.                 // 判断最高位是否为1，可与N/2进行比较
  199.                 // 判断次高位是否为1，可与N/4进行比较
  200.                 k = N >> 1;
  201.                 //如果k<=j,表示j的最高位为1
  202.                 while (k <= j)
  203.                 {
  204.                     //当k<=j时，需要将最高位变为0
  205.                     j = j - k;
  206.                     //判断次高位是否为1,依次类推，逐个比较，直到j某个位为0
  207.                     k >>= 1;
  208.                 }
  209.                 j = j + k;//将0变为1
  210.             }
  211.         }
  212.  
  213.         // 打印
  214.         private static void Print(Complex[] TD2FD)
  215.         {
  216.             for (int i = 0; i < TD2FD.Length; i++)
  217.             {
  218.                 Console.WriteLine("(re={0}, im={1})", TD2FD[i].re, TD2FD[i].im);
  219.             }
  220.             Console.WriteLine();
  221.         }
  222.     }
  223.  
  224.     //定义复数
  225.     public class Complex
  226.     {
  227.         public float re;//实数部
  228.         public float im;//虚数部
  229.  
  230.         // 幅值
  231.         public double Amplitude
  232.         {
  233.             get {
  234.                 return Math.Sqrt(re*re + im*im);
  235.             }
  236.         }
  237.  
  238.         // 相位
  239.         public double Phase
  240.         {
  241.             get {
  242.                 return Math.Atan2(im, re);
  243.             }
  244.         }
  245.  
  246.         public static Complex operator +(Complex lhs, Complex rhs)
  247.         {
  248.             Complex result = new Complex();
  249.             result.re = lhs.re + rhs.re;
  250.             result.im = lhs.im + rhs.im;
  251.             return result;
  252.         }
  253.  
  254.         public static Complex operator -(Complex lhs, Complex rhs)
  255.         {
  256.             Complex result = new Complex();
  257.             result.re = lhs.re - rhs.re;
  258.             result.im = lhs.im - rhs.im;
  259.             return result;
  260.         }
  261.  
  262.         public static Complex operator *(Complex lhs, Complex rhs)
  263.         {
  264.             Complex result = new Complex();
  265.             result.re = lhs.re * rhs.re - lhs.im * rhs.im;
  266.             result.im = lhs.re * rhs.im + lhs.im * rhs.re;
  267.             return result;
  268.         }
  269.  
  270.         public static Complex operator *(float lhs, Complex rhs)
  271.         {
  272.             Complex result = new Complex();
  273.             result.re = lhs * rhs.re;
  274.             result.im = lhs * rhs.im;
  275.             return result;
  276.         }
  277.  
  278.         public static Complex operator *(Complex lhs, float rhs)
  279.         {
  280.             Complex result = new Complex();
  281.             result.re = lhs.re * rhs;
  282.             result.im = lhs.im * rhs;
  283.             return result;
  284.         }
  285.     }
  286.  
  287.     public class Complex2D
  288.     {
  289.         private List<Complex[]> rows = new List<Complex[]>();
  290.         private List<Complex[]> columns = new List<Complex[]>();
  291.         private int m_width;
  292.         private int m_height;
  293.  
  294.         public Complex2D(int width, int height)
  295.         {
  296.             m_width = width;
  297.             m_height = height;
  298.             for (int i = 0; i < height; i++)
  299.                 rows.Add(new Complex[width]);
  300.             for (int i = 0; i < width; i++)
  301.                 columns.Add(new Complex[height]);
  302.         }
  303.  
  304.         public int Width { get { return m_width; } }
  305.         public int Height { get { return m_height; } }
  306.  
  307.         public Complex[] GetRow(int i)
  308.         {
  309.             return rows[i];
  310.         }
  311.  
  312.         public Complex[] GetColumn(int i)
  313.         {
  314.             return columns[i];
  315.         }
  316.  
  317.         public void SetRow(int i, Complex[] src)
  318.         {
  319.             rows[i] = src;
  320.         }
  321.  
  322.         public void SetColumn(int i, Complex[] src)
  323.         {
  324.             columns[i] = src;
  325.         }
  326.  
  327.         public void SetComplex(int i, int j, Complex src)
  328.         {
  329.             rows[i][j] = src;
  330.             columns[j][i] = src;
  331.         }
  332.  
  333.         public void SetComplexs(Complex[][] src)
  334.         {
  335.             for (int i=0; i<src.Length; i++)
  336.             {
  337.                 Complex[] row = src[i];
  338.                 for (int j = 0; j < row.Length; j++)
  339.                     SetComplex(i, j, row[j]);
  340.             }
  341.         }
  342.  
  343.         public void Print()
  344.         {
  345.             for (int i=0; i<rows.Count; i++)
  346.             {
  347.                 Complex[] row = rows[i];
  348.                 for (int j = 0; j < row.Length; j++)
  349.                     Console.Write("{0:G} ", row[j].re.ToString().PadRight(5));
  350.                 Console.WriteLine();
  351.             }
  352.             Console.WriteLine();
  353.         }
  354.     }
  355. }


运行测试
111111.png 


  1. //频谱分析窗函数类型
  2. public enum FFTWindow
  3. {
  4. 	//
  5. 	// 矩形窗 (使用最多的窗)
  6. 	// W[n] = 1.0.
  7. 	Rectangular = 0,
  8. 	//
  9. 	// 三角窗
  10. 	// W[n] = TRI(2n/N).
  11. 	Triangle = 1,
  12. 	//
  13. 	// 汉宁窗(余弦窗)
  14. 	// W[n] = 0.54 - (0.46 * COS(n/N) ).
  15. 	Hamming = 2,
  16. 	//
  17. 	// 海明窗(升余弦窗)
  18. 	// W[n] = 0.5 * (1.0 - COS(n/N) ).
  19. 	Hanning = 3,
  20. 	//
  21. 	// 布拉克曼
  22. 	// W[n] = 0.42 - (0.5 * COS(nN) ) + (0.08 * COS(2.0 * nN) ).
  23. 	Blackman = 4,
  24. 	//
  25. 	// 哈布斯窗
  26. 	// W[n] = 0.35875 - (0.48829 * COS(1.0 * nN)) + (0.14128 * COS(2.0 * nN)) - (0.01168 * COS(3.0 * n/N)).
  27. 	BlackmanHarris = 5
  28. }


更多示例: Matlib——二维傅里叶变换

标签: Algorithms

« C#调用Lua函数 | Lnaguage Intergrated Query(Linq) 集成查询语言»
Powered by emlog  蜀ICP备18021003号-1   sitemap

川公网安备 51019002001593号