数字信号处理C语言（3） ------FFT

快速傅里叶变换

fft.c

#include"math.h"

void fft(double *x,double *y,int n,int sign)
{
    int i,j,k,l,m,n1,n2;
    double c,c1,e,s,s1,t,tr,ti;
    for(j=1,i=1;i<16;i++)
    {
        m=i;
        j=2*j;
        if(j==n)break;
    }
    n1=n-1;
    for(j=0,i=0;i<n1;i++)
    {
        if(i<j)
        {
            tr=x[j];
            ti=y[j];
            x[j]=x[i];
            y[j]=y[i];
            x[i]=tr;
            y[i]=ti;
        }
        k=n/2;
        while(k<(j+1))
        {
            j=j-k;
            k=k/2;
        }
        j=j+k;
    }

    n1=1;
    for(l=1;l<=m;l++)
    {
        n1=2*n1;
        n2=n1/2;
        e=3.14159265359/n2;
        c=1.0;
        s=0.0;
        c1=cos(e);
        s1=-sign*sin(e);
        for(j=0;j<n2;j++)
        {
            for(i=j;i<n;i+=n1)
            {
                k=i+n2;
                tr=c*x[k]-s*y[k];
                ti=c*y[k]+s*x[k];
                x[k]=x[i]-tr;
                y[k]=y[i]-ti;
                x[i]=x[i]+tr;
                y[i]=y[i]+ti;
            }
            t=c;
            c=c*c1-s*s1;
            s=t*s1+s*c1;
        }
    }
    if(sign==-1)
    {
        for(i=0;i<n;i++)
        {
            x[i]/=n;
            y[i]/=n;
        }
    }
}

#include <QCoreApplication>
#include "fft.c"
int main(int argc, char *argv[])
{
    QCoreApplication d(argc, argv);
    int i,j,n;
    double a1,a2,c,c1,c2,d1,d2,q1,q2,w,w1,w2;
    double x[32],y[32],a[32],b[32];
    n=32;
    a1=0.9;
    a2=0.3;
    x[0]=1.0;
    y[0]=0.0;
    for(i=1;i<n;i++)
    {
        x[i]=a1*x[i-1]-a2*y[i-1];
        y[i]=a2*x[i-1]+a1*y[i-1];
    }
    printf("
 INPUT
");
    for(i=0;i<n/2;i++)
    {
        for(j=0;j<2;j++)
        {
            printf("%10.7f+J %10.7f",x[2*i+j],y[2*i+j]);
        }
        printf("
");
    }
    q1=x[n-1];
    q2=y[n-1];
    for(i=0;i<n;i++)
    {
        w=6.28348530718/n*i;
        w1=cos(w);
        w2=-sin(w);
        c1=1.0-a1*w1+a2*w2;
        c2=a1*w2+a2*w1;
        c=c1*c1+c2*c2;
        d1=1.0-a1*q1+a2*q2;
        d2=a1*q2+a2*q1;
        a[i]=(c1*d1+c2*d2)/c;
        b[i]=(c2*d1-c1*d2)/c;
    }
    printf("
 Theoretical DFT
");
    for(i=0;i<n/2;i++)
    {
        for(j=0;j<2;j++)
        {
            printf("%10.7f+J %10.7f",a[2*i+j],b[2*i+j]);
        }
        printf("
");
    }
    fft(x,y,n,1);
    printf("
DFT
");
    for(i=0;i<n/2;i++)
    {
        for(j=0;j<2;j++)
        {
            printf("%10.7f+J %10.7f",x[2*i+j],y[2*i+j]);
        }
        printf("
");
    }
    fft(x,y,n,-1);
    printf("
iDFT
");
    for(i=0;i<n/2;i++)
    {
        for(j=0;j<2;j++)
        {
            printf("%10.7f+J %10.7f",x[2*i+j],y[2*i+j]);
        }
        printf("
");
    }
    return d.exec();
}