C语言版FFT简单测试
本次我们来自己封装一个FFT函数，进行简单的测试。
fft.c
#include "math.h"
#include "fft.h"
//精度0.0001弧度
//复数的交换
void conjugate_complex(int n,complex in[],complex out[])
{
int i = 0;
for(i=0;i<n;i++)
{
out[i].imag = -in[i].imag;
out[i].real = in[i].real;
}
}
//求所有复数的模
void c_abs(complex f[],float out[],int n)
{
int i = 0;
float t;
for(i=0;i<n;i++)
{
t = f[i].real * f[i].real + f[i].imag * f[i].imag;
out[i] = sqrt(t);
}
}
//求复数的和
void c_plus(complex a,complex b,complex *c)
{
c->real = a.real + b.real;
c->imag = a.imag + b.imag;
}
//求复数的差
void c_sub(complex a,complex b,complex *c)
{
c->real = a.real - b.real;
c->imag = a.imag - b.imag;
}
//求复数的积
void c_mul(complex a,complex b,complex *c)
{
c->real = a.real * b.real - a.imag * b.imag;
c->imag = a.real * b.imag + a.imag * b.real;
}
//求复数的商
void c_div(complex a,complex b,complex *c)
{
c->real = (a.real * b.real + a.imag * b.imag)/(b.real * b.real +b.imag * b.imag);
c->imag = (a.imag * b.real - a.real * b.imag)/(b.real * b.real +b.imag * b.imag);
}
#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
void Wn_i(int n,int i,complex *Wn,char flag)
{
Wn->real = cos(2*PI*i/n);
if(flag == 1)
Wn->imag = -sin(2*PI*i/n);
else if(flag == 0)
Wn->imag = -sin(2*PI*i/n);
}
//傅里叶变化
void fft(int N,complex f[])
{
complex t,wn;//中间变量
int i,j,k,m,n,l,r,M;
int la,lb,lc;
/*----计算分解的级数M=log2(N)----*/
for(i=N,M=1;(i=i/2)!=1;M++);
/*----按照倒位序重新排列原信号----*/
for(i=1,j=N/2;i<=N-2;i++)
{
if(i<j)
{
t=f[j];
f[j]=f[i];
f[i]=t;
}
k=N/2;
while(k<=j)
{
j=j-k;
k=k/2;
}
j=j+k;
}
/*----FFT算法----*/
for(m=1;m<=M;m++)
{
la=pow(2,m); //la=2^m代表第m级每个分组所含节点数
lb=la/2; //lb代表第m级每个分组所含碟形单元数
//同时它也表示每个碟形单元上下节点之间的距离
/*----碟形运算----*/
for(l=1;l<=lb;l++)
{
r=(l-1)*pow(2,M-m);
for(n=l-1;n<N-1;n=n+la) //遍历每个分组，分组总数为N/la
{
lc=n+lb; //n,lc分别代表一个碟形单元的上、下节点编号
Wn_i(N,r,&wn,1);//wn=Wnr
c_mul(f[lc],wn,&t);//t = f[lc] * wn复数运算
c_sub(f[n],t,&(f[lc]));//f[lc] = f[n] - f[lc] * Wnr
c_plus(f[n],t,&(f[n]));//f[n] = f[n] + f[lc] * Wnr
}
}
}
}
//傅里叶逆变换
void ifft(int N,complex f[])
{
int i=0;
conjugate_complex(N,f,f);
fft(N,f);
conjugate_complex(N,f,f);
for(i=0;i<N;i++)
{
f[i].imag = (f[i].imag)/N;
f[i].real = (f[i].real)/N;
}
}
fft.h
#ifndef __FFT_H__
#define __FFT_H__
typedef struct complex //复数类型
{
float real; //实部
float imag; //虚部
}complex;
#define PI 3.1415926535897932384626433832795028841971
///
void conjugate_complex(int n,complex in[],complex out[]);
void c_plus(complex a,complex b,complex *c);//复数加
void c_mul(complex a,complex b,complex *c) ;//复数乘
void c_sub(complex a,complex b,complex *c); //复数减法
void c_div(complex a,complex b,complex *c); //复数除法
void fft(int N,complex f[]);//傅立叶变换 输出也存在数组f中
void ifft(int N,complex f[]); // 傅里叶逆变换
void c_abs(complex f[],float out[],int n);//复数数组取模
#endif
测试代码部分
#include "sys.h"
#include "delay.h"
#include "usart.h"
#include "led.h"
#include "math.h"
#include "fft.h"
#define N 256 //采样点数
#define Fs 44800 //采样频率
#define F 175 //分辨率
#define P1 30 //测试相位
#define P2 60
#define P3 90
//FFT测试数据集 输入数组
complex FFT_256PointIn[N];
//FFT测试数据集 输出数组
float FFT_256PointOut[N/2];
//填入数组
void InitBufInArray()
{
unsigned short i;
for(i=0; i<N; i++)
{
FFT_256PointIn[i].real = 1500 * sin(2*PI * i * 350.0 / Fs+(PI*P1/180))
+2700 * sin(2*PI * i * 8400.0 / Fs+(PI*P2/180))
+4000 * sin(2*PI * i * 18725.0 / Fs+(PI*P3/180));
FFT_256PointIn[i].imag = 0;
}
}
/******************************************************************
函数名称:GetPowerMag()
函数功能:计算各次谐波幅值
参数说明:
备 注:先将FFT_256PointIn分解成实部(X)和虚部(Y)，
然后计算幅值:(sqrt(X*X+Y*Y)*2/N
然后计算相位:atan2(Y/X)
作 者:土耳其冰激凌
*******************************************************************/
void GetPowerMag()
{
unsigned short i;
float X,Y,P,Mag;
c_abs(FFT_256PointIn,FFT_256PointOut,N/2);
for(i=0; i<N/2; i++)
{
X = FFT_256PointIn[i].real/N; //计算实部
Y = FFT_256PointIn[i].imag/N; //计算虚部
Mag = FFT_256PointOut[i]*2/N; //计算幅值
P = atan2(Y,X)*180/PI; //计算相位
printf("%d ",i);
printf("%d ",F*i);
printf("%f ",Mag);
printf("%f ",P);
printf("%f ",X);
printf("%f \r\n",Y);
}
}
int main(void)
{
delay_init(); //延时函数初始化
LED_Init(); //初始化与LED连接的硬件接口
uart_init(9600); //初始化串口 9600波特率
printf("这是一个FFT 测试实验\r\n");
InitBufInArray();
fft(N,FFT_256PointIn);
printf("点数 频率 幅值 实部 虚部\n");
GetPowerMag();
while(1)
{
LED0=0;
LED1=1;
delay_ms(300); //延时300ms
LED0=1;
LED1=0;
delay_ms(300); //延时300ms
}
}
输出打印结果：
点数 频率 幅值 相位 实部 虚部
0 0 0.000031 0.000000 0.000015 0.000000
1 175 0.000025 36.658634 0.000010 0.000008
2 350 1500.000000 -60.000000 374.999969 -649.519043
3 525 0.000020 52.221466 0.000006 0.000008
4 700 0.000014 -28.152960 0.000006 -0.000003
5 875 0.000060 -17.167744 0.000029 -0.000009
6 1050 0.000012 -4.051142 0.000006 -0.000000
7 1225 0.000028 27.618505 0.000013 0.000007
8 1400 0.000036 56.088974 0.000010 0.000015
9 1575 0.000013 85.033867 0.000001 0.000007
10 1750 0.000019 -126.592506 -0.000006 -0.000008
11 1925 0.000099 155.458847 -0.000045 0.000021
12 2100 0.000021 -158.704437 -0.000010 -0.000004
13 2275 0.000018 -89.424667 0.000000 -0.000009
14 2450 0.000015 -87.367409 0.000000 -0.000007
15 2625 0.000019 -84.178368 0.000001 -0.000010
16 2800 0.000086 112.500000 -0.000017 0.000040
17 2975 0.000008 -30.908798 0.000003 -0.000002
18 3150 0.000028 -93.026222 -0.000001 -0.000014
19 3325 0.000032 -100.869904 -0.000003 -0.000016
20 3500 0.000028 -35.795570 0.000012 -0.000008
21 3675 0.000248 -9.752425 0.000122 -0.000021
22 3850 0.000036 -59.215248 0.000009 -0.000015
23 4025 0.000010 -76.748405 0.000001 -0.000005
24 4200 0.000050 154.164047 -0.000023 0.000011
25 4375 0.000015 -60.317295 0.000004 -0.000006
26 4550 0.000014 108.588715 -0.000002 0.000006
27 4725 0.000029 98.175682 -0.000002 0.000014
28 4900 0.000016 -124.238319 -0.000004 -0.000007
29 5075 0.000048 108.763763 -0.000008 0.000023
30 5250 0.000047 -21.094650 0.000022 -0.000008
31 5425 0.000009 -38.687473 0.000003 -0.000003
32 5600 0.000113 -87.260597 0.000003 -0.000056
33 5775 0.000016 -63.487057 0.000004 -0.000007
34 5950 0.000022 7.387201 0.000011 0.000001
35 6125 0.000027 63.663502 0.000006 0.000012
36 6300 0.000032 62.835876 0.000007 0.000014
37 6475 0.000076 -38.793385 0.000030 -0.000024
38 6650 0.000009 16.024696 0.000004 0.000001
39 6825 0.000012 -110.747261 -0.000002 -0.000006
40 7000 0.000028 82.453056 0.000002 0.000014
41 7175 0.000016 -27.415680 0.000007 -0.000004
42 7350 0.000027 160.582352 -0.000013 0.000004
43 7525 0.000087 -87.681641 0.000002 -0.000044
44 7700 0.000022 -135.357925 -0.000008 -0.000008
45 7875 0.000014 -169.563293 -0.000007 -0.000001
46 8050 0.000009 1.227116 0.000004 0.000000
47 8225 0.000016 42.631676 0.000006 0.000006
48 8400 2700.000000 -30.000000 1169.134277 -675.000000
49 8575 0.000008 175.232651 -0.000004 0.000000
50 8750 0.000017 147.675629 -0.000007 0.000005
51 8925 0.000056 179.124176 -0.000028 0.000000
52 9100 0.000021 126.577728 -0.000006 0.000008
53 9275 0.000040 28.367682 0.000018 0.000010
54 9450 0.000011 -37.718578 0.000004 -0.000003
55 9625 0.000010 131.671509 -0.000003 0.000004
56 9800 0.000026 172.892471 -0.000013 0.000002
57 9975 0.000018 -143.289139 -0.000007 -0.000005
58 10150 0.000012 -115.207047 -0.000003 -0.000006
59 10325 0.000008 -87.947166 0.000000 -0.000004
60 10500 0.000016 -127.038078 -0.000005 -0.000006
61 10675 0.000024 35.620342 0.000010 0.000007
62 10850 0.000000 0.000000 0.000000 0.000000
63 11025 0.000030 167.302933 -0.000015 0.000003
64 11200 0.000087 -37.874985 0.000034 -0.000027
65 11375 0.000009 -56.846428 0.000002 -0.000004
66 11550 0.000000 0.000000 0.000000 0.000000
67 11725 0.000016 65.517853 0.000003 0.000007
68 11900 0.000020 -62.873825 0.000005 -0.000009
69 12075 0.000049 146.105896 -0.000021 0.000014
70 12250 0.000025 144.522079 -0.000010 0.000007
71 12425 0.000022 83.036545 0.000001 0.000011
72 12600 0.000024 124.945473 -0.000007 0.000010
73 12775 0.000017 150.301559 -0.000008 0.000004
74 12950 0.000006 73.393593 0.000001 0.000003
75 13125 0.000064 101.329811 -0.000006 0.000031
76 13300 0.000022 21.395346 0.000010 0.000004
77 13475 0.000044 -22.719753 0.000020 -0.000008
78 13650 0.000002 67.932175 0.000000 0.000001
79 13825 0.000015 36.952290 0.000006 0.000004
80 14000 0.000061 -90.000000 0.000000 -0.000031
81 14175 0.000024 4.587837 0.000012 0.000001
82 14350 0.000027 -164.107590 -0.000013 -0.000004
83 14525 0.000026 94.922646 -0.000001 0.000013
84 14700 0.000017 113.708290 -0.000003 0.000008
85 14875 0.000049 -94.115646 -0.000002 -0.000025
86 15050 0.000020 3.759978 0.000010 0.000001
87 15225 0.000019 -4.117675 0.000010 -0.000001
88 15400 0.000017 -42.499950 0.000006 -0.000006
89 15575 0.000017 -79.784363 0.000002 -0.000009
90 15750 0.000016 -98.852882 -0.000001 -0.000008
91 15925 0.000054 139.931839 -0.000021 0.000017
92 16100 0.000026 165.618698 -0.000013 0.000003
93 16275 0.000015 -164.399490 -0.000007 -0.000002
94 16450 0.000016 -51.657200 0.000005 -0.000006
95 16625 0.000035 167.911438 -0.000017 0.000004
96 16800 0.000008 -132.260605 -0.000003 -0.000003
97 16975 0.000021 94.842651 -0.000001 0.000011
98 17150 0.000035 -114.642807 -0.000007 -0.000016
99 17325 0.000037 34.428101 0.000015 0.000011
100 17500 0.000016 -93.046417 -0.000000 -0.000008
101 17675 0.000022 32.025196 0.000009 0.000006
102 17850 0.000022 159.313171 -0.000010 0.000004
103 18025 0.000012 137.351868 -0.000005 0.000004
104 18200 0.000034 -31.864140 0.000014 -0.000009
105 18375 0.000028 -157.610428 -0.000013 -0.000005
106 18550 0.000020 56.170975 0.000006 0.000008
107 18725 3999.999756 -0.000001 1999.999878 -0.000040
108 18900 0.000024 70.995369 0.000004 0.000011
109 19075 0.000075 44.852394 0.000026 0.000026
110 19250 0.000026 -126.419708 -0.000008 -0.000010
111 19425 0.000031 -90.581345 -0.000000 -0.000016
112 19600 0.000086 67.500000 0.000017 0.000040
113 19775 0.000032 117.184265 -0.000007 0.000014
114 19950 0.000016 -152.757553 -0.000007 -0.000004
115 20125 0.000047 -12.603742 0.000023 -0.000005
116 20300 0.000009 -68.715248 0.000002 -0.000004
117 20475 0.000067 144.246552 -0.000027 0.000020
118 20650 0.000014 8.900421 0.000007 0.000001
119 20825 0.000013 150.707443 -0.000006 0.000003
120 21000 0.000051 67.017677 0.000010 0.000023
121 21175 0.000016 8.680163 0.000008 0.000001
122 21350 0.000013 -115.376968 -0.000003 -0.000006
123 21525 0.000073 7.219605 0.000036 0.000005
124 21700 0.000011 45.488987 0.000004 0.000004
125 21875 0.000051 79.440254 0.000005 0.000025
126 22050 0.000061 0.000000 0.000031 0.000000
127 22225 0.000012 -156.431137 -0.000006 -0.000002
本次所测得结果更加接近matlab 所测量的结果，而且还加入了相位角，后续的FFT音乐频谱制作，我会使用C语言的FFT来进行操作。
9711
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
