ch3_9_1.m

来自「用MATLAB作的一个滤波程序」· M 代码 · 共 76 行

%生成频率10HZ，占空比为50%的周期方波
f0=10;
Fs=1000
k=0:1000 
ts=k/Fs
w0=2*pi*f0;
y = square(w0*ts, 50); %duty cycle=50%
figure(1)
subplot(221);
plot(ts,y);
title('方波时域图')
xlabel('time ');
ylabel('amplify');
axis([0,0.5,-1.2,1.2]);

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

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

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