数字信号处理C语言——离散傅里叶变换DFT/离散傅里叶反变换IDFT

    #include "math.h"
    double dft(a,b,A,B,N)
    int N;
    double a[],b[],A[],B[];
    {
    double PI=acos(-1.0);
    int n,k;
    double c,d,w0,w,s;
    w0=2*PI/N;
    for(k=0;k<N;k++){
        w=k*w0;
        A[k]=B[k]=0.0;
        for(n=0;n<N;n++){
            d=n*w;
            c=cos(d);
            s=sin(d);
            A[k]+=c*a[n]+s*b[n];
            B[k]+=c*b[n]-s*a[n];
        }
    }
    return PI;
    }

    #include "math.h"
    void idft(A,B,a,b,N)
    int N;
    double A[],B[],a[],b[];
    {
    double PI=acos(-1.0);
    int n,k;
    double c,d,w0,w,s;
    w0=2*PI/N;
    for(n=0;n<N;n++){
        w=n*w0;
        a[n]=b[n]=0.0;
        for(k=0;k<N;k++){
            d=k*w;
            c=cos(d);
            s=sin(d);
            a[n]+=(c*A[k]-s*B[k])/N;
            b[n]+=(c*B[k]+s*A[k])/N;
        }
    }
    return PI;
    }

    #include "stdio.h"
    #include "math.h"
    #include "dft.c"
    #include "idft.c"
    main(){
    int i,j,N;
    printf("N=");
    scanf("%d",&N);
    double cNR,cNI;//complexNumReal,complexNumImaginary
    double a[N],b[N],A[N],B[N];
    printf("complexNumReal=");
    scanf("%lf",&cNR);
    printf("complexNumImaginary=");
    scanf("%lf",&cNI);
    a[0]=1.0;
    b[0]=0.0;
    for(i=1;i<N;i++){
        a[i]=cNR*a[i-1]-cNI*b[i-1];
        b[i]=cNI*a[i-1]+cNR*b[i-1];
    }
    printf("----Original Sequence----\n");
    for(i=0;i<N/2;i++){
        for(j=0;j<2;j++)
            printf("|%d:\t%10.6f + j%10.6f|\t",2*i+j+1,a[2*i+j],b[2*i+j]);
        printf("\n");
    }
    
    double PI=dft(a,b,A,B,N);
    printf("----Discrete Fourier Transform----\n");
    for(i=0;i<N/2;i++){
        for(j=0;j<2;j++)
            printf("|%d:\t%10.6f + j%10.6f |\t",2*i+j+1,A[2*i+j],B[2*i+j]);
        printf("\n");
    }
    
    double nR,nI;//numeratorR,numeratorI
    double dR,dI;//denominatorR,denominatorI
    double w,WkR,WkI;
    printf("----Theoretical Discrete Fourier Transform----\n");
    for(int k=0;k<N;k++){
        w=(2*PI/N)*k;   //W^k=e^-j*k*(2*PI/N)=cos(w)+j*(-1)*sin(w)
        WkR=cos(w);
        WkI=(-1)*sin(w);
        nR=1-cNR*a[N-1]+cNI*b[N-1];
        nI=cNI*a[N-1]+cNR*b[N-1];
        dR=1-cNR*WkR+cNI*WkI;
        dI=cNI*WkR+cNR*WkI;
        A[k]=(nR*dR+nI*dI)/(pow(dR,2)+pow(dI,2));
        B[k]=(nR*dI-nI*dR)/(pow(dR,2)+pow(dI,2));
    }
    for(i=0;i<N/2;i++){
        for(j=0;j<2;j++)
            printf("|%d:\t%10.6f + j%10.6f |\t",2*i+j+1,A[2*i+j],B[2*i+j]);
        printf("\n");
    }
    
    idft(A,B,a,b,N);
    printf("****Inverse Discrete Fourier Transform****\n");
    for(i=0;i<N/2;i++){
        for(j=0;j<2;j++)
            printf("|%d:\t%10.6f + j%10.6f |\t",2*i+j+1,a[2*i+j],b[2*i+j]);
        printf("\n");
    }
    }