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

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

• 数字陷波器函数如下函数：iirnotch
• 功能：数字陷波器设计
• 调用格式：[b,a]= iirnotch(Wo,Bw)

说明：Wo是陷波器的中心频率，Bw是陷波器的带宽，参数b和a是数字滤波器的系数。

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

例、有一个心电信号，但测量时由于设备的问题受到工频信号的干扰，

并完全被干扰所淹没。设计一个数字陷波器来恢复心电信号。

心电信号的数据在noisyecg.mat中，设计数字陷波器滤除干扰噪声。

程序如下：

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

clc; clear; close all;
 
load('iirnotch_noisyecg.mat'); % 读入信号数据和采样频率
x = noisyecg;    % 信号为x
N = length(x);   % 信号长N
t = (0:N-1)./fs; % 时间刻度 
df = fs/(N-1);   % 分辨率
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');