信号滤波及时频特性分析

clc;clear all;
dt=0.01;fs=1/dt;
t = 0:dt:20+dt; %时间序列
a = 0.1*pi;  %周期参数
Y=square(a*t,10); %生成方波,T=2pi/a,a=2*pi*f=2*pi/T,a为角频率,第二参数为占空比，单位为%，如20则为20%
figure(1);plot(t,Y,'r','linewidth',2);
N=length(t);
f=(-N/2:N/2-1)*(fs/N);%自变量频率范围
F1=fft(Y);%快速傅里叶变换
F1=fftshift(F1);%调制频谱图，将零频分量移到频谱中心
F1=abs(F1);
figure(2);plot(f,F1);xlim([0 fs/2]);

wp=0.5;%通带边界频率，模拟频率
ws=3;%阻带边界频率
rp=1;%通带最大衰减dB
rs=30;%阻带最大衰减
[nn,w]=buttord(wp/(fs/2),ws/(fs/2),rp,rs)%确定阶数n和3dB截止频率
[b,a]=butter(nn,w);%确定H（S）的分子分母系数向量
[H,W]=freqz(b,a);%H即频率响应函数，W即数字角频率
figure(3);
plot(W*fs/(2*pi),abs(H));%滤波器幅频特性，自变量为模拟频率
xlabel('f/Hz');ylabel('放大倍数(倍)');title('幅频特性');

Y1=filter(b,a,Y);%低通滤波后的信号
figure(4);plot(t,Y1);

F2=fft(Y1);
F2=fftshift(F2);%调制频谱图，将零频分量移到频谱中心
F2=abs(F2);
figure(5);plot(f,F2);%滤波后的方波信号的FFT
xlim([0 fs/2]);