ch3_9_2.m

用MATLAB作的一个滤波程序

clear all;
%信号的时域及频域分析
Fs=200
k=0:1000 
ts=k/Fs
y=sin(2*pi*5*ts)+sin(2*pi*15*ts)+sin(2*pi*30*ts)
subplot(221)
plot(ts,y);
title('时域图')
xlabel('time ');
ylabel('amplify');
axis([0,1,-3,3]);
grid;

%频谱变换
dfc=fft(y,1024);
L=0:1023;
p=length(L);
X=fftshift(dfc);
w1=-100+(0:p-1)*2*100/p;
subplot(222);
plot(w1,abs(X));
xlabel('frequency');ylabel('magnitude');
axis([-40 40 0 300])
title('DFT变换频谱图')
grid;

%设计IIR高通滤波器滤掉5HZ的频谱
Ap=1;As=60;

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

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

%用IIR高通滤波
figure(2)
z=filter(num,den,y)

%高通滤波后信号的时域图
subplot(221);
plot(ts,z)
axis([0,1,-3,3]);
xlabel('time ');
ylabel('amplify');
title('高通滤波后信号的时域图')
grid;

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

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

g=filter(b,1,z)

%低通滤波后信号的时域图
subplot(223);
plot(ts,g)
axis([0,1,-3,3]);
xlabel('time ');
ylabel('amplify');
title('低通滤波后信号的时域图')
grid;

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

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