dft.cpp

DFT,FFT,IFFT算法

// dft.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include "math.h"
#include <stdio.h>
#include <conio.h>
void dft(double x[],double y[],double a[],double b[],int n,int sign);
void fft(double x[],double y[],int n,int sign);
int main(int argc, char* argv[])
{
	printf("Hello World!\n");
	int i,j,n;
	double a1,a2,c,c1,c2,d1,d2,q1,q2,w,w1,w2;
	double x[32],y[32],a[32],b[32];
	n=32;
	a1 =0.9;
	a2 =0.3;
	x[0] = 1.0;
	y[0]=0.0;
	for (i=1;i<n;i++) {
		x[i]=a1*x[i-1]-a2*y[i-1];
		y[i]=a2*x[i-1]+a1*y[i-1];

	}
	printf("\nOriginal Sequence \n");
	for (i=0;i<n/2;i++) {
		for (j=0;j<2;j++) {
			printf("    %10.7f  +j %10.7f",x[2*i+j],y[2*i+j]);
		}
			printf("\n");

	}
	q1=x[n-1];
	q2=y[n-1];

	dft(x,y,a,b,n,1);
	printf("\n Discrete Fourier Transform \n");
	for (i=0;i<n/2;i++) {
		for (j=0;j<2;j++) {
			printf("	%10.7f + J %10.7f",a[2*i+j],b[2*i+j]);
		}
		printf("\n");

	}
	for (i=0;i<n;i++) {

		w = 6.28318530718/n*i;
		w1 = cos(w);
		w2 = -sin(w);
		c1 = 1.0-a1*w1+a2*w2;
		c2 = a1*w2 +a2*w1;
		c = c1*c1 + c2*c2;
		d1 = 1.0-a1*q1+a2*q2;
		d2 = a1*q2+a2*q1;
		x[i] = (c1*d1+c2*d2)/c;
		y[i]=(c2*d1 - c1*d2)/c;
	}
	printf("\n Theoretical Discrete Fourier Transform \n");
	for (i=0;i<n/2;i++) {
		for (j = 0;j<2;j++) {
			printf("	%10.7f  +J %10.7f",x[2*i+j],y[2*i+j]);
		}
		printf("\n");
	}
	dft(a,b,x,y,n,-1);
	printf("\n Inverse Discrete Fourier Transform \n");
	for (i=0;i<n/2;i++) {
		for (j = 0;j<2;j++) {
			printf("	%10.7f  +J %10.7f",x[2*i+j],y[2*i+j]);
		}
		printf("\n");
	}

		n=32;
	a1 =0.9;
	a2 =0.3;
	x[0] = 1.0;
	y[0]=0.0;
	for (i=1;i<n;i++) {
		x[i]=a1*x[i-1]-a2*y[i-1];
		y[i]=a2*x[i-1]+a1*y[i-1];

	}
	printf("\nOriginal input Sequence \n");
	for (i=0;i<n/2;i++) {
		for (j=0;j<2;j++) {
			printf("    %10.7f  +j %10.7f",x[2*i+j],y[2*i+j]);
		}
			printf("\n");

	}
	q1=x[n-1];
	q2=y[n-1];
	for (i=0;i<n;i++) {

		w = 6.28318530718/n*i;
		w1 = cos(w);
		w2 = -sin(w);
		c1 = 1.0-a1*w1+a2*w2;
		c2 = a1*w2 +a2*w1;
		c = c1*c1 + c2*c2;
		d1 = 1.0-a1*q1+a2*q2;
		d2 = a1*q2+a2*q1;
		a[i] = (c1*d1+c2*d2)/c;
		b[i]=(c2*d1 - c1*d2)/c;
	}
	printf("\n Theoretical DFT \n");
	for (i=0;i<n/2;i++) {
		for (j = 0;j<2;j++) {
			printf("	%10.7f  +J %10.7f",a[2*i+j],b[2*i+j]);
		}
		printf("\n");
	}
	fft(x,y,n,1);
	printf("\n Fast Fourier Transform \n");
	for (i=0;i<n/2;i++) {
		for (j = 0;j<2;j++) {
			printf("	%10.7f  +J %10.7f",x[2*i+j],y[2*i+j]);
		}
		printf("\n");
	}
	fft(x,y,n,-1);
	printf("\n Inverse Fast Fourier Transform \n");
	for (i=0;i<n/2;i++) {
		for (j = 0;j<2;j++) {
			printf("	%10.7f  +J %10.7f",x[2*i+j],y[2*i+j]);
		}
		printf("\n");
	}


	getch();
	return 0;
}

void dft(double x[],double y[],double a[],double b[],int n,int sign)
{
	int i,k;
	double c,d,q,w,s;
	q = 6.28318530718/n;
	for (k=0;k<n;k++) {
		w = k*q;
		a[k]=b[k]=0.0;
		for(i=0;i<n;i++)
		{
			d=i*w;
			c=cos(d);
			s=sin(d)*sign;
			a[k]+=c*x[i]+s*y[i];
			b[k]+=c*y[i]-s*x[i];


		}

	}
	if(sign==-1)
	{
		c = 1.0/n;
		for (k=0;k<n;k++) {
			a[k]=c*a[k];
			b[k]=c*b[k];
		}
	}
}
void fft(double x[],double y[],int n,int sign)
{
	int i,j,k,l,m,n1,n2;
	double c,c1,e,s,s1,t,tr,ti;
	for (j=1,i=1;i<16;i++) {
		m=i;
		j=2*j;
		if(j==n)break;
	}
	n1 = n-1;
	for (j=0,i=0;i<n1;i++) {
		if (i<j) {
			tr = x[j];
			ti =  y[j];
			x[j] = x[i];
			y[j]=y[i];
			x[i]=tr;
			y[i]=ti;
		}
		k=n/2;
		while (k<(j+1)) {
			j = j-k;
			k = k/2;
		}
		j=j+k;

	}
	n1 =1;
	for(l=1;l<=m;l++)
	{
		n1 = 2*n1;
		n2 = n1/2;
		e = 3.14159265359/n2;
		c = 1.0;
		s = 0.0;
		c1 = cos(e);
		s1 =sign*sin(e);
		for (j=0;j<n2;j++) {
			for (i=j;i<n;i+=n1) {
				k = i+n2;
				tr=c*x[k]-s*y[k];
				ti=c*y[k]+s*x[k];
				x[k]=x[i]-tr;
				y[k]=y[i]-ti;
				x[i]=x[i]+tr;
				y[i]=y[i]+ti;
			}
			t=c;
			c = c*c1 - s*s1;
			s = t*s1 + s*c1;
		}
	}
	if(sign==-1)
	{
		for (i=0;i<n;i++) {
			y[i]/=n;
			x[i]/=n;
		}
	}

}
void rfft(double x[],int n)
{
	int i,j,k,m,i1,i2,i3,i4,n1,n2,n4;
	double a,e,cc,ss,xt,t1,t2;
	for (j=1,i=1;i<16;i++) {
		m=i;
		j=2*j;
		if(j==n)
			break;

	}
	n1 = n-1;
	for (j=0,i=0;i<n1;i++) {
		if(i<j)
		{
			xt=x[j];
			x[j]=x[i];
			x[i]=xt;
		}
		k = n/2;
		while (k<j+1) {
			j=j-k;
			k=k/2;
		}
		j=j+k;
	}
	for (i=0;i<n;i+=2) {
		xt = x[i];
		x[i]=xt+x[i+1];
		x[i+1]=xt - x[i+1];
	}
	n2 =1;
	for (k=2;k<=m;k++) {
		n4=n2;
		n2 = 2*n4;
		n1 = 2*n2;
		e = 6.28318530718/n1;
		for (i=0;i<n;i+=n1) {
			xt = x[i];
			x[i]=xt+x[i+n2];
			x[i+n2]=xt-x[i+n2];
			x[i+n2+n4]=-x[i+n2+n4];
			a=e;
			for (j=1;j<=(n4-1);j++) {
				i1= i+j;
				i2 = i-j+n2;
				i3 = i+j+n2;
				i4 = i-j+n1;
				cc = cos(a);
				ss =sin(a);
				a = a+e;
				t1 = cc*x[i3]+ss*x[i4];
				t2 = ss*x[i3]-cc*x[i4];
				x[i4]=x[i2]-t2;
				x[i3]=-x[i2]-t2;
				x[i2]=x[i1]-t1;
				x[i1]=x[i1]+t1;
			}
		}
	}
	
}