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

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

1、概述

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

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

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

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

2、实例

例、有一个心电信号，但测量时由于设备的问题受到工频信号的干扰，并完全被干扰所淹没。设计一个数字陷波器来恢复心电信号。

心电信号的数据在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');

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

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

3、讨论

       有时，在工频干扰中不限于只有基频 50Hz，往往还带有它的谐波。所以在处理之前先做频谱分析，观察除了基波外还有几个谐波。除基波的陷波器外，还可以设计谐波的陷波器，把它们级联在一起，以消除基波和谐波。

      有时处理的不是 50Hz 工频信号，而是在信号中混有其他信号的正弦信号。这时先要做频谱分析，甚至对干扰的正弦信号要用校正法求出该信号的频率，才能进行陷波处理。

========================

从信号中去除 60Hz 杂声_OliverH-yishuihan的博客-CSDN博客

========================

====================另外一种方式===========================================

原文链接：https://blog.csdn.net/limindaihong/article/details/9205211

————————————————

基于matlab的数字陷波器设计

陷波器是一种简单的二阶IIR滤波器，其幅度响应在某一频率上为零，可用来消除某个频率分量，如：滤除信号中由电源引起的50Hz工频干扰。其系统函数为：

其中： ----陷波数字频率(rad)； f0----陷波频率(Hz)；fs----取样频率(Hz)；

r----常数。

实验要求：编程实现以下功能：

1）设陷波频率 f0=50Hz，取样频率 fs=600Hz，r=0.9，画H(z)的幅频和相频特性。

2）画H(z)的零极点图，体会陷波原理。

3）利用该陷波器对信号：， (n=0~599)进行滤波，画出 x(n) 及滤波输出 y(n)。

程序代码：

%数字陷波器
clear
clc
f0 = 50;
fs = 600;
r = 0.9;
w0 = 2*pi*f0/fs;
b = [1 -2*cos(w0) 1];
a = [1 -2*r*cos(w0) r*r];
N = 1024;
[H,w] = freqz(b,a,N);
subplot(221); plot(w,abs(H)); grid on; title('陷波器的幅频响应');
subplot(222); plot(w,angle(H)); grid on; title('陷波器的相频响应');
subplot(223); zplane(b,a); grid on; title('陷波器的零极点图');

n = 0:N-1;
x = sin(2*pi*50*n/fs)+sin(2*pi*100*n/fs);
X = fft(x,N);
y = filter(b,a,x);
Y = fft(y,N);
f = fs/N*(0:N/2-1);
figure;
subplot(221); plot(n,x); grid on; title('原信号x(n)');
subplot(222); plot(f,abs(X(1:N/2))); grid on; title('x(n)的幅频谱');
subplot(223); plot(n,y); grid on; title('陷波器滤波后的信号y(n)');
subplot(224); plot(f,abs(Y(1:N/2))); grid on; title('y(n)的幅频谱');

程序运行结果图：

————————————————

原文链接：https://blog.csdn.net/limindaihong/article/details/9205211

===============================================================