FFT-C语言

http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform#i-21

一、彻底理解傅里叶变换

　　快速傅里叶变换(Fast Fourier Transform)是离散傅里叶变换的一种快速算法，简称FFT，通过FFT可以将一个信号从时域变换到频域。

　　模拟信号经过A/D转换变为数字信号的过程称为采样。为保证采样后信号的频谱形状不失真，采样频率必须大于信号中最高频率成分的2倍，这称之为采样定理。

假设采样频率为fs，采样点数为N，那么FFT结果就是一个N点的复数，每一个点就对应着一个频率点，某一点n(n从1开始)表示的频率为：fn=(n-1)*fs/N。

举例说明：用1kHz的采样频率采样128点，则FFT结果的128个数据即对应的频率点分别是0，1k/128，2k/128，3k/128，…，127k/128 Hz。

这个频率点的幅值为：该点复数的模值除以N/2(n=1时是直流分量，其幅值是该点的模值除以N)。

二、傅里叶变换的C语言编程　　

1、码位倒序。

　　假设一个N点的输入序列，那么它的序号二进制数位数就是t=log2N.
　　码位倒序要解决两个问题：①将t位二进制数倒序;②将倒序后的两个存储单元进行交换。
　　如果输入序列的自然顺序号i用二进制数表示，例如若最大序号为15，即用4位就可表示n3n2n1n0，则其倒序后j对应的二进制数就是n0n1n2n3，那么怎样才能实现倒序呢?利用C语言的移位功能!

/*变址计算，将x(n)码位倒置*/

void Reverse(complex *IN_X, int n)
{
    int row = log(n) / log(2); //row为FFT的N个输入对应的最大列数
    complex temp;  //临时交换变量
    unsigned short i = 0, j = 0, k = 0;
    double t;
    for (i = 0; i < row; i++)  //从第0个序号到第N-1个序号
    {
        k = i;   //当前第i个序号
        j = 0;   //存储倒序后的序号，先初始化为0
        t = row;  //共移位t次
        while ((t--) > 0)    //利用按位与以及循环实现码位颠倒
        {
            j = j << 1;
            j |= (k & 1);    //j左移一位然后加上k的最低位
            k = k >> 1;      //k右移一位，次低位变为最低位
                             /*j为倒序，k为正序，
                               从k右向左的值依次移位，
                               j向左依次填入对应的k的右向位*/
        }
        if (j > i)   //如果倒序后大于原序数，就将两个存储单元进行交换(判断j>i是为了防止重复交换
        {
            temp     = IN_X[i];
            IN_X[i] = IN_X[j];
            IN_X[j] = temp;
        }
    }
}


//复数加法的定义
complex add(complex a, complex b)  
{
    complex c;
    c.real = a.real + b.real;
    c.img = a.img + b.img;
    return c;
}

//复数乘法的定义
complex mul(complex a, complex b)  
{
    complex c;
    c.real = a.real*b.real - a.img*b.img;
    c.img = a.real*b.img + a.img*b.real;
    return c;
}

//复数减法的定义
complex sub(complex a, complex b)  
{
    complex c;
    c.real = a.real - b.real;
    c.img = a.img - b.img;
    return c;
}

2、第二个要解决的问题就是蝶形运算

 

　　对N个输入行，m表示第m列蝶形。

　　①第1级(第1列)每个蝶形的两节点“距离”为1，第2级每个蝶形的两节点“距离”为2，第3级每个蝶形的两节点“距离”为4，第4级每个蝶形的两节点“距离”为8。由此推得，
　　第m级蝶形运算，每个蝶形的两节点“距离”L=2^(m-1)。
　　②对于16点的FFT，第1级有16组蝶形，每组有1个蝶形;第2级有4组蝶形，每组有2个蝶形;第3级有2组蝶形，每组有4个蝶形;第4级有1组蝶形，每组有8个蝶形。由此可推出，
　　对于N点的FFT，第m级有N/(2L)组蝶形，每组有L=2^(m-1)个蝶形。
　　③旋转因子的确定
　　以16点FFT为例，第m级第k个旋转因子为，其中k=0～2^(m-1)-1，即第m级共有2^(m-1)个旋转因子，根据旋转因子的可约性，，所以第m级第k个旋转因子为，其中k=0～2^(m-1)-1。

　　生成一个的序列，代码如下：

/*产生一个周期欧拉形式的Wn的值*/
void InitWn(complex *w,int n)   //n为输入的个数,w为权值Wn
{
    int k;
    for (k = 0; k < n; k++)
    {
        w[k].real = cos(2 * PI / n * k);   //用欧拉公式计算旋转因子
        w[k].img = -1 * sin(2 * PI / n * k);
    }
}

3、算法实现

　　我们已经知道，N点FFT从左到右共有log2N级蝶形，每级有N/2L组，每组有L个。所以FFT的C语言编程只需用3层循环即可实现：最外层循环完成每一级的蝶形运算(整个FFT共log2N级)，中间层循环完成每一组的蝶形运算(每一级有N/2L组)，最内层循环完成单独1个蝶形运算(每一组有L个)。
　

/*快速傅里叶变换*/
void fft(complex *IN_X, int n)
{
    int row = log(n) / log(2);    //row为FFT的N个输入对应的最大列数
    int i = 0, j = 0, k = 0, L = 0;   
    complex top, bottom, product;
    Reverse(IN_X,n);  //调用变址函数
    for (i = 0; i < row ; i++)        /*第i级蝶形运算 */
    {
        L = 1 << i;  //L等于2的i次方,表示每一级蝶形中组间的距离，即一组中蝶形个数
        for (j = 0; j < n; j += 2 * L)    /*j为组偏移地址，每L个蝶形是一组，每级有N/2L组*/
        {
            for (k = 0; k < L; k++)        /*k为一组中蝶形的偏移地址，一个蝶形运算 每个group内的蝶形运算*/
            {
                product = mul(IN_X[j + k + L], Wn[n*k/2/L]);
                top = add(IN_X[j + k], product);
                bottom = sub(IN_X[j + k], product);
                IN_X[j + k] = top;
                IN_X[j + k + L] = bottom;
            }
        }
    }
}

4.全部代码

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <malloc.h>
#define N 1000
#define PI atan(1)* 4

/*定义复数类型*/
typedef struct {
    double real;
    double img;
}complex;

complex x[N], *Wn;                 /*输入序列,复数系数Wn*/
int size_x;                         /*size_x为采样的信号个数，必须为2的倍数*/
void fft(complex *IN_X, int n);         /*对输入IN_X，快速傅里叶变换*/
void InitWn(complex *w, int n);  /*生成一个长度为n的Wn欧拉形式序列*/
void Reverse(complex *IN_X, int n); /*对输入IN_X地址*/
complex add(complex, complex); /*复数加法*/
complex mul(complex, complex); /*复数乘法*/
complex sub(complex, complex); /*复数减法*/
void output();                   /*输出快速傅里叶变换的结果*/


int main()
{
    int i;                             /*输出结果*/
    system("cls"); //调用系统命令，CLS的功能为清除屏幕输出
    printf("输出DIT方法实现的FFT结果
");
    printf("Please input the size of x:
");//输入序列的大小
    scanf_s("%d", &size_x);
    printf("Please input the data in x[N]:
");//输入序列的实部和虚部
    for (i = 0; i < size_x; i++)
    {
        printf("请输入第%d个序列：", i);
        scanf_s("%lf%lf", &x[i].real, &x[i].img);
    }
    printf("输出倒序后的序列
");
    Wn = (complex *)malloc(sizeof(complex) * size_x);  //分配变换后的值的内存空间
    InitWn(Wn , size_x);//调用变换核
    fft(x, size_x);//调用快速傅里叶变换
    printf("输出FFT后的结果
");
    output();//调用输出傅里叶变换结果函数
    scanf_s("%d", &size_x);
    //free(Wn);
    return 0;
}


/*快速傅里叶变换*/
void fft(complex *IN_X, int n)
{
    int row = log(n) / log(2);    //row为FFT的N个输入对应的最大列数
    int i = 0, j = 0, k = 0, L = 0;   
    complex top, bottom, product;
    Reverse(IN_X,n);  //调用变址函数
    for (i = 0; i < row ; i++)        /*第i级蝶形运算 */
    {
        L = 1 << i;  //L等于2的i次方,表示每一级蝶形中组间的距离，即一组中蝶形个数
        for (j = 0; j < n; j += 2 * L)    /*j为组偏移地址，每L个蝶形是一组，每级有N/2L组*/
        {
            for (k = 0; k < L; k++)        /*k为一组中蝶形的偏移地址，一个蝶形运算 每个group内的蝶形运算*/
            {
                product = mul(IN_X[j + k + L], Wn[n*k/2/L]);
                top = add(IN_X[j + k], product);
                bottom = sub(IN_X[j + k], product);
                IN_X[j + k] = top;
                IN_X[j + k + L] = bottom;
            }
        }
    }
}


/*产生一个周期欧拉形式的Wn的值*/
void InitWn(complex *w,int n)   //n为输入的个数,w为权值Wn
{
    int k;
    for (k = 0; k < n; k++)
    {
        w[k].real = cos(2 * PI / n * k);   //用欧拉公式计算旋转因子
        w[k].img = -1 * sin(2 * PI / n * k);
    }
}


/*变址计算，将x(n)码位倒置*/
/*
码位倒序要解决两个问题：
    ①将t位二进制数倒序;②将倒序后的两个存储单元进行交换。
    如果输入序列的自然顺序号i用二进制数表示，
    例如若最大序号为15，即用4位就可表示n3n2n1n0，
    则其倒序后j对应的二进制数就是n0n1n2n3
*/
void Reverse(complex *IN_X, int n)
{
    int row = log(n) / log(2); //row为FFT的N个输入对应的最大列数
    complex temp;  //临时交换变量
    unsigned short i = 0, j = 0, k = 0;
    double t;
    for (i = 0; i < row; i++)  //从第0个序号到第N-1个序号
    {
        k = i;   //当前第i个序号
        j = 0;   //存储倒序后的序号，先初始化为0
        t = row;  //共移位t次
        while ((t--) > 0)    //利用按位与以及循环实现码位颠倒
        {
            j = j << 1;
            j |= (k & 1);    //j左移一位然后加上k的最低位
            k = k >> 1;      //k右移一位，次低位变为最低位
                             /*j为倒序，k为正序，
                               从k右向左的值依次移位，
                               j向左依次填入对应的k的右向位*/
        }
        if (j > i)   //如果倒序后大于原序数，就将两个存储单元进行交换(判断j>i是为了防止重复交换
        {
            temp = IN_X[i];
            IN_X[i] = IN_X[j];
            IN_X[j] = temp;
        }
    }
}



/*输出傅里叶变换的结果*/
void output()
{
    int i;
    printf("The result are as follows：
");
    for (i = 0; i < size_x; i++)
    {
        printf("%.4f", x[i].real);
        if (x[i].img >= 0.0001)printf("+%.4fj
", x[i].img);
        else if (fabs(x[i].img) < 0.0001)printf("
");
        else printf("%.4fj
", x[i].img);
    }
}

//复数加法的定义
complex add(complex a, complex b)
{
    complex c;
    c.real = a.real + b.real;
    c.img = a.img + b.img;
    return c;
}

//复数乘法的定义
complex mul(complex a, complex b)
{
    complex c;
    c.real = a.real*b.real - a.img*b.img;
    c.img = a.real*b.img + a.img*b.real;
    return c;
}

//复数减法的定义
complex sub(complex a, complex b)
{
    complex c;
    c.real = a.real - b.real;
    c.img = a.img - b.img;
    return c;
}