三亩地 三亩地SAN MU DI · CODE DIARY
ARTICLE DETAIL

日记详情

真实记录编程学习的某一天,欢迎挑你感兴趣的翻一翻。

基于傅立叶变换的时序信号去噪实战:从理论到Python实现

基于傅立叶变换的时序信号去噪实战:从理论到Python实现

1. 傅立叶变换:时域与频域的桥梁

第一次接触傅立叶变换时,我也被那一堆数学符号吓到了。直到有天调试传感器数据时,发现传统方法怎么也滤不掉周期性干扰,才真正体会到这个工具的威力。简单来说,它就像给声音装了个"频率显微镜"——把杂乱的时间波形拆解成不同频率的正弦波组合。

想象你在听交响乐录音,所有乐器声音混在一起。傅立叶变换就是那个能把小提琴、大提琴声部分离出来的神奇工具。对于50Hz工频干扰的脑电信号,或是120Hz电机噪声的振动数据,这个特性尤其有用。实际项目中,我常用它处理三类问题:

  • 周期性噪声滤除(如工业设备振动分析)
  • 特征频率提取(如语音识别中的共振峰)
  • 数据压缩(保留主要频率成分)
import numpy as np from scipy.fft import rfft, rfftfreq # 生成含噪信号示例 t = np.linspace(0, 1, 1000) clean = np.sin(2*np.pi*50*t) + 0.5*np.sin(2*np.pi*120*t) noisy = clean + np.random.normal(0, 1, len(t))

2. 实战准备:Python环境与数据模拟

2.1 工具链配置

推荐使用Anaconda创建专属环境:

conda create -n fourier python=3.8 conda install numpy scipy matplotlib

最近发现Scipy的fft模块比Numpy版本更友好,特别是rfft()系列函数会自动忽略负频率镜像,节省了50%计算量。对于超长时序数据,可以试试scipy.fftpack.next_fast_len优化数组长度。

2.2 构造仿真数据

模拟工业场景中常见的"信号+噪声"组合:

def generate_signal(duration=1, sample_rate=1000): t = np.linspace(0, duration, int(duration*sample_rate)) # 真实信号成分 main_component = 2*np.sin(2*np.pi*50*t) # 50Hz主频 harmonic = 0.8*np.sin(2*np.pi*150*t) # 三次谐波 trend = 0.5*t # 线性趋势项 # 噪声成分 random_noise = np.random.normal(0, 1, len(t)) periodic_noise = 1.2*np.sin(2*np.pi*30*t) # 设备振动噪声 return t, main_component + harmonic + trend + random_noise + periodic_noise

这种构造方式比单纯用白噪声更接近真实场景,包含:

  • 主频信号(50Hz)
  • 谐波干扰(150Hz)
  • 低频趋势项
  • 随机噪声
  • 特定频率干扰(30Hz)

3. 频域分析四步法

3.1 快速傅立叶变换实施

使用scipy.fft.rfft获取频谱时,有个容易踩的坑——频率轴的计算:

def compute_fft(signal, sample_rate): n = len(signal) yf = rfft(signal) xf = rfftfreq(n, 1/sample_rate) # 幅度谱修正 magnitude = 2/n * np.abs(yf) return xf, magnitude

注意幅度谱需要做2/N的缩放,这个细节很多教程会忽略。我曾因此浪费半天时间调试为什么幅值对不上。

3.2 频谱诊断技巧

健康的频谱通常呈现几个特征峰+平坦噪声基底。通过观察可以识别:

  1. 突出尖峰 → 主信号成分
  2. 矮胖凸起 → 随机噪声
  3. 规律间隔峰 → 谐波干扰

用Matplotlib做对比可视化:

plt.figure(figsize=(12,6)) plt.plot(xf, magnitude, 'r', label='含噪信号') plt.plot(xf_clean, magnitude_clean, 'b--', label='纯净信号') plt.axvline(50, color='k', linestyle=':', alpha=0.5) plt.axvline(120, color='k', linestyle=':', alpha=0.5) plt.legend()

3.3 动态阈值降噪法

固定阈值在工程中往往不够鲁棒,我改进的滑动窗口阈值算法:

def adaptive_threshold(spectrum, window_size=5): threshold = np.zeros_like(spectrum) for i in range(len(spectrum)): start = max(0, i-window_size//2) end = min(len(spectrum), i+window_size//2+1) local_median = np.median(spectrum[start:end]) threshold[i] = local_median * 3 # 3倍中值作为阈值 return threshold

相比全局阈值,这种方法能更好保留弱信号成分,特别适合信噪比不稳定的场景。

4. 进阶实战:非平稳信号处理

4.1 短时傅立叶变换应用

当信号频率随时间变化时(如变频电机振动),需要scipy.signal.stft:

from scipy.signal import stft f, t, Zxx = stft(signal, fs=sample_rate, nperseg=256) plt.pcolormesh(t, f, np.abs(Zxx), shading='gouraud')

关键参数nperseg控制时间分辨率,通常取2的整数幂。值越小时间分辨率越高,但频率分辨率越低。

4.2 滤波器组实现

对于需要保留多个频段的场景,可以设计带通滤波器组:

from scipy.signal import butter, filtfilt def bandpass_filter(data, lowcut, highcut, fs, order=5): nyq = 0.5 * fs low = lowcut / nyq high = highcut / nyq b, a = butter(order, [low, high], btype='band') return filtfilt(b, a, data)

使用filtfilt实现零相位延迟,比普通lfilter更适合时序分析。

5. 工程经验与调试技巧

5.1 常见问题排查

  • 频谱泄漏:记得加窗函数(汉宁窗效果较好)
window = np.hanning(len(signal)) yf = rfft(signal * window)
  • 频率混叠:确保采样率>2倍最高频率
  • 幅值失真:检查是否做了正确的幅度缩放

5.2 性能优化

处理百万级数据点时:

  1. 使用scipy.fft.next_fast_len优化数组长度
  2. 考虑pyfftw替代方案
  3. 分段处理+重叠保留法
from scipy.fftpack import next_fast_len n_optimized = next_fast_len(len(signal)) yf = rfft(signal, n=n_optimized)

6. 扩展应用案例

6.1 工业振动分析

某风机轴承监测项目中,通过FFT成功分离出:

  • 25Hz的轴转频
  • 107Hz的轴承故障特征频率
  • 300Hz以上的齿轮啮合噪声

6.2 音频处理实例

降噪前后频谱对比显示,300Hz以下的环境噪声被有效抑制:

# 音乐片段处理示例 import librosa y, sr = librosa.load('noisy_audio.wav', sr=None) D = librosa.stft(y) magnitude = np.abs(D) threshold = np.median(magnitude) * 2 D_clean = D * (magnitude > threshold) y_clean = librosa.istft(D_clean)

在最近的心电信号处理项目中,傅立叶变换帮我定位到了50Hz工频干扰和0.5Hz的呼吸伪迹。有意思的是,当尝试用传统滤波器处理时,要么残留噪声,要么损伤QRS波,而频域处理则完美平衡了这两者。这让我想起一位导师说过:"当你卡在时域问题时,不妨换个角度看看频域"。

← 返回列表