iir_lpf.m

this base on DSP a iir design

%writen by Ziwei Chen 
clc;clear all;close all;
%产生模拟信号x(t)
t=[0:0.00001:1];
f1=100;f2=1400;f3=3000;
xt=sin(2*pi*f1*t)+sin(2*pi*f2*t)+sin(2*pi*f3*t);
h1=figure;
set(h1,'name','进行滤波的输入信号');
subplot(2,1,1);
plot(t,xt);
xlabel('时间t （单位：s）');ylabel('幅度x(t)');
title('模拟时间信号x(t)');

%计算采样周期Ts及采样点数N,对原时间信号进行截断，选取Tp时间进行观察
fs=10000;%fs=80000;
Ts=1/fs;
%F=10;
N=128;%fs/F;                               %N=input('请输入采样点数N＝');   
Tp=N*Ts;                              %最小记录时间（对信号的观察时间）

%对模拟信号进行采样，得到序列信号x(n)
n=[0:N-1];
xn=sin(2*pi*f1*n*Ts)+sin(2*pi*f2*n*Ts)+sin(2*pi*f3*n*Ts);
subplot(2,1,2);
stem(n,xn);axis([0 1024 -3 3]);
xlabel('时间n');ylabel('幅度x(n)');
title('离散时间信号x(n)')
name=fopen('inputsignal.dat','w');
count=fwrite(name,xn,'short')
fid=fclose(name)
%对序列信号x(n)进行谱分析，得到信号的频谱
K=[0:N-1];
Nfft=N;
XK=fft(xn,Nfft);
h2=figure;
set(h2,'name','对输入信号进行频谱分析的结果（低通滤波前）');
plot(fs*K/Nfft,abs(XK));grid on;
set(gca,'XTickMode','manual','XTick',[100,1400,3000,4000,5000,6600,7900]);
xlabel('频率f  (单位：Hz)');ylabel('幅度|X(j2\pif)|');
title('输入信号的幅度谱（离散频谱）');
%----------------------------------------------------------------

%数字低通滤波器设计
wp = 2*pi*1000/fs;%fpass=1000Hz %wp =0.25*pi
ws = 2*pi*1500/fs;%fstop=1500Hz deltaf=500Hz %ws = 0.5*pi 
Ap =4; %Ap =0.2;
As =15;%As =40; 
T = 2;
wp = (2/T)*tan(wp/2);
ws = (2/T)*tan(ws/2);
[N,Wc] = buttord(wp,ws,Ap,As,'s');
[z,p,k] = buttap(N);
[b,a] = zp2tf(z,p,k);
[bt,at] = lp2lp(b,a,Wc);
h3=figure;
set(h3,'name','模拟低通巴特沃斯滤波器');
[H,W] = freqs(bt,at);
plot(W*fs/(2*pi),abs(H));grid on;
%set(gca,'XTickMode','manual','XTick',[100,1400,3000,4000,8000]);
set(gca,'XTickMode','manual','XTick',[100,1400,2000,2500,3000,4000]);
title('模拟低通巴特沃斯滤波器')

[bd,ad] = bilinear(bt,at,1/T)

tf2sos(bd,ad);
[H,F] = FREQZ(bd,ad,Nfft,'whole',fs);%[H,W] = freqz(bd,ad);
mag = abs(H);
db = 20*log10((mag+eps)/max(mag));
h4=figure;
set(h4,'name','数字低通巴特沃斯滤波器');
plot(F,db);grid on;%plot(W/(2*pi)*fs,db);
set(gca,'XTickMode','manual','XTick',[100,1400,2000,2500,3000,4000,8000]);
title('数字低通巴特沃斯滤波器')
%----------------------------------------------------------------

%对输入序列进行低通滤波
YK=XK.*H';
%y=xcorr(xn,ifft(H,Nfft));
%YK=fft(y,Nfft);
h5=figure;
set(h5,'name','对输入信号进行低通滤波的结果（低通滤波后）');
plot(fs*K/Nfft,abs(YK));grid on;%axis tight;
set(gca,'XTickMode','manual','XTick',[100,1400,3000,4000,5000,6600,7900]);
xlabel('频率f  (单位：Hz)');ylabel('幅度|Y(j2\pif)|');
title('输出信号的幅度谱（离散频谱）');