FFT C语言 修改了matlab

FFT修改

FFT-dit的算法讨论

级的概念

$M=lo{g}_{2}N$,所以N点DFT可分成M级。

下图为8点FFT时间抽取算法信号流程图，方框分别为m=0,m=1,m=2;

码位倒置

交换后输出序列x(k)依照正序排列，但输入序列x(n)的次序不再是原来的自然次序，这是由于将x(n)按奇,偶分开产生的。

雷德算法： 对于自然顺序（二进制）我们是在低位加 1 得到下一位数，对于倒位序我们是在高位加 1 向低位进位。比如已知一个倒位序数是J求其下一个倒位序数，N位总数 ，把J与N/2比较若J<N/2则J的最高位为 0 ，把最高位置 1 ，就得到了J的下一个倒位序数；若J>=N/2则说明J的最高位为1 ，把最高位置0 ，比较次高位，若次高位为0 ，则把次高位置1，得到J的下一个倒位序，若次高位为1 ， 则把次高位置0,以此类推…

参考：https://blog.csdn.net/corcplusplusorjava/article/details/17119567

Rader参考程序：

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
int main(void)
{
    int array[8]={0,1,2,3,4,5,6,7};
    int i,j,k;
    int N = 8;
    int temp;
    j = 0;
    for(i = 0; i < N -1; i ++)
    {
        if(i < j)
        {
            temp = array[i];
            array[i] = array[j];
            array[j] = temp;
        }
        k = N >> 1;
        while( k <= j)
        {
            j = j - k;
            k >>= 1;
        }

        j = j + k;
    }
    for( i = 0; i < N; i ++)
        printf("%d %d\n",i,array[i]);
    printf("\n");
    return 0;
}

结果：

$W^r$因子的分布

规律：第m级，$W_{2^{m+1}}^r$

在fft中，乘法主要来自旋转因子。因为$W^r=cos(\frac{2\pi r}{N})-jsin(\frac{2\pi r}{N}),所以在对W^r$相乘时，必须产生相应的正，余弦函数。编程时，产生正，余弦函数两个方法，一个是在每一步直接产生，二是在程序开始前预先计算W^r存于一个数组中，等效一个正，余弦函数的“表”。

蝶形单元

第m级，有

$X_{m+1}(p)=x_m(p)+W_N^r{X}_m(q)$

$X_{m+1}(q)=x_m(p)-W_N^r{X}_m(q)$

"组"的概念

每一级的N/2个蝶形单元可以分成若干组，每一组有着相同的结构及W^r因子分布。

第m级的组数是$N/2^{m+1}$.

程序

A.1<<i是把1左移i位，每次左移以为就是乘以2，所以1<<i的结果是1乘以2的i次方 i<<1就是把i左移一位，即i乘以2，假如i=5,最后结果就是5*2=10

B.FFT蝶形算法使用三重循环，下一层的数据都是由上一层计算得到的。

#include <stdio.h>
#include<stdlib.h>
#define pi 3.1415
#define N  64
typedef struct
{
    double real,imag;
} complex; //复数结构
complex fft_outa[100];//a的全部数据
complex fft_outao[100];//a的单个数据
complex fft_outb[100];//b的全部数据
complex fft_outbo[100];//b的单个数据
double x_r[N],x_i[N];
double y_r[N],y_[N];

void Rader(double x_r[N],double x_i[N])//雷德算法排序
{
    int i,j,k;
    j=0;
    double t_r,t_i;
    for(i=0; i<N-1; i++)
    {
        if(i<j)
        {
            t_r=x_r[i];
            t_i=x_i[i];
            x_r[i]=x_r[j];
            x_i[i]=x_i[j];
            x_r[j]=t_r;
            x_i[j]=t_i;
        }
        k=N>>1;
        while(k<=j)
        {
            j=j-k;
            k>>=1;
        }
        j=j+k;
    }
}
void dit(double x_r[N],double x_i[N])//三层循环
{
    int m,r,p,q;
    int step,factor_step;
    int a,b,i;
    double t_real,t_imag;
    double w_re,w_im;
    for(m=0; m<log2(N); m++) //第一层循环
    {
        a=1;
        for(i=0;i<m+1;i++)
            a=a*2;
        b=a/2;
        w_re=cos(pi/b);
        w_im=-sin(pi/b);
       step=1<<(m+1);
       // factor_step=N>>(m+1);//旋转因素变化速度
        //初始化旋转因子
        double factor_real=1.0;
        double factor_imag=0.0;
        for(r=0; r<step/2; r++)
        {
           //   w_re=cos(2*pi*(r+1)*factor_step/N);
           //   w_im=-sin(2*pi*(r+1)*factor_step/N);
            for(p=r; p<N; p+=step)
            {
                q=p+step/2;//蝶形运算的两个输入
                t_real=x_r[q]*factor_real-x_i[q]*factor_imag;
                t_imag=x_r[q]*factor_imag+x_i[q]*factor_real;
                x_r[q]=x_r[p]-t_real;
                x_i[q]=x_i[p]-t_imag;
                x_r[p]=x_r[p]+t_real;
                x_i[p]=x_i[p]+t_imag;
            }

            t_real=factor_real*w_re-factor_imag*w_im;
            t_imag=factor_real*w_im+factor_imag*w_re;
            factor_real=t_real;
            factor_imag=t_imag;

        }
    }

}

int main()
{
    int n,i;
    double y[N];
    for(n=0; n<N; n++)//输入信号
    {
        x_r[n]=cos(n*pi/6);
        x_i[n]=0;
    }
    Rader(x_r,x_i);
    dit(x_r,x_i);
    for(i=0; i<N; i++)//输出fft
    {
        y[i]=sqrt(x_r[i]*x_r[i]+x_i[i]*x_i[i]);
        printf("%d %f\n",i,y[i]);
    }
}

c结果：

matlab的程序结果：

N=64;
n=0:1:N;
x=cos(n*pi/6);
X=fft(x);
figure
stem(n,abs(X))

对比结果还是有差异。

FFT过程中，对C理解还不够，对蝶形运算算法还不理解，对三重循环还不熟悉 。

修改了，matlab的程序

>> n=0:1:63;
>> x=cos(n*pi/6);
>> y=fft(x);
>> figure
>> stem(n,abs(y))

得出和C一样的结果