notch01.m

50hz陷波器

clc
Name1='031015090000-031015090947---上升段1.ECG';
fid1=fopen(Name1);
fseek(fid1,23,'bof');
data1=fread(fid1,'uint16');
data1=data1(1:3500);
data1=data1-mean(data1);
data0=data1;%data0原始信号
fs=250;
fg=50;
k=2*pi*fg/fs;

t=1:3500;
x=100*sin(2*pi*fg*t/250);
data1=data1+x'; %data1加工频干扰信号

b=[1 -2*cos(k) 1];
a=[1 -2*cos(k)*0.99 0.98];
[H,w]=freqz(b,a,1000,'whole');
H=(H(1:501))'; w=(w(1:501))';
mag=abs(H);
db=20*log10((mag+eps)/max(mag));
pha=angle(H);
grd=grpdelay(b,a,w);

y=filter(b,a,data1);
figure(1)
SUBPLOT(3,1,1),plot(data0)
SUBPLOT(3,1,2),plot(data1)
SUBPLOT(3,1,3),plot(y)

figure(2)
%freqz(b,a)
plot(w*125/pi, db)
grid

figure(3)
%功率谱密度计算

Y0=fft(data0,512);
P0=Y0.*conj(Y0)/512;
f0=250*(0:255)/512;

Y1=fft(data1,512);
P1=Y1.*conj(Y1)/512;
f1=250*(0:255)/512;

Y2=fft(y,512);
P2=Y2.*conj(Y2)/512;
f2=250*(0:255)/512;

SUBPLOT(3,1,1),plot(f0,P0(1:256))
SUBPLOT(3,1,2),plot(f1,P1(1:256))
SUBPLOT(3,1,3),plot(f2,P2(1:256))