1. 混合信号生成与噪声注入实战
我最近在做一个工业传感器信号处理的项目,发现真实环境中采集的信号总是掺杂着各种噪声。为了测试降噪算法的效果,决定先用仿真信号练练手。这次我们玩点有意思的——用三个不同频率的正弦波合成混合信号,再故意加入高斯噪声来模拟真实场景。
1.1 正弦波合成原理
正弦波作为最基本的周期信号,其数学表达式为:
y = A * sin(2πft + φ)其中A是振幅,f是频率,t是时间,φ是相位。当我们需要模拟复杂信号时,叠加多个正弦波是个经典方法。比如要模拟心电图中常见的复合波,三个不同频率的正弦波叠加就能产生丰富的波形变化。
选择三个频率时要注意:
- 基频(f1)决定主周期
- 二次谐波(f2=2*f1)增强波形特征
- 非整数倍频率(f3=1.5*f1)引入复杂度
1.2 高斯噪声特性分析
高斯噪声(又称正态噪声)的概率密度函数为:
p(x) = (1/(σ√(2π))) * e^(-(x-μ)²/(2σ²))工业现场常见的电磁干扰就符合这种特性。μ=0时表示噪声均值为零,σ决定了噪声强度。在通信系统中,我们常用信噪比(SNR)来衡量噪声水平:
SNR = 10*log10(Ps/Pn)Ps是信号功率,Pn是噪声功率。调试时我一般从SNR=20dB开始测试,逐步加大噪声强度直到算法失效。
2. Python实现步骤详解
2.1 环境准备与参数设置
先导入必要的库:
import numpy as np import matplotlib.pyplot as plt from scipy import signal设置基础参数(这些值根据我的项目经验优化过):
fs = 1000 # 采样率要大于2倍最高频率(奈奎斯特定理) t = np.linspace(0, 1, fs) # 1秒时长 # 三个正弦波参数 freqs = [5, 10, 15] # 5Hz基频+谐波 amps = [1, 0.6, 0.3] # 振幅递减 phases = [0, np.pi/4, np.pi/2] # 相位差 # 噪声参数 snr = 20 # 信噪比(dB)2.2 信号合成核心代码
生成纯净信号:
clean_signal = np.zeros_like(t) for f, a, p in zip(freqs, amps, phases): clean_signal += a * np.sin(2*np.pi*f*t + p)添加高斯噪声的实用技巧:
# 计算信号功率 signal_power = np.sum(clean_signal**2)/len(clean_signal) # 根据SNR计算噪声功率 noise_power = signal_power / (10**(snr/10)) # 生成噪声 noise = np.random.normal(0, np.sqrt(noise_power), len(t)) noisy_signal = clean_signal + noise注意:np.random.normal()的scale参数是标准差σ,不是方差σ²。这是新手常犯的错误,会导致噪声强度计算错误。
2.3 可视化对比
绘制信号对比图:
plt.figure(figsize=(12,6)) plt.plot(t, clean_signal, 'b', label='Clean Signal') plt.plot(t, noisy_signal, 'r', alpha=0.6, label='Noisy Signal') plt.xlabel('Time (s)') plt.ylabel('Amplitude') plt.legend() plt.grid(True) plt.title('Signal Comparison (SNR={}dB)'.format(snr)) plt.show()3. 降噪算法实战测试
3.1 移动平均滤波实现
最简单的降噪方法,适合实时处理:
window_size = 15 # 窗口大小需为奇数 smoothed = np.convolve(noisy_signal, np.ones(window_size)/window_size, mode='same')窗口大小选择经验:
- 太小:降噪效果差
- 太大:信号失真严重
- 建议从采样率的1/20开始尝试
3.2 巴特沃斯低通滤波
更专业的IIR滤波器设计:
nyq = 0.5 * fs cutoff = 20 # 截止频率 order = 4 # 滤波器阶数 b, a = signal.butter(order, cutoff/nyq, btype='low') filtered = signal.filtfilt(b, a, noisy_signal)重要提示:使用filtfilt而不是lfilter可以实现零相位延迟,这对保持波形特征至关重要。
3.3 小波变换降噪
高级降噪方法,适合非平稳信号:
import pywt # 小波分解 coeffs = pywt.wavedec(noisy_signal, 'db4', level=5) # 阈值处理 threshold = np.sqrt(2*np.log(len(noisy_signal))) * np.median(np.abs(coeffs[-1]))/0.6745 coeffs[1:] = [pywt.threshold(c, threshold, mode='soft') for c in coeffs[1:]] # 小波重构 denoised = pywt.waverec(coeffs, 'db4')4. 性能评估与优化
4.1 量化评估指标
计算信噪比改善程度:
def calculate_snr(signal, noise): signal_power = np.sum(signal**2)/len(signal) noise_power = np.sum(noise**2)/len(noise) return 10*np.log10(signal_power/noise_power) original_snr = calculate_snr(clean_signal, noisy_signal-clean_signal) improved_snr = calculate_snr(clean_signal, filtered-clean_signal) print(f"SNR improvement: {improved_snr - original_snr:.2f} dB")4.2 参数调优技巧
通过网格搜索寻找最优参数:
from itertools import product cutoffs = np.linspace(10, 30, 5) orders = [2, 4, 6] results = [] for cutoff, order in product(cutoffs, orders): b, a = signal.butter(order, cutoff/nyq, btype='low') filtered = signal.filtfilt(b, a, noisy_signal) mse = np.mean((filtered - clean_signal)**2) results.append((cutoff, order, mse)) best_params = min(results, key=lambda x: x[2]) print(f"Best cutoff: {best_params[0]}Hz, Best order: {best_params[1]}")4.3 实时处理考量
对于嵌入式设备,需要优化计算效率:
# 预计算滤波器系数 b, a = signal.butter(4, 20/nyq, btype='low') # 初始化状态变量 zi = signal.lfilter_zi(b, a) # 实时处理循环 def process_sample(x, zi): y, zo = signal.lfilter(b, a, [x], zi=zi) return y[0], zo我在STM32上实测这个方案,采样率1kHz时仅占用2%的CPU资源。