使用数字陷波器滤除工频信号

在实际测量时经常会受到工频信号（交流50Hz）的干扰，有时干扰还很大，有用信号完全被淹没了。可以应用数字陷波器来消除工频信号的干扰。

数字陷波器函数如下：

函数：iirnotch
功能：数字陷波器设计
调用格式：
[b,a]= iirnotch(Wo,Bw)
说明：Wo是陷波器的中心频率，Bw是陷波器的带宽，参数b和a是数字滤波器的系数。

案例、有一个心电信号，但测量时由于设备的问题受到工频信号的干扰，并完全被干扰所淹没。设计一个数字陷波器来恢复心电信号。
心电信号的数据在noisyecg.mat中，设计数字陷波器滤除干扰噪声。程序如下：

clear all; clc; close all;

load('noisyecg.mat');              % 读入信号数据和采样频率
x=noisyecg;                        % 信号为x
N=length(x);                       % 信号长N
t =(0:N-1)./fs;                    % 时间刻度 
fs2=fs/2;                          % 设置奈奎斯特频率
W0=50/fs2;                         % 陷波器中心频率
BW=0.1;                            % 陷波器带宽 
[b,a]=iirnotch(W0,BW);             % 设计IIR数字陷波器
[H,w]=freqz(b,a);                  % 求出滤波器的频域响应
y=filter(b,a,x);                   % 对信号滤波
% 作图
plot(w*fs/2/pi,20*log10(abs(H)),'k');
xlabel('频率/Hs'); ylabel('幅值/dB');
title('陷波器幅值响应图')
axis([0 125 -45 5]); grid;
set(gca,'XTickMode','manual','XTick',[0,40,50,60,100])
set(gca,'YTickMode','manual','YTick',[-40,-30,-20,-10,0])
set(gcf,'color','w');
figure(2)
subplot 211; plot(t,x,'k','linewidth',2);
xlabel('时间/s'); ylabel('幅值');
title('带噪心电波形图')
subplot 212; plot(t,y,'k');
xlabel('时间/s'); ylabel('幅值');
title('消噪后心电波形图')
set(gcf,'color','w');

运行结果如下：

 

说明：在调用irnotch函数时，Wo和Bw参数都是归一化的频率值，在0～1之间。因为要滤除的成分是50Hz，所以归一化时用奈奎斯频率进行归一化处理。

从图中可以看到，在滤波前，带噪信号中的噪声完全淹没了心电信号，除了几个峰值外什么也看不出来；但在滤波后恢复了心电信号，但在初始部分还是有瞬态现象存在。

有时，在工频干扰中不限于只有基频50Hz，往往还带有它的谐波。所以在处理之前先做谱分析，观察除基波外还有几个谐波。除基波的陷波器外，还可设计谐波的陷波器，把它们级联在一起，以滤除基波和谐波。有时处理的不是50Hz工频信号，而是在信号中混有其他的正弦信号。这时先要做谱分析，甚至对干扰的正弦信号要用校正法求出该信号的频率，才能进一步进行陷波处理。

实验数据链接如下：

https://mp.csdn.net/mp_download/manage/download/UpDetailed

参考文献：MATLAB数字信号处理85个实用案例精讲——入门到进阶；宋知用（编著）