ch3_9_12.m

用MATLAB作的一个滤波程序

f0=10;
Fs=8000
k=0:2000 
ts=k/Fs
w0=2*pi*f0;
y=square(w0*ts, 50)+randn(1, length(k)); %duty cycle=50%
figure(1)
subplot(221);
plot(ts,y);
title('方波时域图')
xlabel('time ');
ylabel('amplify');
axis([0,0.2,-3,3]);

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

%设计滤波器滤掉方波中大于400HZ的频率
[N,Wc]=cheb1ord(350/4000,450/4000,1,60);
[numl,denl]=cheby1(N,1,Wc);
[h,w]=freqz(numl,denl,1024,8000);
subplot(223);
plot(w,20*log10(abs(h)))
xlabel(' frequency  HZ');
ylabel('magnitude  DB');
title('cheby1数字低通滤波器')
axis([200 600 -100 1])
grid;

%滤波后方波的频谱
z=filter(numl,denl,y)
fl=fft(z,1024);
L=0:1023;
p=length(L);
X=fftshift(fl);
w1=-4000+(0:p-1)*2*4000./p;
subplot(224);
plot(w1,abs(X));
xlabel(' frequency  HZ');
ylabel('magnitude  ');
axis([-500 500 0 200]);
title('低通滤波后用FFT得到幅频响应')

%滤波后方波的时域图
figure(2)
subplot(221);
plot(ts,z)
axis([0,0.2,-3,3]);
xlabel('time ');
ylabel('amplify');
title('滤波后方波的时域图')
grid;