#include<cmath>
#include<iomanip>
#include<iostream>
#include<time.h>
using namespace std;
const double pi = acos(-1.0);
void FFT(int NUM,double *dataR, double *dataI)
{
int i,j,k;
int r,step;
double *tempR, *tempI;
tempR = new double[NUM];
tempI = new double[NUM];
double sin_tap,cos_tap;
double real1,image1,real2,image2;
int temp = (int)(log(NUM*1.0)/log(2.0)+0.5);
/***************************following code invert sequence**********************/
for(i = 0; i<= NUM-1; i++)
{
int ttemp = 0, k = i;
for(j = 0; j <= temp-1; j++)
{
ttemp*=2;
ttemp += k%2;
k = k >> 1;
}
tempR[ttemp] = dataR[i];
tempI[ttemp] = dataI[i];
}
for(i = 0; i <= NUM-1; i++)
{
dataR[i] = tempR[i];
dataI[i] = tempI[i];
}/*整序完成*/
/***************************following code FFT*********************************/
for(i=1;i<=temp;i++)
{
step=1 << i;
for(k=0; k <= NUM-1; k+=step)
{
for(r=0;r<=step/2-1;r++)
{
cos_tap=cos(2*pi*r/step);
sin_tap=sin(-2*pi*r/step);
real1=dataR[r+k];
image1=dataI[r+k];
real2=dataR[r+k+step/2];
image2=dataI[r+k+step/2];
dataR[r+k]=real1+real2*cos_tap-image2*sin_tap;
dataI[r+k]=image1+real2*sin_tap+image2*cos_tap;
dataR[r+k+step/2]=real1-real2*cos_tap+image2*sin_tap;
dataI[r+k+step/2]=image1-real2*sin_tap-image2*cos_tap;
}
}
}
for(i = 0; i <= NUM-1; i++)
{
// cout<<setw(8)<<i<<setw(16)<<dataR[i]<<setw(16)<<dataI[i]<<endl;
}
}/*end FFT*/
/**********************************begin IFFT***************************************/
void IFFT(int NUM,double *dataR, double *dataI)
{
int i,j,k;
int r,step;
double sin_tap,cos_tap;
double *tempR, *tempI;
tempR = new double[NUM];
tempI = new double[NUM];
double real1,image1,real2,image2;
int temp = (int)(log(NUM*1.0)/log(2.0)+0.5);
/***************************following code invert sequence**********************/
for(i = 0; i<= NUM-1; i++)
{
int ttemp = 0, k = i;
for(j = 0; j <= temp-1; j++)
{
ttemp*=2;
ttemp += k%2;
k = k >> 1;
}
tempR[ttemp] = dataR[i];
tempI[ttemp] = dataI[i];
}
for(i = 0; i <= NUM-1; i++)
{
dataR[i] = tempR[i];
dataI[i] = -tempI[i];
}/*整序完成*/
/***************************following code IFFT*********************************/
for(i=1;i<=temp;i++)
{
step=1 << i;
for(k=0; k<=NUM-1; k+=step)
{
for(r=0;r<=step/2-1;r++)
{
cos_tap=cos(2*pi*r/step);
sin_tap=sin(-2*pi*r/step);
real1=dataR[r+k];
image1=dataI[r+k];
real2=dataR[r+k+step/2];
image2=dataI[r+k+step/2];
dataR[r+k]=(real1+real2*cos_tap-image2*sin_tap)/2;
dataI[r+k]=(image1+real2*sin_tap+image2*cos_tap)/2;
dataR[r+k+step/2]=(real1-real2*cos_tap+image2*sin_tap)/2;
dataI[r+k+step/2]=(image1-real2*sin_tap-image2*cos_tap)/2;
}
}
}
for(i=0;i<=NUM-1;i++)
{
// dataR[i] = dataR[i];
dataI[i] = -dataI[i];
}
}/*end IFFT*/
/****************************begin ZKJJ**********************************/
void ZKJJ(int N, int M, double *dataR1, double *dataI1, double *dataR2, double *dataI2, double *ansR, double *ansI)
{
int NUM;
if(N >= M) NUM=2*N;
else NUM=2*M;
int i;
for(i = 0; i <= NUM-1; i++) ansR[i] = ansI[i] = 0;
double *dataR3 = new double[NUM], *dataI3 = new double[NUM];
memcpy(dataR3, dataR1, sizeof(double)*N);
memcpy(dataI3, dataI1, sizeof(double)*N);
for(i = N; i <= NUM-1; i++)
{
dataR3[i] = dataI3[i] = 0;
}
FFT(NUM, dataR3, dataI3);
int j;
double *dataR4 = new double[NUM], *dataI4 = new double[NUM];
memcpy(dataR4, dataR2, sizeof(double)*M);
memcpy(dataI4, dataI2, sizeof(double)*M);
for(j = M; j <= NUM-1; j++)
{
dataR4[j] = dataI4[j] = 0;
}
FFT(NUM, dataR4, dataI4);
for(i = 0; i <= NUM-1; i++)
{
ansR[i] = dataR3[i]*dataR4[i]-dataI3[i]*dataI4[i];
ansI[i] = dataR3[i]*dataI4[i]+dataI3[i]*dataR4[i];
}
IFFT(NUM, ansR, ansI);
}//end ZKJJ
int main()
{
int Type;
int n,i;
int N,M;
double t,H; //t是周期数；H是占空比；
cout<<"输入序列一采样点数N: ";
cin>>N;
cout<<"输入序列二采样点数M: ";
cin>>M;
double *dataR1 = new double[N], *dataI1 =new double[N];
double *dataR2 = new double[M], *dataI2 =new double[M];
cout<<"第一个序列："<<endl; //选取第一个序列
cout<<"波型选择(1-正弦波 2-矩形波 3-三角波): ";
cin>>Type;
switch(Type)
{
case 1:
cout<<"正弦波"<<endl;
cout<<"周期数: ";
cin>>t;
cout<<endl;
for(n=0; n<=N-1; n++)
{
dataR1[n]=sin(2*pi*t*n/N);
dataI1[n]=0;
}
for(i=0;i<=N-1;i++){
cout<<setw(16)<<dataR1[i];
}
cout<<endl;
break;
case 2:
cout<<"矩形波"<<endl;
cout<<"周期数 占空比: ";
cin>>t;
cin>>H;
for(n=0; n<=N-1; n++)
{
float every_T;
every_T=t*n/N-(int)(t*n/N);
if(every_T<H) dataR1[n]=1;
else dataR1[n]=0;
dataI1[n]=0;
}
cout<<endl;
for(i=0;i<=N-1;i++){
cout<<setw(16)<<dataR1[i];
}
cout<<endl;
break;
case 3:
cout<<"三角波"<<endl;
cout<<"周期数: ";
cin>>t;
for(n=0 ;n<=N-1; n++)
{
float every_T;
every_T=t*n/N-(int)(t*n/N);
if(every_T<0.5) dataR1[n]=2*every_T;
else dataR1[n]=2-2*every_T;
dataI1[n]=0;
}
for(i=0;i<=N-1;i++){
cout<<setw(16)<<dataR1[i];
}
cout<<endl;
break;
default:
cout<<"错误!!(只能输入1-3)"<<endl;
} //取完第一个波
cout<<endl;
/********************************************************************************************/
cout<<"第二个序列: "<<endl; //选取第二个波
cout<<"波型选择(1-正弦波 2-矩形波 3-三角波): ";
cin>>Type;
switch(Type)
{
case 1:
cout<<"正弦波"<<endl;
cout<<"周期数: ";
cin>>t;
cout<<endl;
for(n=0; n<=M-1; n++)
{
dataR2[n]=sin(2*pi*t*n/M);
dataI2[n]=0;
}
for(i=0;i<=M-1;i++){
cout<<setw(16)<<dataR2[i];
}
cout<<endl;
break;
case 2:
cout<<"矩形波"<<endl;
cout<<"周期数 占空比: ";
cin>>t;
cin>>H;
for(n=0; n<=M-1; n++)
{
float every_T;
every_T=t*n/M-(int)(t*n/M);
if(every_T<H) dataR2[n]=1;
else dataR2[n]=0;
dataI2[n]=0;
}
cout<<endl;
for(i=0;i<=M-1;i++){
cout<<setw(16)<<dataR2[i];
}
cout<<endl;
break;
case 3:
cout<<"三角波"<<endl;
cout<<"周期数: ";
cin>>t;
for(n=0 ;n<=M-1; n++)
{
float every_T;
every_T=t*n/M-(int)(t*n/M);
if(every_T<0.5) dataR2[n]=2*every_T;
else dataR2[n]=2-2*every_T;
dataI2[n]=0;
}
for(i=0;i<=M-1;i++){
cout<<setw(16)<<dataR2[i];
}
cout<<endl;
break;
default:
cout<<"错误!!(只能输入1-3)"<<endl;
} //取完第二个序列
cout<<"***********************************************"<<endl;
int NUM;
if(N >= M) NUM=2*N;
else NUM=2*M;
double *ansR = new double[NUM], *ansI = new double[NUM];
for(i=0; i <= NUM-1; i++) ansR[i]=ansI[i]=0;
cout<<"ZKJJ:"<<endl;
double e = clock();
ZKJJ(N, M, dataR1, dataI1, dataR2, dataI2, ansR, ansI);
double s = clock();
for(i = 0; i <= N+M-2; i++)//截断输出
{
cout<<setw(8)<<i<<setw(16)<<ansR[i]<<endl;
}
cout<<"time: ";
printf("%lf",e-s);
cout<<"ms"<<endl;
return 0;
}
评论7
最新资源