ch3_9_3.m

用MATLAB作的一个滤波程序

clear all
%心电图时域图
figure(1)
x1=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84,-90];
x2=[-66,-32,-4,-2,8,12,12,10,6,6,4,0,0,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0];
x=[x1,x2];
subplot(221);
p=0:length(x)-1;
plot(p,x);
title('心电图时域图')
axis([0 50 -100 20])
grid;
%心电图频域图
subplot(222)
fs=1000  %设采样频率1KHZ
k=0:1023
w=-500+k*1000/length(k)
y=fft(x,1024);
y=fftshift(abs(y));
plot(w,y);
title('心电图频域图')
axis([-200 200 0 400])
grid;

%设计IIR低通滤波器滤掉大于50HZ的频谱
Ap=1;As=60;

[N,Wc]=buttord(0.2,0.3,1,60);
disp('滤波器的阶数');
fprintf('%d\n',N);
[num,den]=butter(N,Wc,'low');
[h,w2]=freqz(num,den,256,1000);
subplot(223)
plot(w2,abs(h));
xlabel('normalized frequency');
ylabel('magnitude')
axis([0 200 0 1.2])
title('BW数字低通')
grid;

%设计FIR低通滤波器滤掉大于50HZ的频谱
devp = (10^(Ap/20)-1)/(10^(Ap/20)+1); 
devr = 10^(-As/20);
dev = [devp,devr];
a=[1 0];
f=[100 148];
[n,fo,ao,w]=remezord(f,a,dev,1000);
b=remez(n,fo,ao,w);
[h,w]=freqz(b,1,256,1000);
subplot(224)
plot(w,abs(h));
xlabel('normalized frequency');
ylabel('magnitude')
axis([0 150 0 1.1]);
title('优化设计法数字低通');
grid;

%用IIR低通滤波
figure(2)
z=filter(b,1,x)

%低通滤波后信号的时域图
subplot(221);
plot(p,z)
axis([0,50,-100,20]);
xlabel('time ');
ylabel('amplify');
title('低通滤波后信号的时域图')
grid;

%低通滤波后信号的频谱图

fl=fft(z,1024);
L=0:1023;
p=length(L);
X=fftshift(fl);
w1=-500+(0:p-1)*2*500./p;
subplot(222);
plot(w1,abs(X));
xlabel(' frequency  HZ');
ylabel('magnitude  ');
axis([-500 500 0 400]);
title('低通滤波后用FFT得到幅频响应')
grid;

%功率谱估计原始信号x[k]
subplot(223);
[pxx,f]=periodogram(x)
plot(f,pxx)
title('功率谱估计原始信号x[k]')
grid;

%功率谱估计滤波后的信号y[k]
subplot(224);
[pxx,f]=periodogram(z)
plot(f,pxx)
title('功率谱估计滤波后的信号y[k]')
grid;