0%

用MATLAB做频谱分析

理论准备

上一篇文章介绍了获取模拟信号频谱的过程:先将模拟信号转换为数字信号,然后对数字信号做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) \end{equation} \]

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

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

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

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

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

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

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

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

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

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

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

非周期信号实例

现有一有限长的非周期模拟信号\(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} \end{equation} \]

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

\[ \begin{equation} X_a(j\omega) = \int_{-\infty}^{+\infty} x_a(t) e^{-j\omega t} dt = \int_{-0.5}^{0.5} (1 + \cos 2\pi t) e^{-j\omega t} dt = -\frac{8\pi^{2} \sin{\frac{\omega}{2}}}{\omega^{3} - 4\pi^{2}\omega} \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来分析连续信号的频谱是有效的。