0%

用MATLAB做频谱分析

模拟信号

上一篇文章介绍了获取模拟信号频谱的过程:先将模拟信号转换为数字信号,然后对数字信号做DFT,得到的频谱序列能够反映模拟信号频谱的特征。接下来在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
11
12
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);
n = length(x); % number of samples
f = (0:n-1)/n*fs; % frequency range
power = abs(y).^2/n; % power of the DFT

figure(1); plot(f(1:floor(n/2)), power(1:floor(n/2)));
xlabel('Frequency'); ylabel('Power');

第1到4行代码是采样过程,第6行使用fft计算频谱,第8行的代码是做频率轴尺度换算的,从DFT序列的下标\(k\)换算到实际的频率\(F\),过程有点绕,要详细说明一下。从上一篇文章中我们知道,\(X(k)\)能反映真实的频谱\(X(j\Omega)\)的特性,以\(X(e^{j\omega})\)为桥梁。整个过程记为

\[ X(k) \to X(e^{j\omega}) \to X(j\Omega) \to X(j2\pi F) \]

对应的频率变量换算过程为

\[ k \to \omega \to \Omega \to F \]

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

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

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

\[ \Omega = \frac{\omega}{T_s} = \omega F_s \]

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

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

\(\Omega\)是角频率,一般用频率\(F\)更直观,于是

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

这就是第8行代码的来历。根据傅里叶变换的特性,只显示前一半频率轴的图像就可以了,对应第11到第12行代码。得到的频谱图为

信号的频谱

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

音频信号

分析音频时,使用audioread函数读取音频文件的数据,同时返回采样频率。如下的示例代码分析MATLAB自带音频文件的频谱,采样频率为8000Hz:

1
2
3
4
5
6
7
8
9
10
audioFile = 'dft_voice_8kHz.wav';
[x, fs] = audioread(audioFile);

n = length(x);
y = fft(x);
f = (0:n-1)*/n*fs;
power = abs(y).^2/n;

subplot(2, 1, 1); plot(x);
subplot(2, 1, 2); plot(f(1:floor(n/2)), power(1:floor(n/2)));

音频信号及其频谱如下图所示:

音频信号及其频谱

从频谱可以看出,频率成分大多数分布在500Hz以下。