前言
做OFDM通信少不了频谱分析,基带信号DA后的频谱,以及基带数字上变频后的DA信号都要频谱分析。我觉得其实做任何工程都是这样,先规定实施方案,然后仿真成功,再实际开发,不过也可以一边开发,一边仿真,开发结果要与仿真预期结果一致。
所以分析与仿真工具MATLAB就很重要了,既可以仿真,又可以通过示波器或其他方法把实际信号采下来分析。
matlab使用FFT函数分析信号频谱
一般我使用的FFT分析频谱流程如下:
%% 两个频率分别为15HZ 和 20HZ 的正弦信号[1]
Fs=50;%采样率
f1=15;
f2=20;
t = 0:1/Fs:10-1/Fs; %500个点
x = sin(2*pi*f1*t) + sin(2*pi*f2*t);
figure;
plot(t,x);
y = fft(x);
%将横坐标转化,显示为频率f= n*(fs/N)
f = (0:length(y)-1)*Fs/length(y);
figure;
plot(f,abs(y));
title('Magnitude');
%该变换还会生成尖峰的镜像副本,该副本对应于信号的负频率。
%为了更好地以可视化方式呈现周期性,可以使用 fftshift 函数对变换执行以零为中心的循环平移。
n = length(x);
fshift = (-n/2:n/2-1)*(Fs/n);
yshift = fftshift(y);
figure;
plot(fshift,abs(yshift));
图1时域图像 图2直接FFT的结果 图3 FFT平移后的结果
其中有3个注意的点:
1.FFT的结果看的是频谱,所以怎么把横坐标的值从原来的FFT点数0:N-1转换为频率值呢?
首先要引出频谱分辨率的概念,即分辨两个不同频率信号的最小间隔,FFT结果相邻点间的间隔
。因为N点FFT对应采样率为fs的序列,其频率分辨率为
,其中Ts为采样周期,T为整个序列的时间长度。有关频率分辨率的就不多说了。所以我们横坐标转换为:f = (0:length(y)-1)*Fs/length(y);
2.直接FFT的结果里怎么又多余的信号频率(镜像频率)图2?
DFT具有对称性,因为其是周期序列DFS在一个周期内的点,时域序列是有限长实序列,DFT的结果的实部周期偶对称,虚部周期奇对称,也就是模值周期偶对称,相位周期奇对称。其实从奈奎斯特定律也可以看出,fs>=2f,fs的采样率最多也就显示fs/2的真实频率(感性理解哈哈)。
所以程序处理方式就是周期延拓后取-N/2:N/2-1.用到函数fftshift(),结果如图3.如注释所述:
%该变换还会生成尖峰的镜像副本,该副本对应于信号的负频率。
%为了更好地以可视化方式呈现周期性,可以使用 fftshift 函数对变换执行以零为中心的循环平移。
其实这和设计数字滤波器IIR与FIR也一样,采样率为fs的信号,设计的滤波器的通带阻代也限制在0-fs/2内。
3.程序中的信号幅度值都是1,500点的FFT画出来的幅度值怎么变成了250,应该是1吧?
是的,应该是1。所以怎么变换为1呢,注意到FFT的结果是偶对称的,且其反应的真实频谱是0-fs/2。所以需要的操作是直接取0-N/2的FFT结果,乘以2,然后除以N。即2*abs(y(1:N/2+1))/N,在上面的程序下接着写:
注意到要除以N,也就是FFT的长度,为什么除以N,这个有很多说法,我自己理性理解的也不透彻,所以这里就不解释了自行百度(我只能结合本例子感性理解哈哈)。
%FFT的结果所要展现的真实的频谱幅值[2]
realy=2*abs(y(1:n/2+1))/n;
realf=(0:n/2)*(Fs/n);
figure;
plot(realf,realy);
结果如图4,可见横坐标频率是0-25,纵坐标信号的幅度为1.
然而,有不少人(包括我)平时分析信号都是直接FFT画频谱:图2,且看幅值都是看相对的大小,或者有没有频率分量,就很少做图3与图4的变换。但是我不知道其中的缘由,今天终于明白了。
图4 FFT所反映的真正的频谱幅值