FFT快速傅里叶变化源码

// fft1.cpp : 此文件包含 "main" 函数。程序执行将在此处开始并结束。
//

// #include "pch.h"
#include <iostream>
#include <math.h>
#include <stdio.h>


#define N 64 //64点
#define log2N 6 //log2N=6
//复数类型
typedef struct
{
	float real;
	float img;
}complex;
complex  x[N] = 
{
{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},
{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},
{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},
{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},
{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},{4,0},{1,0},{3,0},{2,0},{5,0},{8,0},
{4,0},{1,0},{3,0},{2,0},{5,0}
}; //输入序列
//复数加法
complex add(complex a, complex b)
{
	complex c;
	c.real = a.real + b.real;
	c.img = a.img + b.img;
	return c;
}
//复数减法
complex sub(complex a, complex b)
{
	complex c;
	c.real = a.real - b.real;
	c.img = a.img - b.img;
	return c;
}
//复数乘法
complex mul(complex a, complex b)
{
	complex c;
	c.real = a.real*b.real - a.img*b.img;
	c.img = a.real*b.img + a.img*b.real;
	return c;
}
//码位倒序函数
void Reverse(void)
{
	unsigned int i, j, k;
	unsigned int t;
	complex temp;//临时交换变量
	for (i = 0;i < N;i++)//从第0个序号到第N-1个序号
	{
		k = i;//当前第i个序号
		j = 0;//存储倒序后的序号，先初始化为0
		for (t = 0;t < log2N;t++)//共移位t次，其中log2N是事先宏定义算好的
		{
			j <<= 1;
			j |= (k & 1);//j左移一位然后加上k的最低位
			k >>= 1;//k右移一位，次低位变为最低位
		}
		if (j > i)//如果倒序后大于原序数，就将两个存储单元进行交换(判断j>i是为了防止重复交换)
		{
			temp = x[i];
			x[i] = x[j];
			x[j] = temp;
		}
	}
}
complex  WN[N] =
{{1.00000,0.00000},{0.99518,-0.09802},{0.98079,-0.19509},{0.95694,-0.29028},
  {0.92388,-0.38268},{0.88192,-0.47140},{0.83147,-0.55557},{0.77301,-0.63439},
  {0.70711,-0.70711},{0.63439,-0.77301},{0.55557,-0.83147},{0.47140,-0.88192},
  {0.38268,-0.92388},{0.29028,-0.95694},{0.19509,-0.98079},{0.09802,-0.99518},
  {0.00000,-1.00000},{-0.09802,-0.99518},{-0.19509,-0.98079},{-0.29028,-0.95694},
  {-0.38268,-0.92388},{-0.47140,-0.88192},{-0.55557,-0.83147},{-0.63439,-0.77301},
  {-0.70711,-0.70711},{-0.77301,-0.63439},{-0.83147,-0.55557},{-0.88192,-0.47140},
  {-0.92388,-0.38268},{-0.95694,-0.29028},{-0.98079,-0.19509},{-0.99518,-0.09802},
  {-1.00000,0.00000},{-0.99518,0.09802},{-0.98079,0.19509},{-0.95694,0.29028},
  {-0.92388,0.38268},{-0.88192,0.47140},{-0.83147,0.55557},{-0.77301,0.63439},
  {-0.70711,0.70711},{-0.63439,0.77301},{-0.55557,0.83147},{-0.47140,0.88192},
  {-0.38268,0.92388},{-0.29028,0.95694},{-0.19509,0.98079},{-0.09802,0.99518},
  {0.00000,1.00000},{0.09802,0.99518},{0.19509,0.98079},{0.29028,0.95694},
  {0.38268,0.92388},{0.47140,0.88192},{0.55557,0.83147},{0.63439,0.77301},
  {0.70711,0.70711},{0.77301,0.63439},{0.83147,0.55557},{0.88192,0.47140},
  {0.92388,0.38268},{0.95694,0.29028},{0.98079,0.19509},{0.99518,0.09802}
};
void FFT(void)
{
	unsigned int i, j, k, l;
	complex top, bottom, xW;
	Reverse(); //码位倒序
	for (i = 0;i < log2N;i++)   /*共log2N级*/
	{ //一级蝶形运算
		l = 1 << i;//l等于2的i次方
		for (j = 0;j < N;j += 2 * l)  /*每L个蝶形是一组，每级有N/2L组*/
		{ //一组蝶形运算
			for (k = 0;k < l;k++)   /*每组有L个*/
			{ //一个蝶形运算
				xW = mul(x[j + k + l], WN[N / (2 * l)*k]); //碟间距为l
				top = add(x[j + k], xW); //每组的第k个蝶形
				bottom = sub(x[j + k], xW);
				x[j + k] = top;
				x[j + k + l] = bottom;
			}
		}
	}
}

int main()
{
	int i;

	FFT();
	for (i = 0;i < 64;i++)
	{
		printf("i is :%d real %f image %f\n",i,x[i].real,x[i].img);
	}




    std::cout << "成功输出\n"; 
}