最近在项目中需要用到FFT，之前对于FFT也只是有一个模糊的印象也并不清楚他的具体物理意义，之前几次想学习都被搁置了，现在项目需要又从新学习，在此把我收获的和大家分享一下：
1- FFT简介
FFT是一种DFT的高效算法，称为快速傅立叶变换（fast Fourier transform）。傅里叶变换是时域--频域变换分析中最基本的方法之一。可以将一个信号变换到频域。有些信号在时域上很难看出什么特征，不利于分析，但是如果变换到频域之后，就很容易看出信号的特征了。这就是很多信号分析采用FFT变换的原因。FFT也可以将一个信号的频谱提取出来，常应用于频谱分析上。
2 -采样定理
采样频率必须是信号的最高频率的两倍及其以上，才能保证被采样的信号不失真。
3- FFT的物理意义
现假设我们需要对一个信号进行采样然后做FFT分析，设定其采样频率为Fs，信号的频率F，采样点数为N。那么FFT之后结果其实就是一个为N点的复数。每一个点就对应着一个频率点。该点的模值，就是该频率值下的幅度特性。
假设经过FFT之后得到的某点n，使用复数表示为a+bi，则该点的参数如下：
模值为A(n) 相位为P(n) = atan2(b, a) 频率表达式为：Fn = (n - 1) * Fs / N 幅度(n ¹ 1) = A(n) / (N / 2) 幅度(n = 1) = A(n) / N (直流分量0Hz) 分辨率 = Fs / N
|
例：
假设现在有一个信号由三个波形组成，分别是幅度为2的直流信号、幅度为3频率为50Hz相位为-30度的正弦波、幅度为1.5频率为75Hz相位为90度的正弦波。使用数学表达式表示为如下所示：
S=Adc+A1*cos(2*pi*F1*t+pi*P1/180)+A2*cos(2*pi*F2*t+pi*P2/180);
假设我们以256Hz采样频率对该信号进行采样，采样点数为256点，即N = 256。由上面总结的公式可是某点n对应的频率为 Fn = (n - 1) * 256 / 256 = n-1，我们使用数学模型模拟的信号频率分别为0Hz、50Hz、75Hz，由于采样的点是从1开始的，因此对应FFT运算之后的数据应该是1点、51点、76点，只要分析以上三点的数据即可知道对应波形的幅度、相位等信息。
1点： 512+0i
2点： -2.6195E-14 - 1.4162E-13i
3点： -2.8586E-14 - 1.1898E-13i
…
50点：-6.2076E-13 - 2.1713E-12i
51点：332.55 - 192i
52点：-1.6707E-12 - 1.5241E-12i
…
75点：-2.2199E-13 -1.0076E-12i
76点：3.4315E-12 + 192i
77点：-3.0263E-14 +7.5609E-13i
由以上数据可以分析出1点、51点、76点对应的模值分别为512、384、192，因此按照以上公式可以得出n=1时频率为0Hz点对应的为直流分量，该点的幅度为：A(1) = 512 / N = 2，51点对应的幅度为：A(51)= 384 / (N / 2)= 3，76点对应的幅度为：A(76)= 192 / (N / 2)= 1.5。
相位计算：对于直流信号来说无相位可言，51点对应的相位为：atan(-192, 332.55) = -0.5236，计算的结果为弧度，转换成角度为：180*（-0.5236）/ p = -30.0001度，76点对应的相位为：atan(192,3.4315E-12) = 1.5708，换算成角度为：90.0002，由此可见，分析后的数据和我们设定的数学模型是相符合的。Matlab仿真结果如下：
Matlab仿真代码如下：
close all; %先关闭所有图片 Adc=2; %直流分量幅度 A1=3; %频率F1信号的幅度 A2=1.5; %频率F2信号的幅度 F1=50; %信号1频率(Hz) F2=75; %信号2频率(Hz) Fs=256; %采样频率(Hz) P1=-30; %信号1相位(度) P2=90; %信号相位(度) N=256; %采样点数 t=[0:1/Fs:N/Fs]; %采样时刻
%信号 S=Adc+A1*cos(2*pi*F1*t+pi*P1/180)+A2*cos(2*pi*F2*t+pi*P2/180); %显示原始信号 plot(S); title('原始信号');
figure; Y = fft(S,N); %做FFT变换 Ayy = (abs(Y)); %取模 plot(Ayy(1:N)); %显示原始的FFT模值结果 title('FFT 模值');
figure; Ayy=Ayy/(N/2); %换算成实际的幅度 Ayy(1)=Ayy(1)/2; F=([1:N]-1)*Fs/N; %换算成实际的频率值 plot(F(1:N/2),Ayy(1:N/2)); %显示换算后的FFT模值结果 title('幅度-频率曲线图'); figure; Pyy=[1:N/2]; for i=1:N/2 Pyy(i)=phase(Y(i)); %计算相位 Pyy(i)=Pyy(i)*180/pi; %换算为角度 end; plot(F(1:N/2),Pyy(1:N/2)); %显示相位图 title('相位-频率曲线图'); |
4- FFT的算法实现
本节将使用STM32官方提供的DSP库运行FFT算法，使用DSP库的好处在于FFT的运算时间很快，在72Mhz的主频下运行256点的FFT仅仅只需要零点几毫秒，能基本满足我们的测试FFT的需求。
void Init_Single(void) { unsigned short i; float fx; for(i = 0; i < N; i++) { fx = 1500 * sin(PI2 * i * 350.0 / Fs) + 2700 * sin(PI2 * i * 8400.0 / Fs) + 4000 * sin(PI2 * i * 18725.0 / Fs); In_Array[i] = ((signed short)fx) << 16; } } void Get_Single_Message() { signed short lX,lY; float X,Y,Mag; unsigned short i; unsigned int f = 0; SEGGER_RTT_printf(0, "Num f(Hz) A X Y\n"); for(i = 0; i < N/2; i++) { lX = (Out_Array[i] << 16) >> 16; lY = (Out_Array[i] >> 16); X = N * ((float)lX) / 32768; Y = N * ((float)lY) / 32768; Mag = sqrt(X * X + Y * Y) / (N / 2); //某点的幅值=该点的模值/（N/2） if(i == 0) MagArray[i] = (unsigned long)(Mag * 32768); else { MagArray[i] = (unsigned long)(Mag * 32768); f += 175; } SEGGER_RTT_printf(0, " %d %d %d %d %d\n", i, f, MagArray[i], Mag, Y); } SEGGER_RTT_printf(0, "-------------------------------------------------------\n"); }
以上若有不严谨之处还请各位指出，谢谢！
