📄 fft.cpp
字号:
#include"COMPLEX_.H"
#include<math.h>
COMPLEX Cadd(COMPLEX z1,COMPLEX z2)
{
z1.x=z1.x+z2.x;
z1.y=z1.y+z2.y;
return z1;
};
// ======define the complex substraction between two complex numbers=====
COMPLEX Csub(COMPLEX z1,COMPLEX z2)
{
z1.x=z1.x-z2.x;
z1.y=z1.y-z2.y;
return z1;
};
// ======define the complex multiplication between two complex numbers=====
COMPLEX Cmul(COMPLEX z1,COMPLEX z2)
{
COMPLEX mul;
mul.x=z1.x*z2.x-z1.y*z2.y;
mul.y=z1.x*z2.y+z1.y*z2.x;
return mul;
} ;
// ======define the complex division between two complex numbers=====
COMPLEX Cdiv(COMPLEX z1,COMPLEX z2)
{
COMPLEX div;
double t;
t=z2.x*z2.x+z2.y*z2.y;
if(t!=0.)
{
div.x=(z1.x*z2.x+z1.y*z2.y)/t;
div.y=(z1.y*z2.x-z1.x*z2.y)/t;
return div;
}
else
{
div.x=0.0;
div.y=0.0;
return div;
}
};
// =====define the conjution of a complex number====
COMPLEX Cconj(COMPLEX z)
{
z.y=-z.y;
return z;
};
// =====define the magnitude of a complex number=====
double Cabs(COMPLEX z)
{
double mod;
mod=sqrt(z.x*z.x+z.y*z.y);
return mod;
};
// =====define the phasor of a complex number=====
double Cpha(COMPLEX z)
{
double phase;
phase=atan2(z.y, z.x);
return phase;
};
// ======define the square of a complex number's magnitude=====
double Csqu(COMPLEX z)
{
double squ;
squ=z.x*z.x+z.y*z.y;
return squ;
};
// =====This function is used to compute the fft of a complex sequence=====
// =====that is shorter than 4096-point=====
void fft1(int N,COMPLEX *x)
{
int i,j,k,l,LE,LE1,IP,M,lll;
COMPLEX T,U,W;
lll=1;
for(i=1;i<16;i++)
{
lll*=2;
if(lll==N)break;
}
M=i;
j=0;
for(i=0;i<N-1;i++)
{
if(i<j)
{
T=x[j];
x[j]=x[i];
x[i]=T;
}
k=N/2;
while(k<=j)
{
j=j-k;
k=k/2;
}
j=j+k;
}
for(l=1;l<=M;l++)
{
LE=1<<l;
LE1=LE/2;
U.x=1.0;
U.y=0.0;
W.x=cos(PI/LE1);
W.y=-sin(PI/LE1);
for(j=0;j<LE1;j++)
{
for(i=j;i<N;i=i+LE)
{
IP=i+LE1;
T=Cmul(x[IP],U);
x[IP]=Csub(x[i],T);
x[i]=Cadd(x[i],T);
}
U=Cmul(U,W);
}
}
}
// =====This function is used to compute the fft of a complex series=====
// =====that equals to 8192-point=====
void fft2(COMPLEX *x,COMPLEX *y)
{
int i,j,k,l,LE,LE1,IP,M,N;
COMPLEX U,W,T;
N=8192;
M=13;
j=0;
for(i=0;i<N-1;i++)
{
if(i<j)
{
if(j>=4096&&i>=4096)
{
T=y[j-4096];
y[j-4096]=y[i-4096];
y[i-4096]=T;
}
else if(j>=4096&&i<4096)
{
T=y[j-4096];
y[j-4096]=x[i];
x[i]=T;
}
else
{
T=x[j];
x[j]=x[i];
x[i]=T;
}
}
k=4096;
while(k<=j)
{
j=j-k;
k=k/2;
}
j=j+k;
}
for(l=1;l<=M;l++)
{
LE=1<<l;
LE1=LE/2;
U.x=1.0;
U.y=0.0;
W.x=cos(PI/LE1);
W.y=-sin(PI/LE1);
for(j=0;j<LE1;j++)
{
for(i=j;i<N;i=i+LE)
{
IP=i+LE1;
if(i>=4096)
{
T=Cmul(y[IP-4096],U);
y[IP-4096]=Csub(y[i-4096],T);
y[i-4096]=Cadd(y[i-4096],T);
}
else if(IP>=4096&&i<4096)
{
T=Cmul(y[IP-4096],U);
y[IP-4096]=Csub(x[i],T);
x[i]=Cadd(x[i],T);
}
else
{
T=Cmul(x[IP],U);
x[IP]=Csub(x[i],T);
x[i]=Cadd(x[i],T);
}
}
U=Cmul(U,W);
}
}
}
// =====The function of the N-point fft of an N-point real series=====
int realfft(int N,COMPLEX *x1,COMPLEX *x2,COMPLEX *x3,COMPLEX *x4)
{
int i,lll;
COMPLEX U,W,T,P,V,G;
T.x=0.5;
T.y=0;
P.x=0;
P.y=-0.5;
V.x=cos(2.0*PI/N);
V.y=-sin(2.0*PI/N);
G.x=1.0;
G.y=0.0;
if(N<=0||N>16384)return 2;
lll=1;
for(i=1;i<15;i++)
{
lll*=2;
if(lll==N)break;
}
if(i>14)return 1;
// ======when series shorter than 4096-point=====
if(N<=4096)
{
for(i=0;i<N/2;i++)
{
x1[i].y=x1[2*i+1].x;
x1[i].x=x1[2*i].x;
}
fft1(N/2,x1);
x1[N/2]=Cmul(Csub(x1[0],Cconj(x1[0])),P);
x1[0]=Cmul(Cadd(x1[0],Cconj(x1[0])),T);
x1[N/4*3]=Cmul(Csub(x1[N/4],Cconj(x1[N/4])),P);
x1[N/4]=Cmul(Cadd(x1[N/4],Cconj(x1[N/4])),T);
for(i=1;i<N/4;i++)
{
x1[i+N/2]=Cmul(Csub(x1[i],Cconj(x1[N/2-i])),P);
x1[N-i]=Cmul(Csub(x1[N/2-i],Cconj(x1[i])),P);
U=Cadd(x1[i],Cconj(x1[N/2-i]));
W=Cadd(x1[N/2-i],Cconj(x1[i]));
x1[i]=Cmul(U,T);
x1[N/2-i]=Cmul(W,T);
}
for(i=0;i<N/2;i++)
{
W=Cmul(x1[N/2+i],G);
x1[N/2+i]=Csub(x1[i],W);
x1[i]=Cadd(x1[i],W);
G=Cmul(G,V);
}
return 0;
}
// ======when series equals to 8192-point=====
if(N==8192)
{
for(i=0;i<4096;i++)
{
if(i<2048)
{
x1[i].y=x1[2*i+1].x;
x1[i].x=x1[2*i].x;
}
else
{
x1[i].y=x2[2*i+1-4096].x;
x1[i].x=x2[2*i-4096].x;
}
}
fft1(4096,x1);
U=Cadd(x1[0],Cconj(x1[0]));
W=Csub(x1[0],Cconj(x1[0]));
x1[0]=Cmul(U,T);
x2[0]=Cmul(W,P);
U=Cadd(x1[2048],Cconj(x1[2048]));
W=Csub(x1[2048],Cconj(x1[2048]));
x1[2048]=Cmul(U,T);
x2[2048]=Cmul(W,P);
for(i=1;i<2048;i++)
{
U=Csub(x1[i],Cconj(x1[4096-i]));
W=Csub(x1[4096-i],Cconj(x1[i]));
x2[i]=Cmul(U,P);
x2[4096-i]=Cmul(W,P);
U=Cadd(x1[i],Cconj(x1[4096-i]));
W=Cadd(x1[4096-i],Cconj(x1[i]));
x1[i]=Cmul(U,T);
x1[4096-i]=Cmul(W,T);
}
for(i=0;i<4096;i++)
{
W=Cmul(x2[i],G);
x2[i]=Csub(x1[i],W);
x1[i]=Cadd(x1[i],W);
G=Cmul(G,V);
}
return 0;
}
// ======when series equals to 16384-point=====
if(N==16384)
{
for(i=0;i<8192;i++)
{
if(i<2048)
{
x1[i].y=x1[2*i+1].x;
x1[i].x=x1[2*i].x;
}
else if(i>=2048&&i<4096)
{
x1[i].y=x2[2*i+1-4096].x;
x1[i].x=x2[2*i-4096].x;
}
else if(i>=4096&&i<6144)
{
x2[i-4096].y=x3[2*i+1-8192].x;
x2[i-4096].x=x3[2*i-8192].x;
}
else
{
x2[i-4096].y=x4[2*i+1-12288].x;
x2[i-4096].x=x4[2*i-12288].x;
}
}
fft2(x1,x2);
U=Cadd(x1[0],Cconj(x1[0]));
W=Csub(x1[0],Cconj(x1[0]));
x1[0]=Cmul(U,T);
x3[0]=Cmul(W,P);
U=Cadd(x2[0],Cconj(x2[0]));
W=Csub(x2[0],Cconj(x2[0]));
x2[0]=Cmul(U,T);
x4[0]=Cmul(W,P);
for(i=1;i<4096;i++)
{
U=Csub(x1[i],Cconj(x2[4096-i]));
W=Csub(x2[4096-i],Cconj(x1[i]));
x3[i]=Cmul(U,P);
x4[4096-i]=Cmul(W,P);
U=Cadd(x1[i],Cconj(x2[4096-i]));
W=Cadd(x2[4096-i],Cconj(x1[i]));
x1[i]=Cmul(U,T);
x2[4096-i]=Cmul(W,T);
}
for(i=0;i<8192;i++)
{
if(i<4096)
{
W=Cmul(x3[i],G);
x3[i]=Csub(x1[i],W);
x1[i]=Cadd(x1[i],W);
G=Cmul(G,V);
}
else
{
W=Cmul(x4[i-4096],G);
x4[i-4096]=Csub(x2[i-4096],W);
x2[i-4096]=Cadd(x2[i-4096],W);
G=Cmul(G,V);
}
}
return 0;
}
}
int mag_pha(int N,COMPLEX *x1,COMPLEX *x2,COMPLEX *x3,COMPLEX *x4)
{
COMPLEX t;
int i;
if(N<=0||N>16384)
{
return 1;;
}
for(i=0;i<N;i++)
{
if(i<4096&&(x1[i].x!=0.||x1[i].y!=0.))
{
t=x1[i];
x1[i].x=Cabs(x1[i]);
x1[i].y=Cpha(t);
}
if((i>=4096&&i<8192)&&(x2[i-4096].x!=0.||x2[i-4096].y!=0.))
{
t=x2[i-4096];
x2[i-4096].x=Cabs(x2[i-4096]);
x2[i-4096].y=Cpha(t);
}
if((i>=8192&&i<12288)&&(x3[i-8192].x!=0.||x3[i-8192].y!=0.))
{
t=x3[i-8192];
x3[i-8192].x=Cabs(x3[i-8192]);
x3[i-8192].y=Cpha(t);
}
if((i>=12288&&i<16384)&&(x4[i-12288].x!=0.||x4[i-12288].y!=0.))
{
t=x4[i-12288];
x4[i-12288].x=Cabs(x4[i-12288]);
x4[i-12288].y=Cpha(t);
}
}
return 0;
}
int maxas1(int N,COMPLEX *x1,COMPLEX *x2,COMPLEX *x3,COMPLEX *x4)
{
COMPLEX t;
float max1,max2;
int i;
int num;
max1=x1[0].x;max2=fabs(x1[0].y);
for(i=0;i<N;i++)
{
if(i<4096)
{
if(x1[i].x>max1)
{
max1=x1[i].x;
num=i;
}
if(fabs(x1[i].y)>max2)max2=fabs(x1[i].y);
}
if(i>=4096&&i<8192)
{
if(x2[i-4096].x>max1)
{
max1=x2[i-4096].x;
num=i;
}
if(fabs(x2[i-4096].y)>max2)max2=fabs(x2[i-4096].y);
}
if(i>=8192&&i<12288)
{
if(x3[i-8192].x>max1)
{
max1=x3[i-8192].x;
num=i;
}
if(fabs(x3[i-8192].y)>max2)max2=fabs(x3[i-8192].y);
}
if(i>=12288&&i<16384)
{
if(x4[i-12288].x>max1)
{
max1=x4[i-12288].x;
num=i;
}
if(fabs(x4[i-12288].y)>max2)max2=fabs(x4[i-12288].y);
}
}
if(max1!=0&&max2!=0)
{
for(i=0;i<N;i++)
{
if(i<4096)
{
x1[i].x=x1[i].x/max1;
x1[i].y=x1[i].y/max2;
}
if(i>=4096&&i<8192)
{
x2[i-4096].x=x2[i-4096].x/max1;
x2[i-4096].y=x2[i-4096].y/max2;
}
if(i>=8192&&i<12288)
{
x3[i-8192].x=x3[i-8192].x/max1;
x3[i-8192].y=x3[i-8192].y/max2;
}
if(i>=12288&&i<16384)
{
x4[i-12288].x=x4[i-12288].x/max1;
x4[i-12288].y=x4[i-12288].y/max2;
}
}
}
return num;
}
// =====This function is used to compute the bessel fuction=====
// =====beesel(x)=1+(1/1!)^2/(x/2.0)^(2*1)+(1/2!)^2/(x/2.0)^(2*2)+....=====
double bessel(double x)
{
int i,j;
double t,y,z;
t=0.0;
y=1.0;
z=1.0;
for(i=1;i<10000;i++)
{
for(j=1;j<=i;j++)
z=1.0/j*z;
y=y+z*z*pow(x/2.0,2.0*i);
if(t==y) break;
t=y;
}
x=x*y;
return x;
}
// =====A function to add Bartlett window to a complex series=====
void addwin2(int N,COMPLEX *x1,COMPLEX *x2,COMPLEX* x3,COMPLEX *x4)
{
int i;
double h;
// ===== while N is shorter than 4096 =====
if(N<=8192)
{
for(i=0;i<N;i++)
{
if(i<N/2)
{
h=i*2.0/N;
x1[i].x=x1[i].x*h;
x1[i].y=x1[i].y*h;
}
else if(i>=N/2&&i<4096)
{
h=(N-i)*2.0/N;
x1[i].x=x1[i].x*h;
x1[i].y=x1[i].y*h;
}
else
{
h=(N-i)*2.0/N;
x2[i-4096].x=x2[i-4096].x*h;
x2[i-4096].y=x2[i-4096].y*h;
}
}
}
// ===== while N is longer than 4096 and shorter than 8192 =====
if(N>8192&&N<=16384)
{
for(i=0;i<N;i++)
{
if(i<4096)
{
h=i*2.0/N;
x1[i].x=x1[i].x*h;
x1[i].y=x1[i].y*h;
}
else if(i>=4096&&i<N/2)
{
h=i*2.0/N;
x2[i-4096].x=x2[i-4096].x*h;
x2[i-4096].y=x2[i-4096].y*h;
}
else if(i>=N/2&&i<8192)
{
h=(N-i)*2.0/N;
x2[i-4096].x=x2[i-4096].x*h;
x2[i-4096].y=x2[i-4096].y*h;
}
else if(i>=8192&&i<12288)
{
h=(N-i)*2.0/N;
x3[i-8192].x=x3[i-8192].x*h;
x3[i-8192].y=x3[i-8192].y*h;
}
else
{
h=(N-i)*2.0/N;
x4[i-12288].x=x4[i-12288].x*h;
x4[i-12288].y=x4[i-12288].y*h;
}
}
}
}
// ===== A function to add all kinds of windows to a complex series =====
int addwindow(int N,int j,COMPLEX *x1,COMPLEX *x2,COMPLEX *x3,COMPLEX *x4)
{
int i,l,m;
COMPLEX T,W,U,V;
double h;
if(N>16384||N<=0)
return 1;
switch(j)
{
// ===== Add a rectangle window =====
case 1:
break;
// ===== Add a Bartlett window =====
case 2:
addwin2(N,x1,x2,x3,x4);
break;
// ===== default Consin type2: w(n)=[sin(n*PI/N)]^2 =====
// ===== It's called Hanning window =====
case 3:
U.x=cos(2*PI/N);
U.y=sin(2*PI/N);
W.x=1.0;
W.y=0.0;
for(i=0;i<N;i++)
{
h=0.5-0.5*W.x;
if(i<4096)
{
x1[i].x=x1[i].x*h;
x1[i].y=x1[i].y*h;
}
if(i>=4096&&i<8192)
{
x2[i-4096].x=x2[i-4096].x*h;
x2[i-4096].y=x2[i-4096].y*h;
}
if(i>=8192&&i<12288)
{
x3[i-8192].x=x3[i-8192].x*h;
x3[i-8192].y=x3[i-8192].y*h;
}
if(i>=12288&&i<16384)
{
x4[i-12288].x=x4[i-12288].x*h;
x4[i-12288].y=x4[i-12288].y*h;
}
W=Cmul(U,W);
}
break;
// ===== Add Windows of modified Hanning type =====
// ===== w(n)=t-(1-t)cos(2*PI*n/N) here 0<t<=1 =====
case 4:
U.x=cos(2*PI/(N-1));
U.y=sin(2*PI/(N-1));
W.x=1.0;
W.y=0.0;
for(i=0;i<N;i++)
{
h=0.54-(1.0-0.54)*W.x;
if(i<4096)
{
x1[i].x=x1[i].x*h;
x1[i].y=x1[i].y*h;
}
else if(i>=4096&&i<8192)
{
x2[i-4096].x=x2[i-4096].x*h;
x2[i-4096].y=x2[i-4096].y*h;
}
else if(i>=8192&&i<12288)
{
x3[i-8192].x=x3[i-8192].x*h;
x3[i-8192].y=x3[i-8192].y*h;
}
else if(i>=12288&&i<16384)
{
x4[i-12288].x=x4[i-12288].x*h;
x4[i-12288].y=x4[i-12288].y*h;
}
W=Cmul(U,W);
}
break;
// ===== Add Windows of Blackman type =====
case 5:
U.x=cos(2*PI/N);
U.y=sin(2*PI/N);
W.x=1.0;
W.y=0.0;
T=W;
V=W;
for(i=0;i<N;i++)
{
h=7938.0/18608.0-9240.0/18608.0*W.x;
h=h+1430.0/18608.0*T.x;
if(i<4096)
{
x1[i].x=x1[i].x*h;
x1[i].y=x1[i].y*h;
}
else if(i>=4096&&i<8192)
{
x2[i-4096].x=x2[i-4096].x*h;
x2[i-4096].y=x2[i-4096].y*h;
}
else if(i>=8192&&i<12288)
{
x3[i-8192].x=x3[i-8192].x*h;
x3[i-8192].y=x3[i-8192].y*h;
}
else if(i>=12288&&i<16384)
{
x4[i-12288].x=x4[i-12288].x*h;
x4[i-12288].y=x4[i-12288].y*h;
}
W=Cmul(U,W);
T=Cmul(Cmul(U,U),T);
}
break;
// ===== Add a Gaussin Window =====
// ===== w(n)=exp(0.5*t^2*(n/(N/2)-1)^2 =====
case 6:
for(i=0;i<N;i++)
{
h=exp(-0.5*9*(i*2.0/N-1.0)*(i*2.0/N-1.0));
if(i<4096)
{
x1[i].x=x1[i].x*h;
x1[i].y=x1[i].y*h;
}
else if(i>=4096&&i<8192)
{
x2[i-4096].x=x2[i-4096].x*h;
x2[i-4096].y=x2[i-4096].y*h;
}
else if(i>=8192&&i<12288)
{
x3[i-8192].x=x3[i-8192].x*h;
x3[i-8192].y=x3[i-8192].y*h;
}
else if(i>=12288&&i<16384)
{
x4[i-12288].x=x4[i-12288].x*h;
x4[i-12288].y=x4[i-12288].y*h;
}
}
break;
// ===== Add a Kaiser Window =====
case 7:
for(i=0;i<N;i++)
{
h=bessel(7.68*sqrt(1-(1.0-2.0*i/N)*(1.0-2.0*i/N)))/bessel(7.68);
if(i<4096)
{
x1[i].x=x1[i].x*h;
x1[i].y=x1[i].y*h;
}
else if(i>=4096&&i<8192)
{
x2[i-4096].x=x2[i-4096].x*h;
x2[i-4096].y=x2[i-4096].y*h;
}
else if(i>=8192&&i<12288)
{
x3[i-8192].x=x3[i-8192].x*h;
x3[i-8192].y=x3[i-8192].y*h;
}
else if(i>=12288&&i<16384)
{
x4[i-12288].x=x4[i-12288].x*h;
x4[i-12288].y=x4[i-12288].y*h;
}
}
break;
default:
break;
}
return 0;
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -