对周期方波信号进行滤波matlab的实现

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('原信号的频谱','滤波后信号的频谱')

%%%%%%%%%%%%%%%%%%%%保留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('原信号的频谱','滤波后信号的频谱')