📄 fft.cpp
字号:
m2_im=c3_2*(zRe[2] - zRe[1]);
s1_re=zRe[0] + m1_re; s1_im=zIm[0] + m1_im;
zRe[1]=s1_re + m2_re; zIm[1]=s1_im + m2_im;
zRe[2]=s1_re - m2_re; zIm[2]=s1_im - m2_im;
break;
case 4 : t1_re=zRe[0] + zRe[2]; t1_im=zIm[0] + zIm[2];
t2_re=zRe[1] + zRe[3]; t2_im=zIm[1] + zIm[3];
m2_re=zRe[0] - zRe[2]; m2_im=zIm[0] - zIm[2];
m3_re=zIm[1] - zIm[3]; m3_im=zRe[3] - zRe[1];
zRe[0]=t1_re + t2_re; zIm[0]=t1_im + t2_im;
zRe[2]=t1_re - t2_re; zIm[2]=t1_im - t2_im;
zRe[1]=m2_re + m3_re; zIm[1]=m2_im + m3_im;
zRe[3]=m2_re - m3_re; zIm[3]=m2_im - m3_im;
break;
case 5 : t1_re=zRe[1] + zRe[4]; t1_im=zIm[1] + zIm[4];
t2_re=zRe[2] + zRe[3]; t2_im=zIm[2] + zIm[3];
t3_re=zRe[1] - zRe[4]; t3_im=zIm[1] - zIm[4];
t4_re=zRe[3] - zRe[2]; t4_im=zIm[3] - zIm[2];
t5_re=t1_re + t2_re; t5_im=t1_im + t2_im;
zRe[0]=zRe[0] + t5_re; zIm[0]=zIm[0] + t5_im;
m1_re=c5_1*t5_re; m1_im=c5_1*t5_im;
m2_re=c5_2*(t1_re - t2_re);
m2_im=c5_2*(t1_im - t2_im);
m3_re=-c5_3*(t3_im + t4_im);
m3_im=c5_3*(t3_re + t4_re);
m4_re=-c5_4*t4_im; m4_im=c5_4*t4_re;
m5_re=-c5_5*t3_im; m5_im=c5_5*t3_re;
s3_re=m3_re - m4_re; s3_im=m3_im - m4_im;
s5_re=m3_re + m5_re; s5_im=m3_im + m5_im;
s1_re=zRe[0] + m1_re; s1_im=zIm[0] + m1_im;
s2_re=s1_re + m2_re; s2_im=s1_im + m2_im;
s4_re=s1_re - m2_re; s4_im=s1_im - m2_im;
zRe[1]=s2_re + s3_re; zIm[1]=s2_im + s3_im;
zRe[2]=s4_re + s5_re; zIm[2]=s4_im + s5_im;
zRe[3]=s4_re - s5_re; zIm[3]=s4_im - s5_im;
zRe[4]=s2_re - s3_re; zIm[4]=s2_im - s3_im;
break;
case 8 : fft8(); break;
case 10 : fft10(); break;
default : fftodd(radix); break;
}
adr=groupOffset;
for (blockNo=0; blockNo<radix; blockNo++)
{
yRe[adr]=zRe[blockNo]; yIm[adr]=zIm[blockNo];
adr=adr+sofarRadix;
}
groupOffset=groupOffset+sofarRadix*radix;
adr=groupOffset;
}
dataOffset=dataOffset+1;
groupOffset=dataOffset;
adr=groupOffset;
}
}
void FFT::initTrig(int radix)
{
int i;
double w,xre,xim;
w=2*pi/radix;
trigRe[0]=1; trigIm[0]=0;
xre=cos(w);
xim=-sin(w);
trigRe[1]=xre; trigIm[1]=xim;
for (i=2; i<radix; i++)
{
trigRe[i]=xre*trigRe[i-1] - xim*trigIm[i-1];
trigIm[i]=xim*trigRe[i-1] + xre*trigIm[i-1];
}
} /* initTrig */
void FFT::fft4(double xRe[], double xIm[])
{
double t1_re,t1_im, t2_re,t2_im;
double m2_re,m2_im, m3_re,m3_im;
t1_re=xRe[0] + xRe[2]; t1_im=xIm[0] + xIm[2];
t2_re=xRe[1] + xRe[3]; t2_im=xIm[1] + xIm[3];
m2_re=xRe[0] - xRe[2]; m2_im=xIm[0] - xIm[2];
m3_re=xIm[1] - xIm[3]; m3_im=xRe[3] - xRe[1];
xRe[0]=t1_re + t2_re; xIm[0]=t1_im + t2_im;
xRe[2]=t1_re - t2_re; xIm[2]=t1_im - t2_im;
xRe[1]=m2_re + m3_re; xIm[1]=m2_im + m3_im;
xRe[3]=m2_re - m3_re; xIm[3]=m2_im - m3_im;
}
void FFT::fft5(double xRe[], double xIm[])
{
double t1_re,t1_im, t2_re,t2_im, t3_re,t3_im;
double t4_re,t4_im, t5_re,t5_im;
double m2_re,m2_im, m3_re,m3_im, m4_re,m4_im;
double m1_re,m1_im, m5_re,m5_im;
double s1_re,s1_im, s2_re,s2_im, s3_re,s3_im;
double s4_re,s4_im, s5_re,s5_im;
t1_re=xRe[1] + xRe[4]; t1_im=xIm[1] + xIm[4];
t2_re=xRe[2] + xRe[3]; t2_im=xIm[2] + xIm[3];
t3_re=xRe[1] - xRe[4]; t3_im=xIm[1] - xIm[4];
t4_re=xRe[3] - xRe[2]; t4_im=xIm[3] - xIm[2];
t5_re=t1_re + t2_re; t5_im=t1_im + t2_im;
xRe[0]=xRe[0] + t5_re; xIm[0]=xIm[0] + t5_im;
m1_re=c5_1*t5_re; m1_im=c5_1*t5_im;
m2_re=c5_2*(t1_re - t2_re); m2_im=c5_2*(t1_im - t2_im);
m3_re=-c5_3*(t3_im + t4_im); m3_im=c5_3*(t3_re + t4_re);
m4_re=-c5_4*t4_im; m4_im=c5_4*t4_re;
m5_re=-c5_5*t3_im; m5_im=c5_5*t3_re;
s3_re=m3_re - m4_re; s3_im=m3_im - m4_im;
s5_re=m3_re + m5_re; s5_im=m3_im + m5_im;
s1_re=xRe[0] + m1_re; s1_im=xIm[0] + m1_im;
s2_re=s1_re + m2_re; s2_im=s1_im + m2_im;
s4_re=s1_re - m2_re; s4_im=s1_im - m2_im;
xRe[1]=s2_re + s3_re; xIm[1]=s2_im + s3_im;
xRe[2]=s4_re + s5_re; xIm[2]=s4_im + s5_im;
xRe[3]=s4_re - s5_re; xIm[3]=s4_im - s5_im;
xRe[4]=s2_re - s3_re; xIm[4]=s2_im - s3_im;
}
void FFT::fft8()
{
double aRe[4], aIm[4], bRe[4], bIm[4], gem;
aRe[0] = zRe[0]; bRe[0] = zRe[1];
aRe[1] = zRe[2]; bRe[1] = zRe[3];
aRe[2] = zRe[4]; bRe[2] = zRe[5];
aRe[3] = zRe[6]; bRe[3] = zRe[7];
aIm[0] = zIm[0]; bIm[0] = zIm[1];
aIm[1] = zIm[2]; bIm[1] = zIm[3];
aIm[2] = zIm[4]; bIm[2] = zIm[5];
aIm[3] = zIm[6]; bIm[3] = zIm[7];
fft4(aRe, aIm); fft4(bRe, bIm);
gem = c8*(bRe[1] + bIm[1]);
bIm[1] = c8*(bIm[1] - bRe[1]);
bRe[1] = gem;
gem = bIm[2];
bIm[2] =-bRe[2];
bRe[2] = gem;
gem = c8*(bIm[3] - bRe[3]);
bIm[3] =-c8*(bRe[3] + bIm[3]);
bRe[3] = gem;
zRe[0] = aRe[0] + bRe[0]; zRe[4] = aRe[0] - bRe[0];
zRe[1] = aRe[1] + bRe[1]; zRe[5] = aRe[1] - bRe[1];
zRe[2] = aRe[2] + bRe[2]; zRe[6] = aRe[2] - bRe[2];
zRe[3] = aRe[3] + bRe[3]; zRe[7] = aRe[3] - bRe[3];
zIm[0] = aIm[0] + bIm[0]; zIm[4] = aIm[0] - bIm[0];
zIm[1] = aIm[1] + bIm[1]; zIm[5] = aIm[1] - bIm[1];
zIm[2] = aIm[2] + bIm[2]; zIm[6] = aIm[2] - bIm[2];
zIm[3] = aIm[3] + bIm[3]; zIm[7] = aIm[3] - bIm[3];
}
void FFT::fft10()
{
double aRe[5], aIm[5], bRe[5], bIm[5];
aRe[0] = zRe[0]; bRe[0] = zRe[5];
aRe[1] = zRe[2]; bRe[1] = zRe[7];
aRe[2] = zRe[4]; bRe[2] = zRe[9];
aRe[3] = zRe[6]; bRe[3] = zRe[1];
aRe[4] = zRe[8]; bRe[4] = zRe[3];
aIm[0] = zIm[0]; bIm[0] = zIm[5];
aIm[1] = zIm[2]; bIm[1] = zIm[7];
aIm[2] = zIm[4]; bIm[2] = zIm[9];
aIm[3] = zIm[6]; bIm[3] = zIm[1];
aIm[4] = zIm[8]; bIm[4] = zIm[3];
fft5(aRe, aIm); fft5(bRe, bIm);
zRe[0] = aRe[0] + bRe[0]; zRe[5] = aRe[0] - bRe[0];
zRe[6] = aRe[1] + bRe[1]; zRe[1] = aRe[1] - bRe[1];
zRe[2] = aRe[2] + bRe[2]; zRe[7] = aRe[2] - bRe[2];
zRe[8] = aRe[3] + bRe[3]; zRe[3] = aRe[3] - bRe[3];
zRe[4] = aRe[4] + bRe[4]; zRe[9] = aRe[4] - bRe[4];
zIm[0] = aIm[0] + bIm[0]; zIm[5] = aIm[0] - bIm[0];
zIm[6] = aIm[1] + bIm[1]; zIm[1] = aIm[1] - bIm[1];
zIm[2] = aIm[2] + bIm[2]; zIm[7] = aIm[2] - bIm[2];
zIm[8] = aIm[3] + bIm[3]; zIm[3] = aIm[3] - bIm[3];
zIm[4] = aIm[4] + bIm[4]; zIm[9] = aIm[4] - bIm[4];
}
void FFT::fftodd(int radix)
{
double rere, reim, imre, imim;
int i,j,k,n,max;
n = radix;
max = (n + 1)/2;
for (j=1; j < max; j++)
{
vRe[j] = zRe[j] + zRe[n-j];
vIm[j] = zIm[j] - zIm[n-j];
wRe[j] = zRe[j] - zRe[n-j];
wIm[j] = zIm[j] + zIm[n-j];
}
for (j=1; j < max; j++)
{
zRe[j]=zRe[0];
zIm[j]=zIm[0];
zRe[n-j]=zRe[0];
zIm[n-j]=zIm[0];
k=j;
for (i=1; i < max; i++)
{
rere = trigRe[k] * vRe[i];
imim = trigIm[k] * vIm[i];
reim = trigRe[k] * wIm[i];
imre = trigIm[k] * wRe[i];
zRe[n-j] += rere + imim;
zIm[n-j] += reim - imre;
zRe[j] += rere - imim;
zIm[j] += reim + imre;
k = k + j;
if (k >= n) k = k - n;
}
}
for (j=1; j < max; j++)
{
zRe[0]=zRe[0] + vRe[j];
zIm[0]=zIm[0] + wIm[j];
}
}
double FFT::leak()
{
int i,k=3;
double sqrsum, sqrmax, sqr;
sqrsum = (yRe[k]-length)*(yRe[k]-length) + yIm[k]*yIm[k];
sqrmax = 0;
for (i=0; i<length; i++)
if (i != k)
{
sqr = yRe[i]*yRe[i] + yIm[i]*yIm[i];
sqrsum = sqrsum + sqr;
if (sqr > sqrmax)
sqrmax = sqr;
}
sqrsum = sqrsum/length; /* division by n*n in two stages to avoid */
sqrsum = sqrsum/length; /* integer overflow !!!! */
sqrmax = sqrmax/length;
sqrmax = sqrmax/length;
if (sqrsum > 0)
sqrsum = 10*log10(sqrsum);
else
sqrsum = -1000;
if (sqrmax > 0)
sqrmax = 10*log10(sqrmax);
else
sqrmax = -1000;
return(sqrsum);
}
void FFT::outputX()
{
int i;
for(i=0;i<length;i++)
cout<<xRe[i]<<" "<<xIm[i]<<endl;
}
void FFT::outputY()
{
int i;
for(i=0;i<length;i++)
cout<<yRe[i]<<" "<<yIm[i]<<endl;
}⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -