基于MATLAB的数字信号处理、数字滤波器设计与实现

最近在调试一个脑电信号采集设备，发现原始信号里总掺杂着50Hz工频干扰。这玩意儿不处理掉根本没法做后续分析，抄起MATLAB就开始折腾滤波器。数字滤波这事儿说简单也简单，但实际操作时各种参数设置能让人抓狂。

先看个活生生的例子。假设我们采样率1000Hz，信号里混着50Hz和120Hz成分：

fs = 1000;
t = 0:1/fs:1;
x = sin(2*pi*50*t) + 0.5*sin(2*pi*120*t);
noise = 0.8*randn(size(t));
signal = x + noise;

这时候频谱分析不能少，FFT走起：

N = length(signal);
f = (-N/2:N/2-1)*(fs/N);
Y = fft(signal);
Y_shift = fftshift(Y);
plot(f, abs(Y_shift)/N);

频谱图上明显看到50Hz和120Hz两个尖峰，中间还夹着噪声的毛刺。我们的目标是用滤波器精准干掉120Hz那个峰，同时尽量保留50Hz信号。这时候FIR滤波器就该出场了——相位特性好，不会把信号波形搞变形。

基于MATLAB的数字信号处理、数字滤波器设计与实现

设计个80Hz截止的低通滤波器试试：

order = 50;
cutoff = 80/(fs/2);
b = fir1(order, cutoff, hamming(order+1));

这里用汉明窗主要是平衡主瓣宽度和旁瓣衰减。不过实际跑起来发现阶数设50还是有点低，120Hz衰减不够。把阶数提到70，效果立竿见影：

freqz(b,1,1024,fs); % 看幅频响应曲线

但阶数太高实时处理会吃力，这时候可以换凯撒窗，用更少阶数达到类似效果：

b = fir1(40, cutoff, kaiser(41,4));

滤波操作直接上卷积：

filtered = conv(signal, b, 'same');

不过要注意边界效应，建议用filtfilt做零相位滤波：

filtered = filtfilt(b, 1, signal);

这时候再看频谱，120Hz成分基本消失了。但实际项目中更常遇到的是带阻滤波需求，比如要专门消除50Hz工频干扰：

wo = 50/(fs/2);
bw = 4/(fs/2);
[b,a] = iirnotch(wo, bw);

这个二阶IIR陷波器简单粗暴，但要注意稳定性和相位失真问题。有个坑是当信号含有多个干扰频率时，可能需要级联多个陷波器，这时候滤波顺序会影响最终效果。

最后来个实用技巧：处理ECG信号时，经常用移动平均滤波消除基线漂移：

windowSize = 20; 
b = ones(1, windowSize)/windowSize;
ecg_filtered = filter(b, 1, ecg_signal);

但这样会产生信号延迟，解决方法同样是使用filtfilt双向滤波。不过要注意，实时系统里可没法用这招，得提前想好延迟补偿方案。