matlabfft变换后的相位精度问题_FFT
最近在项⽬中需要⽤到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+0i2点： -2.6195E-14 - 1.4162E-13i3点： -2.8586E-14 - 1.1898E-13i
…50点：-6.2076E-13 - 2.1713E-12i51点：332.55 - 192i52点：-1.6707E-12 - 1.5241E-12i
…75点：-2.2199E-13 -1.0076E-12i76点：3.4315E-12 + 192i77点：-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傅里叶变换的幅度谱和相位谱
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的需求。
1. void Init_Single(void)
2. {
3. unsigned short i;
4. float fx;
5.
6. for(i = 0; i
7. {
8. fx = 1500 * sin(PI2 * i * 350.0 / Fs) +
9. 2700 * sin(PI2 * i * 8400.0 / Fs) +
10. 4000 * sin(PI2 * i * 18725.0 / Fs);
11. In_Array[i] = ((signed short)fx) <
12. }
13.}
14.
15.void Get_Single_Message()
16.{
17. signed short lX,lY;
18. float X,Y,Mag;
19. unsigned short i;
20. unsigned int f = 0;
21.
22. SEGGER_RTT_printf(0, "Num f(Hz) A X Y\n");
23.
24. for(i = 0; i
25. {
26. lX = (Out_Array[i] <> 16;
27. lY = (Out_Array[i] >> 16);
28. X = N * ((float)lX) / 32768;
29. Y = N * ((float)lY) / 32768;
30. Mag = sqrt(X * X + Y * Y) / (N / 2); //某点的幅值=该点的模值/(N/2)
31. if(i == 0)
32. MagArray[i] = (unsigned long)(Mag * 32768);
33. else
34. {
35. MagArray[i] = (unsigned long)(Mag * 32768);
36. f += 175;
37. }
38.
39. SEGGER_RTT_printf(0, " %d %d %d %d %d\n", i, f, MagArray[i], Mag, Y);
40. }
41.
42. SEGGER_RTT_printf(0, "-------------------------------------------------------\n");
43.}
以上若有不严谨之处还请各位指出，谢谢！
版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。
发表评论