企业官网建设流程全解析

FFT、STFT与DWT时频分析实战：Python代码实现与信号重构误差对比

引言：非平稳信号处理的挑战与机遇

当我们面对现实世界中的振动监测、语音识别或生物医学信号分析时，传统傅里叶变换的局限性变得尤为明显。想象一下工厂里一台运转中的机械设备——它的振动信号可能包含轴承磨损初期产生的高频瞬态冲击，同时伴随着电机运转的低频周期性振动。这种同时包含瞬态和稳态成分的信号，正是典型的非平稳信号。

本文将带您用Python实现三种核心时频分析方法：快速傅里叶变换(FFT)、短时傅里叶变换(STFT)和离散小波变换(DWT)。我们不仅会对比它们的时频表示能力，还将量化评估信号重构误差，为工程实践提供直观的性能参考。所有代码示例均基于SciPy和PyWavelets库，可直接应用于您的实际项目。

1. 实验环境配置与测试信号生成

1.1 基础环境准备

首先确保安装必要的Python库：

pip install numpy scipy matplotlib pywavelets

我们创建一个包含多种频率成分的测试信号，模拟实际工程中的复杂振动：

import numpy as np import matplotlib.pyplot as plt # 信号参数 sample_rate = 1000 # 采样率(Hz) duration = 2.0 # 信号时长(s) t = np.linspace(0, duration, int(sample_rate*duration), endpoint=False) # 生成多组分测试信号 def generate_test_signal(t): signal = np.zeros_like(t) # 0-0.5s: 10Hz正弦波 mask = (t >= 0) & (t < 0.5) signal[mask] = np.sin(2*np.pi*10*t[mask]) # 0.5-1.0s: 50Hz正弦波 + 高斯脉冲 mask = (t >= 0.5) & (t < 1.0) signal[mask] = 0.8*np.sin(2*np.pi*50*t[mask]) signal[mask] += 0.5*np.exp(-((t[mask]-0.7)*100)**2)*np.sin(2*np.pi*120*t[mask]) # 1.0-1.5s: 线性扫频信号 mask = (t >= 1.0) & (t < 1.5) freq = 20 + 60*(t[mask]-1.0)/0.5 # 20Hz到80Hz线性变化 signal[mask] = 0.6*np.sin(2*np.pi*freq*t[mask]) # 1.5-2.0s: 方波信号(30Hz) mask = (t >= 1.5) & (t < 2.0) signal[mask] = 0.4*np.sign(np.sin(2*np.pi*30*t[mask])) return signal signal = generate_test_signal(t)

1.2 信号可视化分析

使用Matplotlib绘制时域波形和频谱图：

plt.figure(figsize=(12, 8)) # 时域波形 plt.subplot(2, 1, 1) plt.plot(t, signal) plt.title('Test Signal in Time Domain') plt.xlabel('Time (s)') plt.ylabel('Amplitude') plt.grid(True) # 频谱分析 plt.subplot(2, 1, 2) freqs = np.fft.rfftfreq(len(signal), 1/sample_rate) fft = np.abs(np.fft.rfft(signal)) plt.plot(freqs, fft) plt.title('Frequency Spectrum') plt.xlabel('Frequency (Hz)') plt.ylabel('Magnitude') plt.grid(True) plt.tight_layout() plt.show()

这个测试信号包含以下特征：

• 稳态正弦波成分(10Hz和50Hz)
• 瞬态高斯调制脉冲(中心频率120Hz)
• 线性变化的扫频信号(20-80Hz)
• 非线性方波信号(30Hz)

2. FFT分析与重构实现

2.1 FFT基本原理与实现

FFT是离散傅里叶变换(DFT)的高效算法，适合分析平稳信号：

def fft_analysis(signal, sample_rate): n = len(signal) freqs = np.fft.rfftfreq(n, 1/sample_rate) fft_coeff = np.fft.rfft(signal) return freqs, fft_coeff def fft_reconstruct(fft_coeff, n): return np.fft.irfft(fft_coeff, n) # 执行FFT分析 freqs, fft_coeff = fft_analysis(signal, sample_rate) reconstructed = fft_reconstruct(fft_coeff, len(signal))

2.2 FFT重构误差评估

计算重构信号与原信号的误差指标：

def calculate_errors(original, reconstructed): mse = np.mean((original - reconstructed)**2) max_error = np.max(np.abs(original - reconstructed)) snr = 10*np.log10(np.var(original)/np.var(original - reconstructed)) return {'MSE': mse, 'Max Error': max_error, 'SNR (dB)': snr} fft_errors = calculate_errors(signal, reconstructed) print("FFT Reconstruction Errors:", fft_errors)

FFT的局限性在时频分析中表现明显：

• 无法定位频率成分的时间信息
• 对瞬态信号分析效果差
• 重构方波信号时出现Gibbs现象

3. STFT分析与参数优化

3.1 STFT核心参数选择

STFT通过加窗分段分析实现时频局部化：

from scipy.signal import stft, istft def stft_analysis(signal, sample_rate, window='hann', nperseg=256, noverlap=None): f, t, Zxx = stft(signal, sample_rate, window=window, nperseg=nperseg, noverlap=noverlap) return f, t, Zxx def stft_reconstruct(Zxx, sample_rate, window='hann', nperseg=256, noverlap=None): _, xrec = istft(Zxx, sample_rate, window=window, nperseg=nperseg, noverlap=noverlap) return xrec[:len(signal)] # 截取与原始信号等长部分

3.2 窗口大小对分辨率的影响

比较不同窗口尺寸下的时频表现：

windows = [64, 128, 256, 512] plt.figure(figsize=(15, 10)) for i, nperseg in enumerate(windows, 1): f, t, Zxx = stft_analysis(signal, sample_rate, nperseg=nperseg) plt.subplot(2, 2, i) plt.pcolormesh(t, f, np.abs(Zxx), shading='gouraud') plt.title(f'STFT (Window Size={nperseg})') plt.ylabel('Frequency [Hz]') plt.xlabel('Time [sec]') plt.colorbar() plt.tight_layout() plt.show()

窗口选择的关键权衡：

• 窄窗口：时间分辨率高，频率分辨率低
• 宽窗口：频率分辨率高，时间分辨率低
• 汉宁窗：平衡主瓣宽度和旁瓣衰减

3.3 STFT重构误差对比

选择最优窗口进行重构误差分析：

# 使用128点汉宁窗 f, t, Zxx = stft_analysis(signal, sample_rate, nperseg=128) stft_recon = stft_reconstruct(Zxx, sample_rate, nperseg=128) stft_errors = calculate_errors(signal, stft_recon) print("STFT Reconstruction Errors:", stft_errors)

STFT在时频分析上的改进：

• 能够定位频率随时间变化
• 对瞬态信号有更好表现
• 但仍受限于固定分辨率

4. 小波变换的多分辨率分析

4.1 DWT基本原理与实现

PyWavelets库提供了完善的DWT实现：

import pywt def dwt_analysis(signal, wavelet='db4', level=5): coeffs = pywt.wavedec(signal, wavelet, level=level) return coeffs def dwt_reconstruct(coeffs, wavelet='db4'): return pywt.waverec(coeffs, wavelet) # 执行5层小波分解 coeffs = dwt_analysis(signal) dwt_recon = dwt_reconstruct(coeffs) dwt_errors = calculate_errors(signal, dwt_recon) print("DWT Reconstruction Errors:", dwt_errors)

4.2 小波基函数选择比较

不同小波基的特性对比：

wavelets = ['db4', 'sym4', 'coif2', 'haar'] errors = {} for wav in wavelets: coeffs = dwt_analysis(signal, wavelet=wav) recon = dwt_reconstruct(coeffs, wavelet=wav) errors[wav] = calculate_errors(signal, recon) # 显示不同小波的误差对比 for wav, err in errors.items(): print(f"{wav} - MSE: {err['MSE']:.2e}, SNR: {err['SNR (dB)']:.1f} dB")

4.3 小波时频图绘制

小波尺度图可视化：

def plot_scalogram(signal, scales, wavelet='morl'): coefficients, frequencies = pywt.cwt(signal, scales, wavelet) plt.figure(figsize=(10, 6)) plt.imshow(np.abs(coefficients), extent=[0, duration, 1, max(scales)], cmap='jet', aspect='auto', vmax=abs(coefficients).max(), vmin=-abs(coefficients).max()) plt.colorbar() plt.title('Wavelet Scalogram') plt.ylabel('Scale') plt.xlabel('Time (s)') plt.show() # 使用墨西哥帽小波绘制尺度图 scales = np.arange(1, 128) plot_scalogram(signal, scales)

小波变换的优势体现：

• 高频成分时间分辨率高
• 低频成分频率分辨率高
• 自适应时频分辨率
• 对瞬态信号捕捉能力强

5. 三种方法综合对比

5.1 重构误差定量比较

汇总三种方法的误差指标：

5.2 时频定位能力对比

三种方法在测试信号各段的表现：

1. 稳态正弦段(0-0.5s)

   • FFT: 频率定位精确，无时间信息
   • STFT: 时频定位良好
   • DWT: 低频分辨率优秀

2. 瞬态脉冲段(0.5-1.0s)

   • FFT: 完全无法定位
   • STFT: 可检测但分辨率有限
   • DWT: 精确捕捉瞬态时刻

3. 扫频段(1.0-1.5s)

   • FFT: 仅显示频带范围
   • STFT: 可跟踪频率变化
   • DWT: 连续尺度变化明显

4. 方波段(1.5-2.0s)

   • FFT: Gibbs现象严重
   • STFT: 谐波随时间分布
   • DWT: 多尺度表征非线性

5.3 工程应用选型建议

根据实际需求选择合适方法：

推荐FFT当：

• 信号严格平稳
• 仅需频率成分分析
• 计算资源有限

推荐STFT当：

• 需要平衡时频分辨率
• 信号变化相对缓慢
• 实现简单性优先

推荐DWT当：

• 信号包含瞬态突变
• 需要多尺度分析
• 计算资源允许

6. 高级应用与实战技巧

6.1 信号去噪实战

利用小波阈值去噪：

def wavelet_denoise(signal, wavelet='db4', level=5, mode='soft'): # 小波分解 coeffs = pywt.wavedec(signal, wavelet, level=level) # 计算阈值(通用阈值法) sigma = np.median(np.abs(coeffs[-level])) / 0.6745 threshold = sigma * np.sqrt(2*np.log(len(signal))) # 应用阈值 new_coeffs = [] new_coeffs.append(coeffs[0]) # 保留近似系数 for i in range(1, len(coeffs)): new_coeffs.append(pywt.threshold(coeffs[i], threshold, mode=mode)) # 重构信号 return pywt.waverec(new_coeffs, wavelet) # 添加高斯噪声 noisy_signal = signal + 0.1*np.random.randn(len(signal)) denoised = wavelet_denoise(noisy_signal) # 计算去噪效果 noise_reduction = 10*np.log10(np.var(signal)/np.var(denoised - signal)) print(f"Noise Reduction: {noise_reduction:.1f} dB")

6.2 实时处理架构建议

对于实时信号处理系统：

from collections import deque class RealTimeProcessor: def __init__(self, window_size=256, wavelet='db4'): self.buffer = deque(maxlen=window_size*2) self.window_size = window_size self.wavelet = wavelet def process_chunk(self, new_data): self.buffer.extend(new_data) if len(self.buffer) >= self.window_size: # 取出最新数据段 segment = list(self.buffer)[-self.window_size:] # 并行执行三种分析 fft_result = np.fft.rfft(segment) _, _, stft_result = stft(segment, fs=1000, nperseg=64) dwt_result = pywt.wavedec(segment, self.wavelet, level=3) return {'fft': fft_result, 'stft': stft_result, 'dwt': dwt_result} return None

6.3 混合分析方法探索

结合STFT和DWT优势：

def hybrid_analysis(signal, sample_rate): # 第一级：STFT粗分析 f, t, Zxx = stft_analysis(signal, sample_rate, nperseg=128) # 识别感兴趣区域 energy = np.sum(np.abs(Zxx), axis=0) roi = t[energy > 0.7*np.max(energy)] # 第二级：对关键区域精细小波分析 if len(roi) > 0: start_idx = int(roi[0]*sample_rate) end_idx = int(roi[-1]*sample_rate)+1 segment = signal[start_idx:end_idx] # 执行多级小波分解 coeffs = pywt.wavedec(segment, 'sym8', level=6) return {'stft': (f, t, Zxx), 'dwt_segment': coeffs, 'roi': (roi[0], roi[-1])} return {'stft': (f, t, Zxx)}

这种方法在机械故障诊断中特别有效，先用STFT定位异常时间段，再用DWT进行精细分析。