fft842.c

基二、基四、基八 FFT变化源代码

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;
}