快速傅里叶算法(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)

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

77.png的对称性可以表示为6666.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

示例


using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace Test11
{
    /// <summary>
    /// 快速傅里叶算法测试(FFT)
    /// </summary>
    class Program
    {
        static void Main(string[] args)
        {
            //初始化时域数据
            Complex[] WT = new Complex[8];
            Complex[] TD2FD = new Complex[8];
            for (int i = 0; i < TD2FD.Length; i++) {
                Complex cpx = new Complex();
                cpx.re = 1;
                cpx.im = 1;
                TD2FD[i] = cpx;

                //初始化旋转因子
                Complex cpx_wt = new Complex();
                float angle = (float)(-i * Math.PI * 2 / TD2FD.Length);
                cpx_wt.re = (float)Math.Cos(angle);
                cpx_wt.im = (float)Math.Sin(angle);
                WT[i] = cpx_wt;
            }
            FFT(TD2FD, WT);
            Console.Read();
        }

        public static void FFT(Complex[] TD2FD, Complex[] WT)
        {
            int power = (int)Math.Log(TD2FD.Length, 2);
            int butterfly;
            int p, s;
            Complex x1, x2, wt;

            for (int k = 0; k < power; k++) //级数
            {
                for (int j = 0; j < 1 << k; j++) //组数
                {
                    butterfly = 1 << (power - k);//每组有几个元素
                    //计算参与蝶形运算的两个元素的索引
                    p = j * butterfly;
                    s = p + butterfly / 2;
                    Console.WriteLine("butterfly={0}, p={1}, s={2}", butterfly, p, s);
                    for (int i = 0; i < butterfly / 2; i++) //蝶形运算次数
                    {
                        Console.Write("  ({0}x{1} wtIdx={2})", i + p, i + s, i * (1 << k));
                        x1 = TD2FD[i + p];
                        x2 = TD2FD[i + s];
                        wt = WT[i*(1<<k)];
                        TD2FD[i + p] = x1 + x2 * wt;
                        TD2FD[i + s] = x1 - x2 * wt;
                    }
                    Console.WriteLine();
                }
            }

            //重新排序
            for (int i = 0; i < TD2FD.Length; i++)
            {
                int r = BitReverse(i);
                if (r > i)
                {
                    Complex t = TD2FD[r];
                    TD2FD[r] = TD2FD[i];
                    TD2FD[i] = t;
                }
            }

            for (int i = 0; i < TD2FD.Length; i++)
            {
                Console.WriteLine("(re={0}, im={1})", TD2FD[i].re, TD2FD[i].im);
            }
        }

        static int BitReverse(int x)
        {
            //0   1   2   3   4   5   6   7   十进制
            //000 001 010 011 100 101 110 111 十进制对应的二进制
            //000 100 010 110 001 101 011 111 码位反转
            //0   4   2   6   1   5   3   7   码位反转后对应的十进制
            int[] table = new int[8] { 0, 4, 2, 6, 1, 5, 3, 7 };
            return table[x];
        }
    }

    //定义复数
    public class Complex
    {
        public float re;//实数部
        public float im;//虚数部

        public static Complex operator +(Complex lhs, Complex rhs)
        {
            Complex result = new Complex();
            result.re = lhs.re + rhs.re;
            result.im = lhs.im + rhs.im;
            return result;
        }

        public static Complex operator -(Complex lhs, Complex rhs)
        {
            Complex result = new Complex();
            result.re = lhs.re - rhs.re;
            result.im = lhs.im - rhs.im;
            return result;
        }

        public static Complex operator *(Complex lhs, Complex rhs)
        {
            Complex result = new Complex();
            result.re = lhs.re * rhs.re;
            result.im = lhs.im * rhs.im;
            return result;
        }
    }
}


运行测试

11111.jpg


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

标签: Algorithms

Powered by emlog  蜀ICP备18021003号   sitemap

川公网安备 51019002001593号