用双线性变换法设计BUTTERWORTH低通IIR数字滤波器

用双线性变换法设计BUTTERWORTH低通IIR数字滤波器：wp=0.2*2*pi;%通带边界频率ws=0.3*2*pi;%阻带截止频率rp=1;%通带最大衰减rs=25;%阻带最小衰减fs=1;%采样频率ts=1/fs;%采样周期wp2=2*fs*tan(wp/2*ts);%预畸变校正ws2=2*fs*tan(ws/2*ts);%预畸变校正[n,wn]=buttord(wp2,ws2,rp,rs,'s')%带入经预畸变后获得的归一化模拟频率参数[z,p,k]=buttap(n)%创建归一化的BUTTERWORTH模拟低通原型滤波器[bap,aap]=zp2tf(z,p,k)%把滤波器零极点模型转化为传递函数模型[b,a]=lp2lp(bap,aap,wn)%把模拟滤波器原型转换成截止频率为Wn的模拟低通滤波器[bz,az]=bilinear(b,a,fs);%用双线性变换法实现模拟滤波器到数字滤波器的转换[h,w]=freqz(bz,az);%绘制频率响应曲线subplot(2,1,1);plot(w/pi,abs(h));gridonxlabel('频率');ylabel('幅度');subplot(2,1,2);plot(w/pi,20*log10(abs(h)));gridonxlabel('频率');ylabel('幅度');用窗口法设计一个线性相位FIR低通滤波器，采用汉宁窗设计：wp=0.2*pi;%通带边界频率ws=0.3*pi;%阻带截止频率wdelta=ws-wp;%过渡带宽N=ceil(8*pi/wdelta);%求出N值wc=(0.2+0.3)*pi/2;%求出截止频率r=(N-1)/2;n=0:N-1;hdn=sin(wc*(n-r))./[pi*(n-r)];%理想的单位冲激响应wn=hanning(N);%求出汉宁窗的表达式h=hdn.*wn';%滤波器加窗H=fft(h,512);w=2*[0:511]/512plot(w,20*log10(abs(H)));用MATLAB实现时间抽选的基2－FFT算法:functiony=myditfft(x)%本程序对输入序列x实现时间抽选的基2－FFT，%点数取大于等于x长度的2的幂次m=nextpow2(x);N=2^m;%求x的长度对应的2的最低幂次miflength(x)Nx=[x,zeros(1,N-length(x))];endnxd=bin2dec(fliplr(dec2bin([1:N]-1,m)))+1;%1：2^m数列的倒位序y=x(nxd);%将x倒位序排列作为y的初始值formm=1:m%将DFT作m次基2分解，从左到右Nmr=2^mm;u=1;%旋转因子u初始化wN^0=1WN=exp(-i*2*pi/Nmr);%当前次分解的基本DFT因子wN=exp(-i*2*pi/Nmr)forj=1:Nmr/2%当前次跨越间隔内的各次蝶形运算fork=j:Nmr:N%当前次蝶形运算的跨越间隔为Nmr＝2^mmkp=k+Nmr/2;%确定蝶形运算的对应单元下标t=y(kp)*u;%蝶形运算的乘积项y(kp)=y(k)-t;%蝶形运算的减法项y(k)=y(k)+t;%蝶形运算的加法项endu=u*WN;%修改旋转因子，多乘一个基本DFT因子wNendend