用MATLAB做频谱分析

发表于 2021-03-10 分类于大学课程 Waline：本文字数： 1.4k 阅读时长 ≈ 5 分钟

理论准备

上一篇文章介绍了获取模拟信号频谱的过程：先将模拟信号转换为数字信号，然后对数字信号做DFT，得到的频谱序列能够反映模拟信号频谱的特征。

具体来说，能反映真实的频谱的特性，以为桥梁。中的自变量是角频率，单位是rad/s，还要转换为更习惯的以Hz为单位的频率变量。整个转换过程记为

得到的值之后，频率变量换算的过程非常关键，这样才知道每个频率成分所占的比重。换算的过程为

的范围是到，记为半开半闭区间，对应的的范围是，它们之间的转换公式为

与之间的关系为

所以，与的转换关系为

角频率与频率的转换关系为

公式就是最终要在代码中实现的频率轴换算公式。频谱的大小则要将DFT的值乘以采样周期。

周期信号

先看周期信号的例子，本例来源于MATLAB官方。假设现有一模拟信号，它由两个不同频率的正弦波叠加得到的，频率分别为和，解析表达式

那么，角频率、，在MATLAB中对该模拟信号进行采样，然后通过fft函数计算其离散傅里叶变换DFT的代码如下：

fs = 100;                                % sample frequency (Hz)
t = 0:1/fs:10-1/fs;                      % 10 second span time vector
x = (1.3)*sin(2*pi*15*t) ...             % 15 Hz component
  + (1.7)*sin(2*pi*40*(t-2));            % 40 Hz component
y = fft(x)*1/fs;        % multiply sample period
n = length(x);          % number of samples
f = (0:n-1)/n*fs;       % frequency range
y = abs(y);             % magnitude of spectrum
figure(1); plot(f(1:floor(n/2)), y(1:floor(n/2)));
xlabel('Frequency (Hz)'); ylabel('Magnitude');

第1到4行代码是采样过程，第6行使用fft计算频谱，第8行根据公式做频率轴尺度换算。DFT的幅值关于中点对称，所以可以只显示前一半频率轴的图像，对应第11到第12行代码。得到的频谱图为

只有15Hz和40Hz两个分量，跟的频率成分是吻合的。本例中的采样频率为100Hz，大于两倍的40Hz，所以能将两个频率分量都计算出来。读者可以尝试减小采样频率，看看会得到什么样的频谱图。

非周期信号

现有一有限长的非周期模拟信号，它的解析表达式为

$其他$

的频谱用傅里叶变换描述，借助软件工具，可以计算得到

该频谱是个实函数。信号及其频谱的图像如图2所示，频谱图横轴的单位为Hz，可以发现大于8Hz的频率分量趋于0，可以认为它是个带宽有限的信号，带宽小于8Hz。

现在对这个持续时间为2s的信号做频谱分析，根据采样定理，采样频率取16Hz即可。全部代码如下：

T = 2;       % 信号的持续时间为2s，也是做周期延拓时的周期
fs = 16;     % 采样频率
N = T*fs;    % 采样点个数

% 采样过程，得到离散信号
t = -T/2:T/N:T/2-T/N;
xn = 1 + cos(2*pi*t);
k1 = find(t < -0.5); xn(k1) = 0;
k2 = find(t > 0.5); xn(k2) = 0;

x = [xn(N/2+1:N), xn(1:N/2)];  % 将负半轴的数据接到正半轴数据的后面
y = fft(x)*1/fs;               % 将DFT的值乘以采样周期
yshift = fftshift(real(y));    % 移位，使得以频率0为中心
fshift = (-N/2:N/2-1)/N*fs;    % 构造以频率0为中心的频率轴

figure(2)
subplot(2, 1, 1); stem(0:N-1, xn, '.'); xlim([0, N]); grid on
subplot(2, 1, 2); stem(fshift, yshift, '.'); grid on

第5到9行代码是根据信号的解析表达式得到离散信号的过程。第11行代码将负半轴的数据接到正半轴数据的后面，然后交给第12行计算DFT。因为DFT的本质是信号的周期延拓的DFS，将原信号周期延拓之后，再从下标到取一个周期，即为第11行代码得到的序列。第13行代码使用fftshift函数将数据移位，得到关于频率0对称的数据。第14行代码再构造一个关于频率0对称的频率轴，频率转换关系同样遵循公式，接下来就可以画图显示。代码的运行结果如图3所示，分别为采样序列与它的DFT的图像。

很明显，图2和图3的结果是一致的。这就证明了用DFT来分析连续信号的频谱是有效的。