用MATLAB做频谱分析

发表于 分类于 Waline： 本文字数： 1.4k 阅读时长 ≈ 5 分钟

理论准备

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

具体来说，\(X(k)\)能反映真实的频谱\(X(j\Omega)\)的特性，以\(X(e^{j\omega})\)为桥梁。\(X(j\Omega)\)中的自变量\(\Omega\)是角频率，单位是rad/s，还要转换为更习惯的以Hz为单位的频率变量\(F\)。整个转换过程记为

\[ \begin{equation} X(k) \to X(e^{j\omega}) \to X(j\Omega) \to X(j2\pi F) \tag{1} \end{equation} \]

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

\[ \begin{equation} k \to \omega \to \Omega \to F \tag{2} \end{equation} \]

\(k\)的范围是\(0\)到\(N-1\)，记为半开半闭区间\([0, N)\)，对应的\(\omega\)的范围是\([0, 2\pi)\)，它们之间的转换公式为

\[ \begin{equation} \omega = \frac{2\pi}{N}k \tag{3} \end{equation} \]

\(\Omega\)与\(\omega\)之间的关系为

\[ \begin{equation} \Omega = \frac{\omega}{T_s} = \omega F_s \tag{4} \end{equation} \]

所以，\(\Omega\)与\(k\)的转换关系为

\[ \begin{equation} \Omega = \frac{2\pi k}{N} F_s \tag{5} \end{equation} \]

角频率\(\Omega\)与频率\(F\)的转换关系为

\[ \begin{equation} F = \frac{\Omega}{2\pi} = \frac{k}{N} F_s \tag{6} \end{equation} \]

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

周期信号

先看周期信号的例子，本例来源于MATLAB官方。假设现有一模拟信号\(x(t)\)，它由两个不同频率的正弦波叠加得到的，频率分别为\(f_1=15 \mathrm{Hz}\)和\(f_2=40 \mathrm{Hz}\)，解析表达式

\[ x(t) = 1.3\sin \omega_1 t + 1.7\sin \omega_2 (t-2) \]

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

1
2
3
4
5
6
7
8
9
10
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行根据公式\((6)\)做频率轴尺度换算。DFT的幅值关于中点\(N/2\)对称，所以可以只显示前一半频率轴的图像，对应第11到第12行代码。得到的频谱图为

图1 连续周期信号的频谱

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

非周期信号

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

\[ \begin{equation} x_a(t) = \begin{cases} 1 + \cos 2\pi t, & -0.5 \le t \le 0.5 \\ 0, & 其他 \end{cases} \tag{7} \end{equation} \]

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

\[ \begin{equation} X_a(j\omega) = \int_{-\infty}^{+\infty} x_a(t) e^{-j\omega t} \mathrm{d}t = \int_{-0.5}^{0.5} (1 + \cos 2\pi t) e^{-j\omega t} \mathrm{d}t = -\frac{8\pi^{2} \sin{\frac{\omega}{2}}}{\omega^{3} - 4\pi^{2}\omega} \tag{8} \end{equation} \]

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

图2 连续非周期信号及其频谱

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
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，将原信号周期延拓之后，再从下标\(0\)到\(N-1\)取一个周期，即为第11行代码得到的序列。第13行代码使用fftshift函数将数据移位，得到关于频率0对称的数据。第14行代码再构造一个关于频率0对称的频率轴，频率转换关系同样遵循公式\((6)\)，接下来就可以画图显示。代码的运行结果如图3所示，分别为采样序列\(x_a(n)\)与它的DFT的图像。

图3 采样序列及其DFT

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