4fft.cpp

基于VC的基4FFT 变换 ,蝶形表单独生成

#include <iostream.h>
#include <math.h>
#define N 1024
#define pi 3.14159
void main(){
	double mod[N],source[N],im[3*N/4],re[3*N/4],tempre2[N],tempim2[N],tempre1[N],tempim1[N],temp;
	int i,j;
	for(i=0;i<N;i++){
		source[i]=4*sin(2*pi*100*i/1024+0.33);
		tempre1[i]=source[i];
		tempim1[i]=0;
	}
	//蝶形表
	for(i=0;i<(N*3/4);i++){
		im[i]=-sin(2*pi*i/N);
		re[i]=cos(2*pi*i/N);
	}
	int j1,j11,m,n,k=1,ii=0;
	int flag=1;
	n=0;
	for(i=N/4;i>=1;i/=4){
		ii++;
			if(flag==1){
				flag=2;
			    j1=0;
				for(m=0;m<k;m++){
					j11=0;
				for(j=m*i*4;j<i+m*i*4;j++){
					tempre2[j1]=tempre1[j]+tempre1[j+i]+tempre1[j+2*i]+tempre1[j+3*i];
					tempim2[j1]=tempim1[j]+tempim1[j+i]+tempim1[j+2*i]+tempim1[j+3*i];
					j1+=i;
					tempre2[j1]=tempre1[j]+tempim1[j+i]-tempre1[j+2*i]-tempim1[j+3*i];
					tempim2[j1]=tempim1[j]-tempre1[j+i]-tempim1[j+2*i]+tempre1[j+3*i];
					//用蝶形表re(j11*N/(i*4))im// 查表时控制步长即可
					temp=tempre2[j1]*re[j11*N/(i*4)]-tempim2[j1]*im[j11*N/(i*4)];
					tempim2[j1]=tempim2[j1]*re[j11*N/(i*4)]+tempre2[j1]*im[j11*N/(i*4)];
					tempre2[j1]=temp; 
					j1+=i;
					tempre2[j1]=tempre1[j]-tempre1[j+i]+tempre1[j+2*i]-tempre1[j+3*i];
					tempim2[j1]=tempim1[j]-tempim1[j+i]+tempim1[j+2*i]-tempim1[j+3*i];
					temp=tempre2[j1]*re[2*j11*N/(i*4)]-tempim2[j1]*im[2*j11*N/(i*4)];
					tempim2[j1]=tempim2[j1]*re[2*j11*N/(i*4)]+tempre2[j1]*im[2*j11*N/(i*4)];
					tempre2[j1]=temp;
					j1+=i;
					tempre2[j1]=tempre1[j]-tempim1[j+i]-tempre1[j+2*i]+tempim1[j+3*i];
					tempim2[j1]=tempim1[j]+tempre1[j+i]-tempim1[j+2*i]-tempre1[j+3*i];
					temp=tempre2[j1]*re[3*j11*N/(i*4)]-tempim2[j1]*im[3*j11*N/(i*4)];
					tempim2[j1]=tempim2[j1]*re[3*j11*N/(i*4)]+tempre2[j1]*im[3*j11*N/(i*4)];
					tempre2[j1]=temp;
					j1=j1-3*i+1;
					j11++;
				}
				j1=j1+3*i;
				}
			}else{
				flag=1;
			    j1=0;
				for(m=0;m<k;m++){
					j11=0;
				for(j=m*i*4;j<i+m*i*4;j++){
					tempre1[j1]=tempre2[j]+tempre2[j+i]+tempre2[j+2*i]+tempre2[j+3*i];
					tempim1[j1]=tempim2[j]+tempim2[j+i]+tempim2[j+2*i]+tempim2[j+3*i];
					j1+=i;
					tempre1[j1]=tempre2[j]+tempim2[j+i]-tempre2[j+2*i]-tempim2[j+3*i];
					tempim1[j1]=tempim2[j]-tempre2[j+i]-tempim2[j+2*i]+tempre2[j+3*i];
					//用蝶形表
					temp=tempre1[j1]*re[j11*N/(i*4)]-tempim1[j1]*im[j11*N/(i*4)];
					tempim1[j1]=tempim1[j1]*re[j11*N/(i*4)]+tempre1[j1]*im[j11*N/(i*4)];
					tempre1[j1]=temp;
					j1+=i;
					tempre1[j1]=tempre2[j]-tempre2[j+i]+tempre2[j+2*i]-tempre2[j+3*i];
					tempim1[j1]=tempim2[j]-tempim2[j+i]+tempim2[j+2*i]-tempim2[j+3*i];
					temp=tempre1[j1]*re[2*j11*N/(i*4)]-tempim1[j1]*im[2*j11*N/(i*4)];
					tempim1[j1]=tempim1[j1]*re[2*j11*N/(i*4)]+tempre1[j1]*im[2*j11*N/(i*4)];
					tempre1[j1]=temp;
					j1+=i;
					tempre1[j1]=tempre2[j]-tempim2[j+i]-tempre2[j+2*i]+tempim2[j+3*i];
					tempim1[j1]=tempim2[j]+tempre2[j+i]-tempim2[j+2*i]-tempre2[j+3*i];
					temp=tempre1[j1]*re[3*j11*N/(i*4)]-tempim1[j1]*im[3*j11*N/(i*4)];
					tempim1[j1]=tempim1[j1]*re[3*j11*N/(i*4)]+tempre1[j1]*im[3*j11*N/(i*4)];
					tempre1[j1]=temp;
					j1=j1-3*i+1;
					j11++;
				}
				j1=j1+3*i;
				}
		}
		k*=4;
		n++;
	}
	//////变址运算///////////////
	m=1;
	n=N/16;
	int i1=0;
	for(m=1;m<N/4;){
		m*=4;	
		i1++;
		int j1=0;
		if(flag==1){	
			flag=2;
		    for(k=0;k<n;k++){	
			for(j=m*k*4;j<m+m*k*4;j++){			
					tempre2[j1]=tempre1[j];
					tempim2[j1]=tempim1[j];
					j1++;
					tempre2[j1]=tempre1[j+m];
					tempim2[j1]=tempim1[j+m];
					j1++;
					tempre2[j1]=tempre1[j+2*m];
					tempim2[j1]=tempim1[j+2*m];
					j1++;
					tempre2[j1]=tempre1[j+3*m];   
					tempim2[j1]=tempim1[j+3*m];
					j1++;
			}
			}
		}else{
					flag=1;
                    for(k=0;k<n;k++){
						for(j=m*k*4;j<m+m*k*4;j++){	
					tempre1[j1]=tempre2[j];
					tempim1[j1]=tempim2[j];
					j1++;
					tempre1[j1]=tempre2[j+m];
					tempim1[j1]=tempim2[j+m];
					j1++;
					tempre1[j1]=tempre2[j+2*m];
					tempim1[j1]=tempim2[j+2*m];
					j1++;
					tempre1[j1]=tempre2[j+3*m];
					tempim1[j1]=tempim2[j+3*m];
					j1++;					
						}
					}
		}
		n/=4;
		cout<<j1;
	}
	j=0;
	for(i=0;i<N;i++){
		if(flag==1){
			mod[i]=sqrt(tempre1[i]*tempre1[i]+tempim1[i]*tempim1[i]);
		}else{
			mod[i]=sqrt(tempre2[i]*tempre2[i]+tempim2[i]*tempim2[i]);
		}
		if(mod[i]>200)
		cout<<i<<"m"<<mod[i]/512<<endl;
	}
}