%%%%%%%%%% 数据读取
ECG=load('1record.dat');
l=length(ECG);
t=[150/l:150/l:30]
x=ECG(4,:); %可取不同的组
y=x(1:5:end)
figure(1)
plot(t, 20*log(abs(y)));
title('含噪心电信号');
xlabel('时间(s)');
ylabel('幅度(dB)');
axis([0 15 179 183])
%%%%%%%%%% 频域分析
n=4000
m=abs(fft(y,n));
fs=200;
f=fs/n*(0:n-1);
figure(2)
plot(f, m);
title('心电信号的频谱图');
xlabel('频率 f/Hz');ylabel('幅值/db');
axis([0 100 -100 40000])
%%%%%%%%%% 加入工频噪声 时域及频域分析
x=20*sin(2*50*pi*t);
y1=y+x;
figure(3);
subplot(211);
plot(t, 20*log(abs(y1)));
axis([0 15 179 183])
title('加工频噪声时域波形');
xlabel('时间t/s');ylabel('幅值/db');
k=abs(fft(y1,n));
subplot(212);
plot(f,k);
axis([0 100 0 40000])
title('加工频噪声频谱图');
xlabel('频率 /Hz');ylabel('幅值')
%%%%%%%%%% 工频噪声的滤除 巴特沃思
fs=200;
wp=[47 53]*2/fs;
ws=[42 58]*2/fs;
Rp=3;Rs=45;
[N,wn]=buttord(wp,ws,Rp,Rs) % N=6
[b,a]=butter(N,wn,'stop')
figure(4)
[H,W]=freqz(b,a,1024);%生成滤波器的幅频响应
k=0:1023;
plot((fs/2)/1024 *k,abs(H));%输出滤波器的频率响应
axis([0 100 0 1.2])
grid on
title('巴特沃思陷波器频率响应')
y1n=filter(b,a, y1); %滤波
Yk=fft(y1n,n);
Ykl=abs(Yk)/n;
figure(5)
subplot(2,1,1)
plot(t, 20*log(abs(y1n)));
axis([0 20 179 183])
title('波形通过陷波器(巴特沃思)后的波形');
subplot(2,1,2)
plot(f, 20*log(abs(Ykl)));
title('波形通过陷波器(巴特沃思)后的频谱')
axis([0 100 0 100])
%%%%%%%%%% 工频噪声的滤除 整系数滤波器 每隔50Hz
b=[1 1 1 1];
figure(6)
freqz(b,1);
title('整系数陷波器频率响应');
y1n=filter(b,1, y1); %滤波
Yk=fft(y1n,n);
Ykl=abs(Yk)/n;
figure(7)
subplot(2,1,1)
plot(t, 20*log(abs(y1n)));
axis([0 20 206 210])
title('波形通过整系数陷波器后的波形');
subplot(2,1,2)
plot(f, 20*log(abs(Ykl)));
title('波形通过整系数陷波器后的频谱')
axis([0 100 0 100])
%%%%%%%%%% 工频噪声的滤除 切比雪夫1 IIR
fs=200;
wp=[47 53]*2/fs;
ws=[42 58]*2/fs;
Rp=3;Rs=45;
[N,wn]= cheb1ord(wp,ws,Rp,Rs) % N=4
[b,a]=cheby1(N,Rp,wn,'stop')
figure(8)
[H,W]=freqz(b,a,1024);%生成滤波器的幅频响应
k=0:1023;
plot((fs/2)/1024 *k,abs(H));%输出滤波器的频率响应
axis([0 100 0 1.2])
grid on
title('切比雪夫1陷波器频率响应')
y1n=filter(b,a, y1); %滤波
Yk=fft(y1n,n);
Ykl=abs(Yk)/n;
figure(9)
subplot(2,1,1)
plot(t, 20*log(abs(y1n)));
axis([0 20 172 175])
title('波形通过切比雪夫1陷波器后的波形');
subplot(2,1,2)
plot(f, 20*log(abs(Ykl)));
title('波形通过切比雪夫1陷波器后的频谱')
axis([0 100 0 100])
%%%%%%%%%% 工频噪声的滤除 切比雪夫II IIR
fs=200;
wp=[47 53]*2/fs;
ws=[42 58]*2/fs;
Rp=3;Rs=45;
[N,wn]= cheb2ord(wp,ws,Rp,Rs) % N=4
[b,a]=cheby2(N,Rs ,wn,'stop')
figure(10)
[H,W]=freqz(b,a,1024);%生成滤波器的幅频响应
k=0:1023;
plot((fs/2)/1024 *k,abs(H));%输出滤波器的频率响应
axis([0 100 0 1.2])
grid on
title('切比雪夫II陷波器频率响应')
y1n=filter(b,a, y1); %滤波
Yk=fft(y1n,n);
Ykl=abs(Yk)/n;
figure(11)
subplot(2,1,1)
plot(t, 20*log(abs(y1n)));
axis([0 20 178 182])
title('波形通过切比雪夫II陷波器后的波形');
subplot(2,1,2)
plot(f, 20*log(abs(Ykl)));
title('波形通过切比雪夫II陷波器后的频谱')
axis([0 100 0 100])
%%%%%%%%%% 工频噪声的滤除 椭圆滤波器 IIR
fs=200;
wp=[47 53]*2/fs;
ws=[42 58]*2/fs;
Rp=3;Rs=45;
[N,wn]= ellipord(wp,ws,Rp,Rs) % N=3
[b,a]=ellip (N,Rp,Rs ,wn,'stop')
figure(12)
[H,W]=freqz(b,a,1024);%生成滤波器的幅频响应
k=0:1023;
plot((fs/2)/1024 *k,abs(H));%输出滤波器的频率响应
title('椭圆陷波器频率响应')
axis([0 100 0 1.2])
grid on
y1n=filter(b,a, y1); %滤波
Yk=fft(y1n,n);
Ykl=abs(Yk)/n;
figure(13)
subplot(2,1,1)
plot(t, 20*log(abs(y1n)));
axis([0 20 178 182])
title('波形通过椭圆陷波器后的波形');
subplot(2,1,2)
plot(f, 20*log(abs(Ykl)));
title('波形通过椭圆陷波器后的频谱')
axis([0 100 0 100])
%%%%%%%%%% 滤除工频干扰 FIR 布莱克曼窗
wp1=0.46;wp2=0.54;
ws1=0.49;ws2=0.51;
wp=[wp1 wp2]*pi;
ws=[ws1 ws2]*pi;
wdel=min((ws1-wp1),( wp2-ws2));
N=ceil(8*pi/wdel); % N=838
Wn=(wp+ws) /2;
window=blackman (N+1);
b=fir1(N,Wn/pi, 'stop' ,window);
figure(14)
freqz(b,1,1024)
title('FIR陷波器频率响应')
y1n=filter(b,1, y1); %滤波
Yk=fft(y1n,n);
Ykl=abs(Yk)/n;
figure(15)
subplot(2,1,1)
plot(t, 20*log(abs(y1n)));
axis([0 20 178 182])
title('波形通过FIR陷波器后的波形');
subplot(2,1,2)
plot(f, 20*log(abs(Ykl)));
title('波形通过FIR陷波器后的频谱')
axis([0 100 0 100])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% 去除基线漂移 巴特沃斯
fs=200;
wp=1*2/fs;
ws=0.1 *2/fs;
rp=3;
rs=45;
n=4000;
[N,wn]=buttord(wp,ws,rp,rs); % N=3
[b,a]=butter(N,wn,'high');
figure(16)
freqz(b,a,n,fs)
title('高通滤波器(巴特沃思)');
axis([0 100 -40 40 ])
yn=filter(b,a, y); %滤波
Yk=fft(yn,n);
Ykl=abs(Yk)/n;
figure(17)
subplot(2,1,1)
plot(t, yn);
axis([0 30 -600 600])
title('波形通过高通滤波器(巴特沃思)后的波形');
subplot(2,1,2)
plot(f, Ykl);
title('波形通过高通滤波器(巴特沃思)后的频谱')
axis([0 100 0 90])
%%%%%%%%%% 去除基线漂移 简单极点法
b=[0.9876 -0.9876];
a=[1 -0.9752];
figure(18)
freqz(b,a,4000,200)
title('基于简单极点法的高通滤波器频率响应')
yn=filter(b,a, y); %滤波
Yk=fft(yn,n);
Ykl=abs(Yk)/n;
figure(19)
subplot(2,1,1)
plot(t, yn);
axis([0 30 -600 600])
title('波形通过高通滤波器后的波形');
subplot(2,1,2)
plot(f, Ykl);
title('波形通过高通滤波器后的频谱')
axis([0 100 0 90])
%%%%%%%%%% 去除基线漂移 切比雪夫1 IIR
fs=200;
wp=1*2/fs;
ws=0.1 *2/fs;
rp=3;
rs=45;
n=4000;
[N,wn]=cheb1ord (wp,ws,rp,rs); % N=2
[b,a]=cheby1(N,rp,wn,'high');
figure(20)
freqz(b,a,n,fs)
title('切比雪夫1高通滤波器频率响应');
axis([0 100 -40 40 ])
yn=filter(b,a, y); %滤波
Yk=fft(yn,n);
Ykl=abs(Yk)/n;
figure(21)
subplot(2,1,1)
plot(t, yn);
axis([0 30 -400 400])
title('波形通过切比雪夫1高通滤波器后的波形');
subplot(2,1,2)
plot(f, Ykl);
title('波形通过切比雪夫1高通滤波器后的频谱')
axis([0 100 0 60])
%%%%%%%%%% 去除基线漂移 切比雪夫II IIR
fs=200;
wp=1*2/fs;
ws=0.1 *2/fs;
rp=3;
rs=45;
n=4000;
[N,wn]=cheb2ord (wp,ws,rp,rs); % N=2
[b,a]=cheby2(N,rs,wn,'high');
figure(22)
freqz(b,a,n,fs)
title('切比雪夫II高通滤波器频率响应');
axis([0 100 -40 40 ])
yn=filter(b,a, y); %滤波
Yk=fft(yn,n);
Ykl=abs(Yk)/n;
figure(23)
subplot(2,1,1)
plot(t, yn);
axis([0 30 -600 600])
title('波形通过切比雪夫II高通滤波器后的波形');
subplot(2,1,2)
plot(f, Ykl);
title('波形通过切比雪夫II高通滤波器后的频谱')
axis([0 100 0 60])
%%%%%%%%%% 去除基线漂移 椭圆滤波器 IIR
fs=200;
wp=1*2/fs;
ws=0.1 *2/fs;
rp=3;
rs=45;
n=4000;
[N,wn]=ellipord (wp,ws,rp,rs); % N=2
[b,a]=ellip (N,rp,rs,wn,'high');
figure(24)
freqz(b,a,n,fs)
title('椭圆高通滤波器频率响应')
axis([0 100 -40 40 ])
yn=filter(b,a, y); %滤波
Yk=fft(yn,n);
Ykl=abs(Yk)/n;
figure(25)
subplot(2,1,1)
plot(t, yn);
axis([0 30 -500 500])
title('波形通过椭圆高通滤波器后的波形');
subplot(2,1,2)
plot(f, Ykl);
title('波形通过椭圆高通滤波器后的频谱')
axis([0 100 0 60])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% 功率谱计算 Welch y(初始信号,含基线漂移,不含工频)
fs=200;
t=[150/l:150/l:30]
nfft=4096;
window=hamming(512);
overlap=256;
f1=fs/nfft*(0:nfft-1);
[Pyy,f1]=pwelch(y,window,overlap,nfft,fs, 'onesided');
figure(26)
plot(f1,10*log10(Pyy));
xlabel(' Frequency(Hz) '); ylabel(' Power Spectral(db/Hz) ');
title('去噪的的心电信号的Welch功率谱估计');
%% y (初始信号,含基线漂移,无工频干扰)
%% y1 (有基线漂移,有工频干扰)
%% y1n (y1信号通过陷波器,含基线漂移)
%% yn (y信号通过高通滤波器)
- 1
- 2
- 3
前往页