1. 从时域到频域信号频谱分析的核心价值信号处理工程师拿到一段时域波形就像拿到一张未经冲洗的底片。我们能看到信号在时间轴上的起伏变化知道它的振幅、频率大概是多少但更精细的“成分”却隐藏在背后。比如一段混杂了50Hz工频干扰、1kHz正弦波和随机噪声的音频信号在时域图上可能就是一团复杂的、周期性不明显的振荡。这时候频谱分析就成为了我们的“显影液”它能将时域信号这个“混合物”分解成不同频率的正弦波“成分”并以一种直观的图形方式——频谱图——展示出来。在Matlab这个强大的计算与可视化平台上实现这一过程不仅高效而且能让我们深入理解信号的频域特性包括幅度谱和相位谱这两个关键维度。幅度谱告诉我们每个频率分量有多“强”而相位谱则揭示了这些分量在时间上的“对齐”关系。无论是通信系统中的调制解调、音频处理中的滤波降噪还是振动分析中的故障诊断频谱分析都是不可或缺的基石。2. 理论基础FFT、幅度谱与相位谱的数学内涵在动手写代码之前我们必须先厘清几个核心概念否则很容易陷入“知其然不知其所以然”的境地导致对结果解读错误。2.1 快速傅里叶变换从连续到离散的桥梁我们理论上学习的傅里叶变换Fourier Transform, FT是针对连续时间、能量有限的信号。但在计算机中我们处理的是经过采样和量化的离散时间信号。离散傅里叶变换Discrete Fourier Transform, DFT就是为离散信号设计的。然而直接计算DFT的复杂度是O(N²)对于长序列来说计算量巨大。快速傅里叶变换FFT不是一种新的变换而是计算DFT的一种高效算法它将复杂度降到了O(N log N)。Matlab中的fft函数就是基于FFT算法实现的。理解这一点至关重要fft(x)得到的结果就是信号x的DFT结果它是一个复数序列包含了我们需要的全部频域信息。2.2 幅度谱信号能量的频率分布图fft输出的结果X fft(x)是一个复数数组。幅度谱Magnitude Spectrum就是这个复数序列的模绝对值。它直观地反映了信号中各个频率分量的强度大小。计算方式为magnitude abs(X)。但这里有几个关键细节单边谱与双边谱由于实数信号的频谱具有共轭对称性其FFT结果的前半部分0到奈奎斯特频率和后半部分是对称的。因此在绘图时我们通常只显示前半部分这就是单边幅度谱。为了保持能量守恒单边谱的幅度除直流分量外需要乘以2。频率轴的映射FFT结果数组的索引k对应的实际频率f需要根据采样频率Fs和点数N来计算f (k-1) * Fs / N对于从1开始索引的Matlab。其中k1对应直流分量0 HzkN/21对应奈奎斯特频率Fs/2。2.3 相位谱频率分量的时间关系图相位谱Phase Spectrum是复数序列X的辐角argument它描述了各频率正弦波分量在时间起点t0处的初始相位。计算方式为phase angle(X)结果单位是弧度。相位信息在信号重构、滤波器设计、通信系统同步等方面至关重要。一个常见的误解是忽略相位谱认为只要幅度谱对就能恢复信号。实际上没有正确的相位信息恢复出的信号波形会面目全非。例如一个方波和它的时移版本幅度谱完全一样但相位谱不同。注意angle函数返回的相位值范围在[-π, π]之间。当相位跨越这个边界时即“相位卷绕”直接绘制的相位谱会出现剧烈的跳变不利于分析。有时需要使用unwrap函数来解卷绕获得平滑的相位变化曲线。3. Matlab实战一步步绘制信号的频谱图理论清晰后我们用一个综合例子来演示全过程。假设我们有一个信号由5Hz、20Hz的两个正弦波和随机噪声叠加而成采样频率为100Hz。3.1 信号生成与参数设置首先我们生成这个测试信号并定义基本参数。清晰的参数设置是后续所有计算的基础。% 参数设置 Fs 100; % 采样频率 (Hz) T 1; % 信号总时长 (秒) N Fs * T; % 采样点数 t (0:N-1)/Fs; % 时间向量 (秒) % 生成信号两个正弦波 噪声 f1 5; % 第一个正弦波频率 (Hz) A1 1.0; % 第一个正弦波幅度 f2 20; % 第二个正弦波频率 (Hz) A2 0.5; % 第二个正弦波幅度 signal A1 * sin(2*pi*f1*t) A2 * sin(2*pi*f2*t pi/4); % f2分量有pi/4的初始相位 noise 0.1 * randn(1, N); % 加入高斯白噪声 x signal noise; % 待分析的合成信号 % 绘制原始时域信号 figure(Position, [100, 100, 800, 400]); subplot(2,1,1); plot(t, x, b-, LineWidth, 1.2); xlabel(时间 (秒)); ylabel(幅度); title(原始时域信号 (含噪声)); grid on; xlim([0, 1]);这段代码生成了我们的测试信号。注意我们在20Hz的分量上特意添加了一个π/445度的初始相位以便在后续的相位谱中观察。加入噪声是为了模拟真实情况让频谱分析更贴近实际应用场景。3.2 执行FFT与计算单边频谱接下来是核心的FFT计算和频谱映射。这里需要仔细处理频率轴和幅度缩放。% 执行FFT X fft(x); % 计算信号的FFT得到复数频谱 % 准备频率向量 (双边谱) f_double (0:N-1) * (Fs/N); % 双边频率轴范围从0到Fs % 计算双边幅度谱和相位谱 magnitude_double abs(X); % 双边幅度谱 phase_double angle(X); % 双边相位谱 (弧度) % 转换为单边谱 (通常更常用) N_half floor(N/2) 1; % 单边谱的点数包含直流和奈奎斯特频率点 f_single f_double(1:N_half); % 单边频率轴 magnitude_single magnitude_double(1:N_half); phase_single phase_double(1:N_half); % 对单边幅度谱进行校正除直流和奈奎斯特频率点外幅度乘2 magnitude_single(2:end-1) 2 * magnitude_single(2:end-1); % 注意如果N是偶数最后一个点奈奎斯特频率点是唯一的不乘2。 % 如果N是奇数则没有单独的奈奎斯特频率点校正范围是2:end。 % 绘制单边幅度谱 subplot(2,1,2); stem(f_single, magnitude_single, r, LineWidth, 1.5, Marker, none); xlabel(频率 (Hz)); ylabel(幅度); title(单边幅度谱); grid on; xlim([0, Fs/2]);这里有几个极易出错的坑点频率向量计算f (0:N-1) * (Fs/N)是正确的。如果错误地写成f (0:N-1) * Fs频率轴会完全错乱。单边谱校正校正因子2只应用于除直流0Hz和奈奎斯特频率Fs/2仅当N为偶数时存在以外的所有频率分量。忘记乘2会导致幅度谱的幅值只有真实值的一半。绘图选择对于频谱这类离散数据使用stem茎状图比plot连线图更合适因为它能清晰表明每个频率点上的谱线。3.3 深入分析相位谱与解卷绕现在我们来观察相位谱并处理可能出现的相位卷绕问题。% 绘制相位谱 figure(Position, [100, 100, 800, 400]); % 绘制原始卷绕相位谱 subplot(2,1,1); stem(f_single, phase_single, b, LineWidth, 1.5, Marker, none); xlabel(频率 (Hz)); ylabel(相位 (弧度)); title(单边相位谱 (卷绕)); grid on; xlim([0, Fs/2]); ylim([-pi, pi]); % 计算并绘制解卷绕后的相位谱 phase_unwrapped unwrap(phase_single); % 解相位卷绕 subplot(2,1,2); plot(f_single, phase_unwrapped, g-, LineWidth, 1.5); xlabel(频率 (Hz)); ylabel(相位 (弧度)); title(单边相位谱 (解卷绕后)); grid on; xlim([0, Fs/2]);unwrap函数通过检测相邻相位差是否超过π默认阈值来消除±2π的跳变从而获得连续的相位变化。这在分析系统相频特性或需要精确相位值时非常有用。在我们的例子中由于噪声的存在非信号频率点的相位是随机的但在5Hz和20Hz附近你应该能看到相对稳定的相位值20Hz处接近我们设定的π/4。3.4 功率谱密度另一种视角除了幅度谱功率谱密度Power Spectral Density, PSD在工程中同样重要它反映了信号功率在频域的分布。Matlab提供了pwelch函数来估计PSD这是一种基于Welch方法的更稳健的估计方式尤其适用于随机信号或噪声。% 使用pwelch方法估计功率谱密度 figure; [pxx, f_psd] pwelch(x, hamming(256), 128, 256, Fs); % 使用汉明窗256点分段50%重叠 plot(f_psd, 10*log10(pxx), m-, LineWidth, 1.5); % 转换为dB单位 xlabel(频率 (Hz)); ylabel(功率/频率 (dB/Hz)); title(Welch方法估计的功率谱密度 (PSD)); grid on; xlim([0, Fs/2]);pwelch函数通过将数据分段、加窗、计算周期图再平均来减少直接使用abs(fft).^2带来的方差得到更平滑、统计特性更好的谱估计。参数选择窗函数、窗长、重叠率会影响频率分辨率和估计方差需要根据实际信号权衡。4. 关键技巧与深度避坑指南掌握了基本流程后一些进阶技巧和常见陷阱决定了你分析的可靠性与专业性。4.1 窗函数的选择与应用艺术我们之前计算FFT时默认对整段数据进行了矩形窗截断。这相当于在时域乘以一个矩形窗在频域会导致频谱泄漏Spectral Leakage——即一个频率的能量“泄漏”到其他频率点上表现为频谱图中的“拖尾”现象会模糊靠近的频谱峰并抬高噪声基底。为了抑制泄漏需要在FFT前对信号加窗。常用窗函数有汉明窗最通用的选择能有效抑制旁瓣主瓣宽度适中。汉宁窗旁瓣衰减比汉明窗更好但主瓣稍宽。布莱克曼窗旁瓣抑制最好主瓣最宽频率分辨率最低。凯泽窗可通过参数β灵活调节主瓣宽度与旁瓣衰减的权衡。% 加窗处理对比 window_hamming hamming(N); % 生成汉明窗转置为行向量 x_windowed x .* window_hamming; % 加窗 X_nowin fft(x); X_win fft(x_windowed); % 计算加窗带来的幅度补偿因子窗函数的相干增益 coherent_gain sum(window_hamming)/N; magnitude_compensated abs(X_win) / coherent_gain; % 对比加窗前后的幅度谱单边已校正 figure; hold on; stem(f_single, magnitude_single, b, DisplayName, 无窗); stem(f_single, magnitude_compensated(1:N_half), r, DisplayName, 汉明窗(补偿后)); xlabel(频率 (Hz)); ylabel(幅度); title(加窗对幅度谱的影响对比); legend; grid on; xlim([0, 50]);重要提示加窗会损失信号能量因此计算出的幅度谱幅值会变小。为了得到正确的幅度估计必须对结果进行幅度补偿除以窗函数的平均高度相干增益。sum(window)/length(window)。pwelch等函数内部已经自动处理了这些补偿。4.2 频谱分辨率、补零与频率插值频谱分辨率Δf Fs / N它决定了你能区分开两个频率的最小间隔。增加信号时长T从而增加N是提高分辨率的唯一根本方法。但有时数据长度固定我们想通过“补零”来获得更平滑的频谱图。% 补零操作 N_zeropad 4 * N; % 补零到原长度的4倍 x_zeropad [x, zeros(1, N_zeropad - N)]; % 在信号后补零 X_zeropad fft(x_zeropad, N_zeropad); f_zeropad (0:N_zeropad-1) * (Fs/N_zeropad); magnitude_zeropad abs(X_zeropad); % 绘制对比 figure; subplot(2,1,1); stem(f_single, magnitude_single, b., MarkerSize, 10); hold on; plot(f_zeropad(1:2*N_half), magnitude_zeropad(1:2*N_half), r-, LineWidth, 0.5); xlabel(频率 (Hz)); ylabel(幅度); title(补零效果对比 (红线为补零后插值曲线)); legend(原始FFT点, 补零后插值); grid on; xlim([0, 30]); subplot(2,1,2); % 局部放大观察 idx find(f_single 4 f_single 6); stem(f_single(idx), magnitude_single(idx), b., MarkerSize, 15); hold on; idx_z find(f_zeropad 4 f_zeropad 6); plot(f_zeropad(idx_z), magnitude_zeropad(idx_z), r-, LineWidth, 1); xlabel(频率 (Hz)); ylabel(幅度); title(5Hz峰值的局部放大图); legend(原始点, 补零插值); grid on;必须明确补零不能提高频谱分辨率它只是对已有的频谱进行了插值让曲线看起来更平滑有助于更精确地通过观察找到谱峰的大致位置例如通过findpeaks函数。分辨率依然由原始数据长度决定。上图中红线的峰值点更多但峰值的实际宽度由主瓣决定并没有变窄。4.3 使用 findpeaks 函数精确提取谱峰对于自动化的频谱分析我们需要从幅度谱中精确提取峰值频率和幅度。Matlab的findpeaks函数非常强大。% 使用findpeaks在单边幅度谱中寻找峰值 [magnitude_peaks, peak_locs] findpeaks(magnitude_single, f_single, ... MinPeakHeight, 0.2, ... % 最小峰高用于抑制噪声 MinPeakDistance, 2); % 最小峰间距(Hz)避免在同一个峰附近检测到多个点 % 绘制频谱并标记峰值 figure; stem(f_single, magnitude_single, b, LineWidth, 1, Marker, none); hold on; plot(peak_locs, magnitude_peaks, rv, MarkerSize, 10, MarkerFaceColor, r); text(peak_locs, magnitude_peaks, ... cellstr(num2str([peak_locs, magnitude_peaks], (%.1fHz, %.2f))), ... VerticalAlignment, bottom, HorizontalAlignment, center, FontSize, 9); xlabel(频率 (Hz)); ylabel(幅度); title(使用findpeaks自动检测频谱峰值); grid on; xlim([0, Fs/2]); disp(检测到的峰值频率和幅度); disp([peak_locs, magnitude_peaks]);findpeak的参数调优是关键MinPeakHeight根据信号的噪声水平设置可以有效滤除噪声引起的假峰。MinPeakDistance根据你预期的信号最小频率间隔设置可以防止在一个较宽的谱峰上检测到多个点。其他参数如MinPeakProminence最小峰突出度在复杂频谱中也非常有用它能更好地识别“真正”的峰。5. 综合案例分析一段实际音频信号让我们将所学应用于一个更接近实际的场景分析一段包含特定频率蜂鸣声的音频文件。% 假设我们有一个音频文件 beep_signal.wav 其中包含1kHz的蜂鸣声 % [audio_data, Fs_audio] audioread(beep_signal.wav); % x_audio audio_data(:,1); % 取单声道 % 为演示我们模拟生成一段类似的音频信号 Fs_audio 8000; t_audio 0:1/Fs_audio:0.5; % 0.5秒时长 f_beep 1000; % 1kHz蜂鸣声 x_audio 0.8 * sin(2*pi*f_beep*t_audio); % 加入一些背景噪声和谐波 x_audio x_audio 0.1*sin(2*pi*2*f_beep*t_audio) 0.05*randn(size(t_audio)); N_audio length(x_audio); % 1. 绘制时域波形 figure(Position, [50, 50, 1200, 600]); subplot(2,2,1); plot(t_audio, x_audio, b-); xlabel(时间 (秒)); ylabel(幅度); title(音频信号时域波形); grid on; xlim([0, 0.05]); % 只看前50ms % 2. 计算并绘制幅度谱 (加汉宁窗) window_audio hanning(N_audio); x_win_audio x_audio .* window_audio; X_audio fft(x_win_audio); f_audio (0:N_audio-1) * (Fs_audio/N_audio); magnitude_audio abs(X_audio); % 单边谱处理 N_half_audio floor(N_audio/2)1; f_audio_single f_audio(1:N_half_audio); magnitude_audio_single magnitude_audio(1:N_half_audio); magnitude_audio_single(2:end-1) 2 * magnitude_audio_single(2:end-1) / (sum(window_audio)/N_audio); % 幅度补偿 subplot(2,2,2); plot(f_audio_single, magnitude_audio_single, r-, LineWidth, 1.5); xlabel(频率 (Hz)); ylabel(幅度); title(音频信号单边幅度谱 (加汉宁窗)); grid on; xlim([0, 4000]); % 显示到4kHz % 标记1kHz主峰和2kHz谐波 hold on; plot([1000, 2000], interp1(f_audio_single, magnitude_audio_single, [1000, 2000]), ko, MarkerSize, 8, MarkerFaceColor, y); % 3. 使用pwelch估计PSD subplot(2,2,3); [pxx_audio, f_pwelch] pwelch(x_audio, hamming(512), 256, 512, Fs_audio); plot(f_pwelch, 10*log10(pxx_audio), m-, LineWidth, 1.5); xlabel(频率 (Hz)); ylabel(功率谱密度 (dB/Hz)); title(Welch方法估计的PSD); grid on; xlim([0, 4000]); % 4. 绘制相位谱 (解卷绕) phase_audio angle(X_audio); phase_audio_single phase_audio(1:N_half_audio); phase_audio_unwrapped unwrap(phase_audio_single); subplot(2,2,4); plot(f_audio_single, phase_audio_unwrapped, g-, LineWidth, 1.5); xlabel(频率 (Hz)); ylabel(相位 (弧度)); title(音频信号相位谱 (解卷绕)); grid on; xlim([0, 4000]);这个综合案例展示了从时域观察、加窗幅度谱分析、功率谱估计到相位谱查看的完整流程。在实际工作中你可能需要根据信号特性如是否平稳和分析目的如测量频率、观察谐波、分析噪声来选择不同的谱估计方法和参数。例如对于瞬态信号如一个脉冲可能更适合直接观察其FFT幅度谱而对于平稳随机噪声pwelch估计的PSD更能反映其统计特性。
)
