求脉冲响应不变法设计Elliptic型IIR数字低通滤波器
发布网友
发布时间：2022-04-20 08:14
共2个回答
热心网友
时间：2023-10-04 05:35
代码为：
fp = 30;
fs = 35;
Fs = 100;
wp = 2*pi*fp/Fs;
ws = 2*pi*fs/Fs;
wp = tan(wp/2);
ws = tan(ws/2); % 通带最大衰减为0.5dB，阻带最小衰减为40dB
[N, wn] = buttord(wp, ws, 0.5, 40, 's'); % 模拟低通滤波器极零点
[z, p, k] = buttap(N); % 由极零点获得转移函数参数
[b, a] = zp2tf(z, p, k); % 由原型滤波器获得实际低通滤波器
[B, A] = lp2lp(b, a, wp);
[bz, az] = bilinear(B, A, .5);
还有三句。。。给分后给你
热心网友
时间：2023-10-04 05:35
wp=0.4*pi;ws=0.5*pi;
rp=1;as=30;
Fs=20000;T=1/Fs;
OmgP=wp/T;
OmgS=ws/T;
ep=sqrt(10^(rp/10)-1);
ripple=sqrt(1/(1+ep*ep));
attn=1/(10^(as/20));
[cs,ds]=afd_elip(OmgP,OmgS,rp,as);
[b,a]=impinvar(cs,ds,Fs);
[db,mag,pha,grd,w]=freqz_m(b,a);
%
%绘出各条曲线
subplot(2,2,1);plot(w/pi,mag/max(mag));title('幅频特性');
xlabel('w(/pi)');ylabel('|H(jw)|');
axis([0,1,0,1.1]);
set(gca,'XTickMode','manual','XTick',[0 0.4 0.5 1]);
set(gca,'YTickMode','manual','YTick',[0 attn ripple 1]);grid
subplot(2,2,2);plot(w/pi,db);title('幅频特性(db)');
xlabel('w(/pi)');ylabel('dB');
axis([0,1,-50,5]);
set(gca,'XTickMode','manual','XTick',[0 0.4 0.5 1]);
set(gca,'YTickMode','manual','YTick',[-60 -as -rp 0]);grid
subplot(2,2,3);plot(w/pi,pha/pi);title('相频特性');
xlabel('w(/pi)');ylabel('pha(/pi)');
axis([0,1,-1,1]);
set(gca,'XTickMode','manual','XTick',[0 0.4 0.5 1]);grid
subplot(2,2,4);plot(w/pi,grd);title('群延时');
xlabel('w(/pi)');ylabel('Sample');
axis([0,1,0,20]);
set(gca,'XTickMode','manual','XTick',[0 0.4 0.5 1]);grid
其中用到freqz_m函数和afd_elip函数，你的matlab里可能没有，我把代码给出，你把下面的代码复制进txt文本里，分别重命名为freqz_m.m和afd_elip.m并放到work文件夹里。再运行上面的主程序。
函数一
function [db,mag,pha,grd,w] = freqz_m(b,a);
% Modified version of freqz subroutine
% ------------------------------------
% [db,mag,pha,grd,w] = freqz_m(b,a);
% db = Relative magnitude in dB computed over 0 to pi radians
% mag = absolute magnitude computed over 0 to pi radians
% pha = Phase response in radians over 0 to pi radians
% grd = Group delay over 0 to pi radians
% w = 501 frequency samples between 0 to pi radians
% b = numerator polynomial of H(z) (for FIR: b=h)
% a = denominator polynomial of H(z) (for FIR: a=[1])
%
[H,w] = freqz(b,a,1000,'whole');
H = (H(1:1:501))'; w = (w(1:1:501))';
mag = abs(H);
db = 20*log10((mag+eps)/max(mag));
pha = angle(H);
% pha = unwrap(angle(H));
grd = grpdelay(b,a,w);
% grd = diff(pha);
% grd = [grd(1) grd];
% grd = [0 grd(1:1:500); grd; grd(2:1:501) 0];
% grd = median(grd)*500/pi;
函数二
function [b,a] = afd_elip(Wp,Ws,Rp,As);
% Analog Lowpass Filter Design: Elliptic
% --------------------------------------
% [b,a] = afd_elip(Wp,Ws,Rp,As);
% b = Numerator coefficients of Ha(s)
% a = Denominator coefficients of Ha(s)
% Wp = Passband edge frequency in rad/sec; Wp > 0
% Ws = Stopband edge frequency in rad/sec; Ws > Wp > 0
% Rp = Passband ripple in +dB; (Rp > 0)
% As = Stopband attenuation in +dB; (As > 0)
%
if Wp <= 0
error('Passband edge must be larger than 0')
end
if Ws <= Wp
error('Stopband edge must be larger than Passband edge')
end
if (Rp <= 0) | (As < 0)
error('PB ripple and/or SB attenuation ust be larger than 0')
end
ep = sqrt(10^(Rp/10)-1);
A = 10^(As/20);
OmegaC = Wp;
k = Wp/Ws;
k1 = ep/sqrt(A*A-1);
capk = ellipke([k.^2 1-k.^2]); % Version 4.0 code
capk1 = ellipke([(k1 .^2) 1-(k1 .^2)]); % Version 4.0 code
N = ceil(capk(1)*capk1(2)/(capk(2)*capk1(1)));
fprintf('\n*** Elliptic Filter Order = %2.0f \n',N)
[b,a]=u_elipap(N,Rp,As,OmegaC);
