（1）生成一个基频为10Hz的周期方波信号；
（2）设计一个滤波器，滤去该周期信号中40Hz以后的频率成分，观察滤波前后的信号波形和频谱。
%采样频率取200Hz，可以改变看看。输入、输出信号频谱同时绘制在一起对比
%采用带通滤波器，只保留40Hz频率成分
clear
fs=200;
t=0:1/fs:1;
x=square(2*pi*10*t);
wp=[35 45]*2/fs;
ws=[30 50]*2/fs;
Rp=3;Rs=45;Nn=128;
[N,wn]=buttord(wp,ws,Rp,Rs)
[b,a]=butter(N,wn,'bandpass') %只保留40Hz频率成分
y=filter(b,a,x);
figure(1)
plot(t,x,'r-',t,y)
grid on
axis([0 1.2 -1.2 1.2])
title('红色代表原信号，蓝色代表只保留40Hz频率成分')
figure(2)
[H,W]=freqz(b,a);
k=0:511;
plot((fs/2)/512*k,abs(H));
grid on
title('滤波器频率响应')
T=1/fs;
N=4*(fs/10);
n=0:N-1;
xn=square(2*pi*10*n*T);
X=fftshift(fft(xn,512));
xk=1/N*X;
Y=fftshift(fft(y,512));
yk=1/N*Y;
figure(3)
plot(-fs/2+fs/512*k,abs(xk),'b',-fs/2+fs/512*k,abs(yk),'r')
grid on
legend('原信号的频谱','滤波后信号的频谱')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%采样频率取200Hz，可以改变看看。输入、输出信号频谱同时绘制在一起对比
%采用低通滤波器，保留40Hz以下频率成分
clear
fs=200;
t=0:1/fs:1;
x=square(2*pi*10*t);
wp=40*2/fs;
ws=45*2/fs;
Rp=3;Rs=45;Nn=128;
[N,wn]=buttord(wp,ws,Rp,Rs)
[b,a]=butter(N,wn,'low')
y=filter(b,a,x);
figure(1)
plot(t,x,'r-',t,y)
%grid on
axis([0 1.2 -2 2])
title('红色代表原信号，蓝色代表只保留40Hz以下频率成分')
figure(2)
[H,W]=freqz(b,a);
k=0:511;
plot((fs/2)/512*k,abs(H));
grid on
title('滤波器频率响应')
T=1/fs;
N=4*(fs/10);
n=0:N-1;
xn=square(2*pi*10*n*T);
X=fftshift(fft(xn,512));
xk=1/N*X;
Y=fftshift(fft(y,512));
yk=1/N*Y;
figure(3)
plot(-fs/2+fs/512*k,abs(xk))
grid on
legend('方波信号的频谱')
figure(4)
plot(-fs/2+fs/512*k,abs(xk),'b',-fs/2+fs/512*k,abs(yk),'r')
grid on
legend('原信号的频谱','滤波后信号的频谱')
