使用数字陷波器滤除工频信号
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
===============================================================