matlab使用教程(39)—基本频谱分析

        傅里叶变换是用于对时域信号执行频率和功率谱分析的工具。

1.频谱分析数量

        频谱分析研究非均匀采样的离散数据中包含的频谱。傅里叶变换是通过在频率空间表示基于时间或空间的信号来揭示该信号的频率分量的工具。下表列出了用于描述和解释信号属性的常用量。

2含噪信号分析

        傅里叶变换可以计算被随机噪声破坏的信号的频率分量。创建具有 15 Hz 和 40 Hz 分量频率的信号，并插入随机高斯噪声。

rng('default')
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
 + 2.5*randn(size(t)); % Gaussian noise;

        信号的傅里叶变换可确定其频率分量。在 MATLAB® 中，fft 函数使用快速傅里叶变换算法计算傅里叶变换。使用 fft 计算信号的离散傅里叶变换。

y = fft(x);

        将功率谱绘制为频率的函数。尽管噪声在基于时间的空间内伪装成信号的频率分量，但傅里叶变换将其显现为功率尖峰。

n = length(x); % number of samples
f = (0:n-1)*(fs/n); % frequency range
power = abs(y).^2/n; % power of the DFT
plot(f,power)
xlabel('Frequency')
ylabel('Power')

        在许多应用中，查看以 0 频率为中心的功率谱更加方便，因为它能更好地显示信号的周期性。使用fftshift 函数对 y 执行循环平移，并绘制以 0 为中心的功率。

y0 = fftshift(y); % shift y values
f0 = (-n/2:n/2-1)*(fs/n); % 0-centered frequency range
power0 = abs(y0).^2/n; % 0-centered power
plot(f0,power0)
xlabel('Frequency')
ylabel('Power')

3音频信号分析

        您可以使用傅里叶变换来分析音频数据的频谱。

        文件 bluewhale.au 包含水下麦克风记录的加利福尼亚海岸的太平洋蓝鲸发声的音频数据。此文件来自于康奈尔大学生物声学研究项目保存的动物发声库。

        由于蓝鲸的叫声频率如此之低，以至人类几乎听不到。数据中的时间标度压缩了 10 倍，以便提高音调并使叫声更清晰可闻。读取并绘制音频数据。可使用命令 sound(x,fs) 来收听音频。

whaleFile = 'bluewhale.au';
[x,fs] = audioread(whaleFile);
plot(x)
xlabel('Sample Number')
ylabel('Amplitude')

        第一个声音为“颤音”，之后是三个“呻吟音”。本示例将分析单个呻吟音。指定大致包含第一个呻吟音的新数据，并校正时间数据以体现 10 部的加速。将截断的信号绘制为时间的函数。

moan = x(2.45e4:3.10e4);
t = 10*(0:1/fs:(length(moan)-1)/fs);
plot(t,moan)
xlabel('Time (seconds)')
ylabel('Amplitude')
xlim([0 t(end)])

        数据的傅里叶变换确定了音频信号的频率分量。在一些使用 fft 处理大量数据的应用中，通常需要调整输入，使样本数量为 2 的幂。这样可以大幅提高变换计算的速度，对于具有较大质因数的样本大小更是如此。指定新的信号长度 n （2 的幂），并使用 fft 函数计算信号的离散傅里叶变换。 fft 会自动使用零来填充原始数据，以增加样本大小。

m = length(moan); % original sample length
n = pow2(nextpow2(m)); % transform length
y = fft(moan,n); % DFT of signal

        根据加速因子调整频率范围，并计算和绘制信号的功率谱。绘图指示，呻吟音包含约 17 Hz 的基本频率和一系列谐波（其中强调了第二个谐波）。

f = (0:n-1)*(fs/n)/10;
power = abs(y).^2/n; 
plot(f(1:floor(n/2)),power(1:floor(n/2)))
xlabel('Frequency')
ylabel('Power')