fft.cpp

计算快速fouries变换的c++源程序

#include <iostream>
#include <cmath>
#include <iomanip>
#include <stdlib.h>
#include <fstream>
#include <string>
#include <process.h>
#include<stdio.h>
using namespace std;


void four1(double data[], int nn, int isign);
void realft(double data[], int &n, int &isign);
int main()
{	
	int i,np,npp2,n,temp,temp1;
    double data[35], size1[33],eps,width1,per,big,small,scal1,nlim;
	temp=1;
	temp1=-1;
    eps = 0.001;
    np = 32;
    npp2 = np + 2;
    width1 = 50.0;
    const double pi = 3.14159;
    n = np/2;
	per = 32;
/*
    cout<<"Period of sinusoid in channels (2-"<<np<<")"<<endl;
	cout<<per<<endl;
*/
    if (per <= 0.0)
		exit(1);
    for (i = 0; i<np; i++)
	{
        data[i] = cos(2.0 * pi * i / per);
    }
    realft(data, n, temp);
    big = -10000000000.0;
    for( i = 0; i<n; i++)
	{
        size1[i] = sqrt(data[2 * i + 1]*data[2 * i + 1] + data[2 * i]*data[2 * i]);
        if (i == 0)
		    size1[i] = data[i];
        if (size1[i] > big )
			big = size1[i];
    }
    scal1 = width1 / big;
/*
	int j;
    for (i = 0; i<n; i++)
	{
        nlim = scal1 * size1[i] + eps;
        cout<<setw(4)<<i<<"  ";
        for (j = 0; j<nlim + 1; j++)
		{
            cout<<"*";
		}
		cout<<endl;
    }
*/
    realft(data, n, temp1);
    big = -10000000000.0;
    small = 10000000000.0;
    for (i = 0;i<np;i++)
	{
        if (data[i] < small)
			small = data[i];
        if (data[i] > big)
			big = data[i];
    }
    scal1 = width1 / (big - small);
    for (i = 0; i<np; i++)
	{
        nlim = int(scal1 * (data[i] - small) + eps);
/*
        cout<<setw(4)<<i<<"  ";
        for (j = 1; j<=nlim + 1; j++)
		{
            cout<<"*";
        }
		cout<<endl;
*/
    }
	return 0;
}

void realft(double data[], int &n, int &isign)
{
	int n2p3,i,i1,i2,i3,i4;
	double theta,c1,c2,wpr,wpi,wr,wi,wrs,wis,h1r,h1i,h2r;
	double h2i,wtemp,wri;
	wri=0.0;
    theta = 6.28318530717959 / 2.0 / n;
    c1 = 0.5;
    if (isign == 1)
	{
        c2 = -0.5;
        four1(data, n, 1);
	}
    else
	{
        c2 = 0.5;
        theta = -theta;
	}
    wpr = -2.0 * sin(0.5 * theta) * sin(0.5 * theta);
    wpi = sin(theta);
    wr = 1.0 + wpr;
    wi = wpi;
    n2p3 = 2 * n + 1;
    for (i = 1; i<n / 2 + 1; i++)
	{
        i1 = 2 * i;
        i2 = i1 + 1;
        i3 = n2p3 - i2;
        i4 = i3 + 1;
        wrs = float(wr);
        wis = float(wi);
        h1r = c1 * (data[i1] + data[i3]);
        h1i = c1 * (data[i2] - data[i4]);
        h2r = -c2 * (data[i2] + data[i4]);
        h2i = c2 * (data[i1] - data[i3]);
        data[i1] = h1r + wrs * h2r - wis * h2i;
        data[i2] = h1i + wrs * h2i + wis * h2r;
        data[i3] = h1r - wrs * h2r + wis * h2i;
        data[i4] = -h1i + wrs * h2i + wis * h2r;
        wtemp = wr;
        wr = wr * wpr - wi * wri + wr;
        wi = wi * wpr + wtemp * wpi + wi;
    }
    if( isign == 1 )
	{
        h1r = data[0];
        data[0] = h1r + data[1];
        data[1] = h1r - data[1];
	}
    else
	{
        h1r = data[0];
        data[0] = c1 * (h1r + data[1]);
        data[1] = c1 * (h1r - data[1]);
        four1(data, n, -1);
    }
}


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-1];
			data[j-1] = data[i-1];
            data[i-1] = tempr;
			tempi = data[j];
            data[j] = data[i];
			data[i] = 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-1]-double(wi)*data[j];
                tempi=double(wr)*data[j]+double(wi)*data[j-1];
                data[j-1] = data[i-1] - tempr;
                data[j] = data[i] - tempi;
                data[i-1] = data[i-1] + tempr;
                data[i] = data[i] + tempi;
            }
            wtemp = wr;
            wr = wr * wpr - wi * wpi + wr;
            wi = wi * wpr + wtemp * wpi + wi;
        }
        mmax = istep;
    }
}