1. 声学信号分析基础入门第一次接触声学信号分析时我被各种专业术语搞得晕头转向。直到把麦克风录到的声音信号导入MATLAB看到时域波形变成频谱图的那一刻才真正理解频域分析的魅力。简单来说时域信号告诉我们声音随时间的变化而频域分析则揭示了声音由哪些频率组成 - 就像把一道复合光分解成不同颜色的光谱。频谱分析的核心是傅里叶变换这个数学工具能将时域信号转换为频域表示。我常跟学生打比方时域信号像是一杯混合果汁而傅里叶变换就是榨汁机的过滤网能把橙子、苹果、西瓜的成分分别分离出来。在声学分析中我们最常用的是快速傅里叶变换FFT它就像是个超级高效的频率分解器。实际工程中会遇到几个关键参数采样频率决定了能分析的最高频率Nyquist频率FFT点数影响频率分辨率点数越多分辨率越高窗函数选择常用汉宁窗来减少频谱泄漏平均次数通过多次平均提高信噪比% 基本FFT分析示例 fs 44100; % 采样率 t 0:1/fs:1; % 1秒时间向量 x 0.7*sin(2*pi*1000*t) 0.3*cos(2*pi*5000*t); % 合成信号 nfft 2^14; % FFT点数 window hann(length(x)); % 汉宁窗 [Pxx,f] periodogram(x,window,nfft,fs); % 功率谱估计 plot(f,10*log10(Pxx)); % 绘制频谱图2. 从原始信号到频谱的完整处理流程拿到一段录音或振动信号后我通常会按照这个标准化流程处理2.1 数据预处理原始信号往往包含直流偏移和高频噪声。有次分析工厂噪声时就因为没有做去趋势处理导致低频部分出现异常峰值。现在我的预处理清单必含这三步去趋势消除信号基线漂移带通滤波保留感兴趣频段如20Hz-20kHz数据分段对长信号分帧处理x detrend(x); % 去趋势 [b,a] butter(4,[20 20000]/(fs/2),bandpass); % 4阶带通滤波器 x_filtered filtfilt(b,a,x); % 零相位滤波2.2 频谱计算技巧直接做FFT常会遇到频谱泄漏问题。经过多次测试我总结出几个实用技巧窗函数选择汉宁窗适合大多数情况平顶窗适合幅值测量重叠分段75%重叠率能提高数据利用率补零处理适当补零可以平滑频谱曲线nfft 2^16; % 补零到65536点 window hann(4096); % 汉宁窗 noverlap 3072; % 75%重叠 [Pxx,f] pwelch(x,window,noverlap,nfft,fs); % Welch法功率谱估计 semilogx(f,10*log10(Pxx)); % 对数频率坐标2.3 声压级计算声学分析中常用声压级(SPL)表示信号强度计算公式为SPL 20*log10(p/p0)其中p020μPa是基准声压。在MATLAB中实现时要注意单位转换p0 20e-6; % 基准声压 p_rms sqrt(mean(x.^2)); % 有效值 SPL 20*log10(p_rms/p0); % 总声压级3. 1/3倍频程分析的工程实践3.1 倍频程概念解析第一次接触1/3倍频程时我困惑为什么要用这种非均匀频带。后来在噪声评估项目中才明白人耳对频率的感知是对数关系而1/3倍频程正好模拟了这种特性。标准中心频率按公式计算fc 1000 * 10^(3n/30)其中n为频带编号从-16到13对应20Hz到20kHz。3.2 MATLAB实现细节写1/3倍频程分析代码时我踩过三个坑频带边界计算不精确导致能量泄漏未考虑FFT的频率分辨率对数坐标显示处理不当经过优化后的核心代码如下% 标准1/3倍频程中心频率 fc [20 25 31.5 40 50 63 80 100 125 160 200 ... 250 315 400 500 630 800 1000 1250 1600 2000 ... 2500 3150 4000 5000 6300 8000 10000 12500 16000]; % 计算上下限频率 fl fc/(2^(1/6)); fu fc*(2^(1/6)); % 匹配FFT频率点 f_resolution fs/nfft; % 频率分辨率 bin_low round(fl/f_resolution)1; bin_high round(fu/f_resolution)1; % 计算各频带能量 for i 1:length(fc) band_energy(i) sum(Pxx(bin_low(i):bin_high(i))); end % 转换为声压级 SPL_1_3 10*log10(band_energy/(p0^2));3.3 结果可视化技巧好的可视化能事半功倍我常用的绘图设置包括对数频率坐标参考等高线专业配色方案figure(Color,w); semilogx(fc,SPL_1_3,-o,LineWidth,1.5); xlabel(中心频率 (Hz)); ylabel(声压级 (dB)); title(1/3倍频程分析); set(gca,XScale,log,XTick,fc(1:3:end)); grid on;4. 典型工程案例解析去年参与的一个工业风机噪声分析项目完整展示了这套方法的实用价值。客户提供的原始噪声信号时域波形杂乱无章但经过频谱和1/3倍频程分析后我们清晰地识别出以下几个特征4.1 异常频率定位通过频谱分析发现基频125Hz对应电机转速250Hz处谐波异常升高800Hz宽带噪声突显% 异常频率检测 [f_peaks, locs] findpeaks(Pxx,MinPeakHeight,max(Pxx)/10); significant_peaks f(locs(f_peaks mean(Pxx)3*std(Pxx)));4.2 1/3倍频程诊断1/3倍频程分析更直观地显示125Hz频带超标15dB高频段整体抬升特定频带出现凹陷这引导我们检查风机的叶片结构和轴承状况最终发现叶片动平衡失调和轴承磨损问题。4.3 报告生成自动化为提升效率我开发了自动生成分析报告的脚本% 创建Word报告 import mlreportgen.dom.*; doc Document(噪声分析报告,docx); append(doc,Heading1(频谱分析结果)); append(doc,Image(which(spectrum_plot.png))); append(doc,Heading2(1/3倍频程评估)); table Table({频带(Hz),声压级(dB); fc, SPL_1_3}); append(doc,table); close(doc);这套方法后来被客户采纳为标准分析流程我也因此养成了保存完整分析脚本的习惯。每次遇到新项目只需调整输入参数就能快速得到专业级的分析结果。
