fft.cpp

几个FFT算法

#include <iostream.h>
#include <math.h>
#include <iomanip.h>
#include <stdlib.h>
#include <fstream.h>
#include <string>
#include <process.h>
#include <stdio.h>
 
void four1(double data[], int nn, int isign)
{
    int n,j,i,m,mmax,istep;
    double tempr,tempi,theta,wpr,wpi,wr,wi,wtemp;
    n = 2 * nn;
    j = 1;
    for (i = 1; i<=n; i=i+2)
    {
        if(j > i)
        {
            tempr = data[j];
            tempi = data[j + 1];
            data[j] = data[i];
            data[j + 1] = data[i + 1];
            data[i] = tempr;
            data[i + 1] = tempi;
        }
        m = n / 2;
        while (m >= 2 && j > m)
        {
            j = j - m;
            m = m / 2;
        }
        j = j + m;
    }
    mmax = 2;
    while(n > mmax)
    {
        istep = 2 * mmax;
        theta = 6.28318530717959 / (isign * mmax);
        wpr = -2.0 * sin(0.5 * theta)*sin(0.5 * theta);
        wpi = sin(theta);
        wr = 1.0;
        wi = 0.0;
        for(m = 1; m<=mmax; m=m+2)
        {
            for (i = m; i<=n; i=i+istep)
            {
                j = i + mmax;
                tempr=double(wr)*data[j]-double(wi)*data[j+1];
                tempi=double(wr)*data[j+1]+double(wi)*data[j];
                data[j] = data[i] - tempr;
                data[j + 1] = data[i + 1] - tempi;
                data[i] = data[i] + tempr;
                data[i + 1] = data[i + 1] + tempi;
            }
            wtemp = wr;
            wr = wr * wpr - wi * wpi + wr;
            wi = wi * wpr + wtemp * wpi + wi;
        }
        mmax = istep;
    }
}
 
/*
void twofft(double data1[], double data2[], double fft1[], double fft2[], int& n)
{
    int j,n2,j2;
    double c1r,c1i,c2r,c2i,conjr,conji,h1r,h1i,h2r,h2i;
    c1r = 0.5;
    c1i = 0.0;
    c2r = 0.0;
    c2i = -0.5;
    for (j = 1; j<=n; j++)
    {
        fft1[2 * j - 1] = data1[j];
        fft1[2 * j] = data2[j];
    }
    four1(fft1, n, 1);
    fft2[1] = fft1[2];
    fft2[2] = 0.0;
    fft1[2] = 0.0;
    n2 = 2 * (n + 2);
    for (j = 2; j<=n / 2 + 1; j++)
    {
        j2 = 2 * j;
        conjr = fft1[n2 - j2 - 1];
        conji = -fft1[n2 - j2];
        h1r=c1r*(fft1[j2-1]+conjr)-c1i*(fft1[j2]+conji);
        h1i=c1i*(fft1[j2-1]+conjr)+c1r*(fft1[j2]+conji);
        h2r=c2r*(fft1[j2-1]-conjr)-c2i*(fft1[j2]-conji);
        h2i=c2i*(fft1[j2-1]-conjr)+c2r*(fft1[j2]-conji);
        fft1[j2 - 1] = h1r;
        fft1[j2] = h1i;
        fft1[n2 - j2 - 1] = h1r;
        fft1[n2 - j2] = -h1i;
        fft2[j2 - 1] = h2r;
        fft2[j2] = h2i;
        fft2[n2 - j2 - 1] = h2r;
        fft2[n2 - j2] = -h2i;
    }
}
 

int cint(double x)
{
    int temp;
    double iprt;
    if (x>0)
    {
        x=modf(x,&iprt);
        if(fabs(x)<0.5)
            temp=int(iprt);
        else
            temp=int(iprt+1);
    }
    else if(x==0)
        temp=0;
    else
    {
        x=modf(x,&iprt);
        if(fabs(x)<0.5)
            temp=int(iprt);
        else
            temp=int(iprt)-1;
    }
    return temp;
}
 
void prntft(double data[], double nn2)
{
    int n,m,mm;
    cout<<"n   real(n)   imag.(n)   real(N-n)   imag.(N-n)"<<endl;
    cout<<setw(1)<<"0";
    cout<<setw(14)<<data[1];
    cout<<setw(14)<<data[2];
    cout<<setw(14)<<data[1];
    cout<<setw(14)<<data[2]<<endl;
    for (n = 3; n<=(nn2 / 2) + 1; n=n+2)
    {
        m = (n - 1) / 2;
        mm = nn2 + 2 - n;
        cout<<setiosflags(ios::fixed);
        cout<<setprecision(0)<<setw(1)<<m;
        cout<<setprecision(6)<<setw(14)<<data[n];
        cout<<setprecision(6)<<setw(14)<<data[n+1];
        cout<<setprecision(6)<<setw(14)<<data[mm];
        cout<<setprecision(6)<<setw(14)<<data[mm+1]<<endl;
    }
}
 
void main()
{
    //program d12r2
    //driver for routine twofft
    int n,i,n2,isign;
    double data1[33], data2[33], fft1[65], fft2[65],per,x;
    n = 32;
    n2 = 2 * n;
    per = 8.0;
    const double pi = 3.1415926;   
    for(i = 1; i<=n; i++)
    {
        x = 2.0 * pi * i / per;
        data1[i] = cint(cos(x));
        data2[i] = cint(sin(x));
    }
    twofft(data1, data2, fft1, fft2, n);
 
    cout<<setw(1)<<"Fourier transform of first function:"<<endl;
    prntft(fft1, n2);
    cout<<setw(1)<<"Fourier transform of second function:"<<endl;
    prntft(fft2, n2);
    //invert transform
    isign = -1;
    four1(fft1, n, isign);
    cout<<setw(1)<<"Inverted transform = first function:"<<endl;
    prntft(fft1, n2);
    four1(fft2, n, isign);
    cout<<setw(1)<<"Inverted transform = second function:"<<endl;
    prntft(fft2, n2);
}
*/