使用C语言实现FFT算法（快速傅里叶变换）
文章目录(一)FFT的基本原理(二)FFT代码(三)使用典型函数进行测试(一)FFT的基本原理FFT的总思路第一步，将长序列DFT分解成短序列的DFT（本来是N点的DFT，被分解成两个 N/2 的DFT）第二步，分别将 N/2 的DFT分解成 N/4 的DFT第三步，继续分。。。最后被分成 2点 的DFT注意事项只要开始的序列长度是2的正整数次方，就可...
·
文章目录
(一)FFT的基本原理
(二)FFT代码
(三)使用典型函数进行测试
(一)FFT的基本原理
-
FFT的总思路
第一步，将长序列DFT分解成短序列的DFT（本来是N点的DFT，被分解成两个 N/2 的DFT） 第二步，分别将 N/2 的DFT分解成 N/4 的DFT 第三步，继续分。。。最后被分成 2点 的DFT -
注意事项
只要开始的序列长度是2的正整数次方，就可以最终被分解成2点的DFT 分解原则是：对输入进行偶奇分解，对输出进行前后分解（关于偶奇分解可以看上一篇雷德算法） -
前后分解
-
蝶形运算（可以体现出偶奇分解和前后分解）
(二)FFT代码
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define N 1024
#define PI 3.1415
/*复数结构体类型*/
typedef struct{
double realPart;
double imaginaryPart;
}complexNumber;
complexNumber x[N], W[N/2]; //序列数组以及旋转因子数组
int size_x=0; //输入序列的长度，只能为2的次幂
void fastFourierOperation();//快速傅里叶变换算法
void initRotationFactor(); //生成旋转因子数组
void changePosition(); //偶奇变换算法
void outputArray(); //遍历输出数组
void add(complexNumber ,complexNumber ,complexNumber*); //复数加法
void mul(complexNumber ,complexNumber ,complexNumber*); //复数乘法
void sub(complexNumber ,complexNumber ,complexNumber*); //复数减法
int main()
{
printf("请输入序列的长度："); //输入序列的大小
scanf("%d",&size_x);
printf("\n请分别输入序列的实部和虚部\n"); //输入序列的实部和虚部
for(int i=0;i<size_x;i++)
{
printf("请输入第%d个序列的实部：",i);
scanf("%lf",&x[i].realPart);
printf("请输入第%d个序列的虚部：",i);
scanf("%lf",&x[i].imaginaryPart);
}
printf("\n输出原序列：\n");
outputArray();
initRotationFactor(); //生成旋转因子数组
fastFourierOperation();//调用快速傅里叶变换
printf("\n输出FFT后的结果：\n");
outputArray(); //遍历输出数组
return 0;
}
/*快速傅里叶变换*/
void fastFourierOperation()
{
int l=0;
complexNumber up,down,product;
changePosition(); //偶奇变换算法
for(int i=0;i<size_x/2;i++) //一级蝶形运算
{
l=1<<i;
for(int j=0;j<size_x;j+= 2*l ) //一组蝶形运算
{
for(int k=0;k<l;k++) //一个蝶形运算
{
mul(x[j+k+l],W[size_x*k/2/l],&product);
add(x[j+k],product,&up);
sub(x[j+k],product,&down);
x[j+k]=up;
x[j+k+l]=down;
}
}
}
}
/*生成旋转因子数组*/
void initRotationFactor()
{
for(int i=0;i<size_x/2;i++)
{
W[i].realPart=cos(2*PI/size_x*i);//用欧拉公式计算旋转因子
W[i].imaginaryPart=-1*sin(2*PI/size_x*i);
}
}
/*偶奇变换算法*/
void changePosition()
{
int j=0,k;//待会儿进行运算时需要用到的变量
complexNumber temp;
for (int i=0;i<size_x-1;i++)
{
if(i<j)
{
//若倒序数大于顺序数，进行变址（换位）；否则不变，防止重复交换，变回原数
temp=x[i];
x[i]=x[j];
x[j]=temp;
}
k=size_x/2; //判断j的最高位是否为0
while(j>=k)
{ //最高位为1
j=j-k;
k=k/2;
}
j=j+k; //最高位为0
}
printf("\n输出倒序后的结果：\n");
outputArray();
}
/*遍历输出数组*/
void outputArray()
{
for(int i=0;i<size_x;i++)
{
if(x[i].imaginaryPart<0){
printf("%.4f-j%.4f\n",x[i].realPart,abs(x[i].imaginaryPart));
}
else{
printf("%.4f+j%.4f\n",x[i].realPart,x[i].imaginaryPart);
}
}
}
/*复数加法的定义*/
void add(complexNumber a,complexNumber b,complexNumber *c)
{
c->realPart=a.realPart+b.realPart;
c->imaginaryPart=a.imaginaryPart+b.imaginaryPart;
}
/*复数乘法的定义*/
void mul(complexNumber a,complexNumber b,complexNumber *c)
{
c->realPart=a.realPart*b.realPart - a.imaginaryPart*b.imaginaryPart;
c->imaginaryPart=a.realPart*b.imaginaryPart + a.imaginaryPart*b.realPart;
}
/*复数减法的定义*/
void sub(complexNumber a,complexNumber b,complexNumber *c)
{
c->realPart=a.realPart-b.realPart;
c->imaginaryPart=a.imaginaryPart-b.imaginaryPart;
}
(三)使用典型函数进行测试
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define N 8
/*复数结构体类型*/
typedef struct{
double realPart;
double imaginaryPart;
}complexNumber;
complexNumber x[N], W[N/2]; //序列数组以及旋转因子数组
int size_x=N;
double PI=atan(1)*4; //圆周率
void fastFourierOperation();//快速傅里叶变换算法
void initRotationFactor(); //生成旋转因子数组
void changePosition(); //偶奇变换算法
void outputArray(); //遍历输出数组
void add(complexNumber ,complexNumber ,complexNumber*); //复数加法
void mul(complexNumber ,complexNumber ,complexNumber*); //复数乘法
void sub(complexNumber ,complexNumber ,complexNumber*); //复数减法
/*长度为1024的冲激序列 */
void init1()
{
printf("Example：度为1024的冲激序列\n");
x[0].realPart=1;
x[0].imaginaryPart=0;
for(int i=1;i<N;i++)
{
x[i].realPart=0;
x[i].imaginaryPart=0;
}
}
/*长度为1024的矩形序列 */
void init2()
{
printf("Example：长度为1024的矩形序列\n");
for(int i=0;i<N;i++)
{
x[i].realPart=1;
x[i].imaginaryPart=0;
}
}
/*长度为1024的sin[2πk/1024]序列 */
void init3()
{
printf("Example：长度为1024的sin[2πk/1024]序列\n");
for(int i=0;i<N;i++)
{
x[i].realPart=sin(2*PI*i/N);
x[i].imaginaryPart=0;
}
}
/*长度为1024的cos[2πk/1024]序列 */
void init4()
{
printf("Example：长度为1024的cos[2πk/1024]序列\n");
for(int i=0;i<N;i++)
{
x[i].realPart=cos(2*PI*i/N);
x[i].imaginaryPart=0;
}
}
int main()
{
init3(); //调用不同的初始化函数，使用不同的例程进行测试
printf("=============================================================================================================================================================================================\n");
printf("输出原序列：\n");
outputArray();
initRotationFactor(); //生成旋转因子数组
fastFourierOperation();//调用快速傅里叶变换
printf("\n\n=============================================================================================================================================================================================\n");
printf("输出FFT后的结果：\n");
outputArray(); //遍历输出数组
return 0;
}
/*快速傅里叶变换*/
void fastFourierOperation()
{
int l=0;
complexNumber up,down,product;
changePosition(); //偶奇变换算法
for(int i=0;i<log(size_x)/log(2);i++) //一级蝶形运算
{
l=1<<i;
for(int j=0;j<size_x;j+= 2*l ) //一组蝶形运算
{
for(int k=0;k<l;k++) //一个蝶形运算
{
mul(x[j+k+l],W[size_x*k/2/l],&product);
add(x[j+k],product,&up);
sub(x[j+k],product,&down);
x[j+k]=up;
x[j+k+l]=down;
}
}
}
}
/*生成旋转因子数组*/
void initRotationFactor()
{
for(int i=0;i<size_x/2;i++)
{
W[i].realPart=cos(2*PI/size_x*i);//用欧拉公式计算旋转因子
W[i].imaginaryPart=-1*sin(2*PI/size_x*i);
}
}
/*偶奇变换算法*/
void changePosition()
{
int j=0,k;//待会儿进行运算时需要用到的变量
complexNumber temp;
for (int i=0;i<size_x-1;i++)
{
if(i<j)
{
//若倒序数大于顺序数，进行变址（换位）；否则不变，防止重复交换，变回原数
temp=x[i];
x[i]=x[j];
x[j]=temp;
}
k=size_x/2; //判断j的最高位是否为0
while(j>=k)
{ //最高位为1
j=j-k;
k=k/2;
}
j=j+k; //最高位为0
}
printf("\n\n=============================================================================================================================================================================================\n");
printf("输出倒序后的结果：\n");
outputArray();
}
/*遍历输出数组*/
void outputArray()
{
int num=0;
for(int i=0;i<size_x;i++)
{
if(x[i].imaginaryPart<0)
{
printf("%.4f-j%.4f ",x[i].realPart,fabs(x[i].imaginaryPart));
}
else
{
printf("%.4f+j%.4f ",x[i].realPart,x[i].imaginaryPart);
}
num++;
if(num==10){
printf("\n");
num=0;
}
}
}
/*复数加法的定义*/
void add(complexNumber a,complexNumber b,complexNumber *c)
{
c->realPart=a.realPart+b.realPart;
c->imaginaryPart=a.imaginaryPart+b.imaginaryPart;
}
/*复数乘法的定义*/
void mul(complexNumber a,complexNumber b,complexNumber *c)
{
c->realPart=a.realPart*b.realPart - a.imaginaryPart*b.imaginaryPart;
c->imaginaryPart=a.realPart*b.imaginaryPart + a.imaginaryPart*b.realPart;
}
/*复数减法的定义*/
void sub(complexNumber a,complexNumber b,complexNumber *c)
{
c->realPart=a.realPart-b.realPart;
c->imaginaryPart=a.imaginaryPart-b.imaginaryPart;
}
部分运行效果如下：
开放原子开发者工作坊旨在鼓励更多人参与开源活动，与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动，如meetup、训练营等，主打技术交流，干货满满，真诚地邀请各位开发者共同参与！
更多推荐
41
0- 0
已为社区贡献1条内容
所有评论(0)