?? fft842.c
字號:
void FFT842(float *pReal, float *pImag, int Nthpo, int TransformSign)
{
// 快速傅里葉變換:8 4 2
// 若 TransformSign = 1 則置換后的pReal、pImag[0...Nthpo-1]是輸入的離散傅里葉變換;
// 若 TransformSign = 0 則置換后的pReal、pImag[0...Nthpo-1]是輸入的離散傅里逆葉變換;
// pReal、pImag是界長為Nthpo的復數組。
// Nthpo必須是2的整數冪(這點未檢驗!)
int i, j, k, n, length;
int N2pow=1; while((1<<(++N2pow)) < Nthpo);
int N8pow=(N2pow/3)*3;
if(TransformSign) for(i=0; i<Nthpo; i++) pImag[i] = -pImag[i];
double p7 =0.5*sqrt(2.0), scale;
double c1, c2, c3, c4, c5, c6, c7;
double s1, s2, s3, s4, s5, s6, s7;
float r0, r1, r2, r3, r4, r5, r6, r7;
float i0, i1, i2, i3, i4, i5, i6, i7;
double tr, ti;
float *R0, *R1, *R2, *R3, *R4, *R5, *R6, *R7;
float *I0, *I1, *I2, *I3, *I4, *I5, *I6, *I7;
for(i=3; i <= N8pow; i+=3)
{
// 基8
length = (n = 1<<(N2pow-i))<<3;
scale = atan(1.0)/n;
for(j=0; j<n; j++)
{
c1=cos(scale*j), s1=sin(scale*j);
c2=c1*c1-s1*s1, s2=c1*s1+c1*s1;
c3=c1*c2-s1*s2, s3=c2*s1+s2*c1;
c4=c2*c2-s2*s2, s4=c2*s2+c2*s2;
c5=c2*c3-s2*s3, s5=c3*s2+s3*c2;
c6=c3*c3-s3*s3, s6=c3*s3+c3*s3;
c7=c3*c4-s3*s4, s7=c4*s3+s4*c3;
for(k=j; k<Nthpo; k+=length)
{
R7=(R6=(R5=(R4=(R3=(R2=(R1=(R0=pReal+k)+n)+n)+n)+n)+n)+n)+n;
I7=(I6=(I5=(I4=(I3=(I2=(I1=(I0=pImag+k)+n)+n)+n)+n)+n)+n)+n;
r0=(*R0+*R4)+(*R2+*R6), i0=(*I0+*I4)+(*I2+*I6);
r1=(*R1+*R5)+(*R3+*R7), i1=(*I1+*I5)+(*I3+*I7);
r2=(*R0+*R4)-(*R2+*R6), i2=(*I0+*I4)-(*I2+*I6);
r3=(*R1+*R5)-(*R3+*R7), i3=(*I1+*I5)-(*I3+*I7);
r4=(*R0-*R4)-(*I2-*I6), i4=(*I0-*I4)+(*R2-*R6);
r5=(*R1-*R5)-(*I3-*I7), i5=(*I1-*I5)+(*R3-*R7);
r6=(*R0-*R4)+(*I2-*I6), i6=(*I0-*I4)-(*R2-*R6);
r7=(*R1-*R5)+(*I3-*I7), i7=(*I1-*I5)-(*R3-*R7);
*R0=r0+r1, *I0=i0+i1;
if(j)
{
*R1=(float)(c4*(r0-r1)-s4*(i0-i1)), *I1=(float)(c4*(i0-i1)+s4*(r0-r1));
*R2=(float)(c2*(r2-i3)-s2*(i2+r3)), *I2=(float)(c2*(i2+r3)+s2*(r2-i3));
*R3=(float)(c6*(r2+i3)-s6*(i2-r3)), *I3=(float)(c6*(i2-r3)+s6*(r2+i3));
tr = p7*(r5-i5), ti = p7*(r5+i5);
*R4=(float)(c1*(r4+tr)-s1*(i4+ti)), *I4=(float)(c1*(i4+ti)+s1*(r4+tr));
*R5=(float)(c5*(r4-tr)-s5*(i4-ti)), *I5=(float)(c5*(i4-ti)+s5*(r4-tr));
tr = -p7*(r7+i7), ti = p7*(r7-i7);
*R6=(float)(c3*(r6+tr)-s3*(i6+ti)), *I6=(float)(c3*(i6+ti)+s3*(r6+tr));
*R7=(float)(c7*(r6-tr)-s7*(i6-ti)), *I7=(float)(c7*(i6-ti)+s7*(r6-tr));
}else
{
*R1 = r0-r1, *I1 = i0-i1;
*R2 = r2-i3, *I2 = i2+r3;
*R3 = r2+i3, *I3 = i2-r3;
tr = p7*(r5-i5), ti = p7*(r5+i5);
*R4=(float)(r4+tr), *I4 = (float)(i4+ti);
*R5=(float)(r4-tr), *I5 = (float)(i4-ti);
tr = -p7*(r7+i7), ti = p7*(r7-i7);
*R6=(float)(r6+tr); *I6=(float)(i6+ti);
*R7=(float)(r6-tr); *I7=(float)(i6-ti);
}
}
}
}
if(N2pow-N8pow == 1)
{
// 基2
R1=(R0=pReal)+1, I1=(I0=pImag)+1;
for(k=0; k < Nthpo; k += 2)
{
r0 = *R0+*R1; *R1 = *R0-*R1; *R0++ = r0; R1=(++R0)+1;
i0 = *I0+*I1; *I1 = *I0-*I1; *I0++ = i0; I1=(++I0)+1;
}
}else if(N2pow-N8pow == 2)
{
// 基4
for(k=0; k < Nthpo; k+= 4)
{
R3=(R2=(R1=(R0=pReal+k)+1)+1)+1, I3=(I2=(I1=(I0=pImag+k)+1)+1)+1;
r1 = *R0+*R2, r2 = *R0-*R2, r3 = *R1+*R3, r4 = *R1-*R3;
i1 = *I0+*I2, i2 = *I0-*I2, i3 = *I1+*I3, i4 = *I1-*I3;
*R0 = r1+r3, *R1 = r1-r3, *R2 = r2-i4, *R3 = r2+i4;
*I0 = i1+i3, *I1 = i1-i3, *I2 = i2+r4, *I3 = i2-r4;
}
}
j = 0;
for (i=0; i<Nthpo; i++)
{
if (j>i)
{
r0 = pReal[j]; pReal[j] = pReal[i]; pReal[i] = r0;
i0 = pImag[j]; pImag[j] = pImag[i]; pImag[i] = i0;
}
k = Nthpo;
while( (k>>=1) > 1 && j >= k) j -= k;
j += k;
}
r0 = (float)(1.0/sqrt(double(Nthpo)));
if(TransformSign) for(i=0; i<Nthpo; i++) *pReal++ *= r0, *pImag++ *= -r0;
else for(i=0; i<Nthpo; i++) *pReal++ *= r0, *pImag++ *= r0;
}
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -