d12r2.cpp

傅立叶变换是信号处理时最常用的算法之一

#include <iostream.h>
#include <math.h>
#include <iomanip.h>
#include <stdlib.h>
#include <fstream.h>
#include <string>
#include <process.h>
#include<stdio.h>

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 four1(double data[65], 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;
    }
}

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);
}