#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[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 prntft(double data[], double nn2)
{
	int n,m,mm;
    cout<<"n       Real H(n)     Imag H(n)     Real H(N-n)    Imag H(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<<setprecision(0)<<setiosflags(ios::fixed)<<setw(1)<<m;
		cout<<setprecision(6)<<setiosflags(ios::fixed)<<setw(14)<<data[n];
		cout<<setprecision(6)<<setiosflags(ios::fixed)<<setw(14)<<data[n+1];
		cout<<setprecision(6)<<setiosflags(ios::fixed)<<setw(14)<<data[mm];
		cout<<setprecision(6)<<setiosflags(ios::fixed)<<setw(14)<<data[mm+1]<<endl;
    }
}

void main()
{
    //program d12r1
    //driver for routine four1
	int nn,nn2,i,isign;
    double data[65], dcmp[65],aaa,j;
    nn = 32;
    nn2 = 2 * nn;   
    cout<<setw(1)<<"h(t)=real-valued even-function"<<endl;
    cout<<setw(1)<<"H(n)=H(N-n) and real?"<<endl;
    for (i = 1; i<=2 * nn - 1; i+=2)
	{
        data[i] = 1.0 / (((i - nn - 1.0) / nn)*((i - nn - 1.0) / nn)  + 1.0);
        data[i + 1] = 0.0;
    }
    isign = 1;
    four1(data, nn, isign);
    prntft(data, nn2);
    cout<<setw(1)<<"h(t)=imagnary-valued even-function"<<endl;
    cout<<setw(1)<<"H(n)=H(N-n) and imaginary?"<<endl;
    for (i = 1; i<= 2 * nn - 1; i+=2)
	{
        data[i + 1] = 1.0 / (((i - nn - 1.0) / nn)*((i - nn - 1.0) / nn) + 1.0);
        data[i] = 0.0;
    }
    isign = 1;
    four1(data, nn, isign);
    prntft(data, nn2);
	cout<<setw(1)<<"h(t)=real-valued odd-function"<<endl;
    cout<<setw(1)<<"H(n)=-H(N-n) and imaginary?"<<endl;
    for (i = 1; i<=2 * nn - 1; i+=2)
	{
        data[i] = (i - nn - 1.0) / nn / (((i - nn - 1.0) / nn)*((i - nn - 1.0) / nn) + 1.0);
        data[i + 1] = 0.0;
    }
    data[1] = 0.0;
    isign = 1;
    four1(data, nn, isign);
    prntft(data, nn2);
	cout<<setw(1)<<"h(t)=imagnary-valued odd-function"<<endl;
    cout<<setw(1)<<"H(n)=-H(N-n) and real?"<<endl;
    for( i = 1; i<= 2 * nn - 1; i+=2)
	{
        aaa = (i - nn - 1.0);
        data[i + 1] = aaa / nn / (((i - nn - 1.0) / nn)*((i - nn - 1.0) / nn) + 1.0);
        data[i] = 0.0;
    }
    data[2] = 0.0;
    isign = 1;
    four1(data, nn, isign);
    prntft(data, nn2);
    //transrorm, inverse-transform test
    for (i = 1; i<=2 * nn - 1; i+=2)
	{
        data[i] = 1.0 / ((0.5 * (i - nn - 1) / nn)*(0.5 * (i - nn - 1) / nn) + 1.0);
        dcmp[i] = data[i];
        data[i + 1] = (0.25 * (i - nn - 1) / nn);
        data[i + 1] = data[i + 1] * exp(-(0.5 * (i - nn - 1) / nn)*(0.5 * (i - nn - 1) / nn));
        dcmp[i + 1] = data[i + 1];
    }
    isign = 1;
    four1(data, nn, isign);
    isign = -1;
    four1(data, nn, isign);
	cout<<endl;
	cout<<"     Double Fourier transform:        Original data:"<<endl;
    cout<<"k       Real h(k)     Imag h(k)     Real h(k)    Imag h(k)"<<endl;
    for (i = 1; i<=nn; i+=2)
	{
        j = (i + 1) / 2;
        cout<<setprecision(0)<<setw(1)<<j;
        cout<<setprecision(6)<<setw(14)<<dcmp[i];
        cout<<setprecision(6)<<setw(14)<<dcmp[i + 1];
        cout<<setprecision(6)<<setw(14)<<data[i]/nn;
        cout<<setprecision(6)<<setw(14)<<data[i + 1] / nn<<endl;
    }
}