fft源码，VS2010上可编译运行

最新推荐文章于 2026-05-01 10:28:08 发布

原创最新推荐文章于 2026-05-01 10:28:08 发布 · 1.9k 阅读

5 ·

本内容遵循CC 4.0 BY-SA版权协议

标签

#c语言 #fft #傅立叶变换快速傅立变换蝶形

源码专栏收录该内容

2 篇文章

订阅专栏

本文介绍了如何在Visual Studio 2010环境下使用C语言编写并编译运行快速傅立叶变换(FFT)程序，重点讲解了关键的蝶形运算流程。

#include<math.h>
#include<stdio.h>
#include <memory>

#define FFT_COUNT	256
#define PI	3.1415926f
int FFT_R = 0;

double val[FFT_COUNT][2];
double val2[FFT_COUNT][2];
double fft_w[FFT_COUNT][2];
double fft_x1[FFT_COUNT][2];
int pos[FFT_COUNT];
double anm[FFT_COUNT];

//fft 如果isr==1表示傅氏变换时域->频域，否则表示反变换频域->时域
void fft2(double (*val_in)[2], double (*X1)[2], int isr)
{
	double  rcos = 0.0f, isin  = 0.0f;	//中间临时变量
	int mmax = 0;
	int n = 0, i=0, k=0, j=0;
	
	//重新排序输入量
	for(i = 0; i<FFT_COUNT; i++)
	{
		X1[pos[i]][0] = val_in[i][0];
		X1[pos[i]][1] = val_in[i][1] * isr; //0;//
	}

	//第一层，表示几个蝶形
	for(k=0; k<FFT_R; k++ )
	{
		mmax = 1 << ( k + 1 ); //2,4,8
		for( j = 0; j <(mmax>>1) ; j++)		//	
		{
			int wk = j * ( 1 << ( FFT_R - k - 1) );
			for( i = j; i < FFT_COUNT; i+= mmax )
			{
				n = i + ( 1 << k );
				rcos = X1[n][0] * fft_w[wk][0] - X1[n][1] * fft_w[wk][1];
				isin = X1[n][0] * fft_w[wk][1] + X1[n][1] * fft_w[wk][0];

				X1[n][0] = X1[i][0] - rcos;
				X1[n][1] = X1[i][1] - isin;
				X1[i][0] += rcos;
				X1[i][1] += isin;
			}
		}
	}
	if(isr==1)
	{
	}
	else
	{
		for(i=0; i<FFT_COUNT; i++)
		{
			X1[i][0] = X1[i][0] / FFT_COUNT;
			X1[i][1] = X1[i][1] / (-FFT_COUNT);
		}
	}
}

void main()
{
	int i = 0;
	int j = 0;
	int k = 0;
	
	FFT_R = 0;
	while( (1 << FFT_R ) <  FFT_COUNT ) FFT_R++;	//变换次数是2的多少次方

	// 倒序排序 如二进制的 011 变成 110
	for(i=0; i<FFT_COUNT; i++)
	{
		k = 0;
		for(j=0; j<FFT_R; j++)
		{
			if( ( ( i >> j ) & 1 ) == 1 )  k |= ( 1 << ( FFT_R - j-1 ) );
		}
		pos[i] = k;
	}

	//生成变换因子w
	for(i = 0; i < FFT_COUNT / 2; i++)
	{
		fft_w[i][0] = cos( -2.0f * PI * i / FFT_COUNT );
		fft_w[i][1] = sin( -2.0f * PI * i / FFT_COUNT );
	}

	//生成需要变换的模拟向量
	for(i = 0; i<FFT_COUNT; i++)
	{
		val[i][0] = sqrt(2) * cos(2 * PI * i / FFT_COUNT ) + 0.0*sqrt(2) * cos(2 * PI * i *50.0f/49.0f/ FFT_COUNT );
		val[i][1] = 0;
	}
	
	fft2(val,fft_x1, 1 );
	fft2(fft_x1, val2, -1 );		//反变换

	char bufff[10];
	float anlog = 0.0f;
	for(i=0; i<FFT_COUNT>>1; i++)
	{
		fft_x1[i][0] /= FFT_COUNT/1.41421356f;
		fft_x1[i][1] /= (-FFT_COUNT/1.41421356f);
		
		//幅值
		anm[i] = sqrt(fft_x1[i][0] * fft_x1[i][0] + fft_x1[i][1] * fft_x1[i][1]);
		anlog += fft_x1[i][0] * fft_x1[i][0] + fft_x1[i][1] * fft_x1[i][1];
		memset(bufff, 0, 10);
		sprintf(bufff, "  %3.*f",  3, anm[i]);
		for(j=9; j>0; j--)
		{
			if( bufff[j]=='0' ) bufff[j] = 0;
			if( bufff[j]=='.')
			{
				bufff[j]=0; break;
			}
			if( bufff[j]!=0) break;
		}
		 printf(bufff);
		if( (i+1) % 16==0 ) printf("\r\n");
	}
	printf("\r\n");
	printf("%.4f", (float)sqrt(anlog) );

}