low_mufilter.m

信号处理－滤波器设计（基于matlab和Mathmatica的设计方法）中关于无乘法器椭圆IIR滤波器设计的原创程序。

function [num_result,den_result]=Low_Mufilter(Fp,Fs,Ap,As)
%  Low_Mufilter.m  设计无乘法器滤波器
%  $date 4/10/2005
%
%   作者：赵泽平
%      [e-mail]
%   
%
%   References:
%        Miroslav D. Lutovac, Dejan V. Tosic, Brian L. Evans
%           Advanced Filter Design for Signal Processing
%                   Using MATLAB and Mathematica

%% calls: 
%        Low_Mu_fp.m,Low_Mu_e.m,Low_Mu_Bmax.m,Low_Mu_L.m,Low_Mu_sort.m,
%        Low_Mu_dX.m,Low_Mul_D.m,Low_Mul_Pha_plot.m,Low_Mul_quan

global n f3dB Bmax 
%求n
if nargin ~=4
   error('Requires 4 inputs.')
end
n=ellipord(2*Fp,2*Fs,Ap,As);
if abs(n-2*fix(n/2)) == 0
    n=n+1;
end
disp(sprintf('最小阶数为n=%d',n))
n=input('滤波器的阶数（必须为奇数）');
if abs(n-2*fix(n/2)) ~=1
   error('滤波器的阶数必须是奇数')
end

%求f3dB
amin=((tan(pi*Fp))^2-1)/(1+(tan(pi*Fp))^2);
amax=((tan(pi*Fs))^2-1)/(1+(tan(pi*Fs))^2);
a=(tan(pi*Fp)*tan(pi*Fs)-1)/(1+tan(pi*Fp)*tan(pi*Fs));
disp(sprintf('a值的范围：%f<α<%f',amin,amax))
disp(sprintf('α＝%f(参考值)',a))
a=input('选择： α=');
f3dB=atan(((1+a)/(1-a))^0.5)/pi;
%求Bmax1
fp1=Low_Mu_fp(f3dB,Fs);
e1=Low_Mu_e(fp1,Fs);
L=Low_Mu_eL(n,e1);
ap1=10*log10(1+1/L);
if ap1>Ap
    error(sprintf('通带衰减的要求太严格，它必须大于%f',ap1));
end
as1=10*log10(1+L);
if as1<As
    error(sprintf('阻带衰减的要求太严格，它必须小于%f',as1));
end
[z,p,k]=ellip(n,ap1,as1,2*fp1);
Bmax1=max(abs(p).^2);
%求Bmax2
L=10^(As/10)-1;
e2=Low_Mu_Le(n,L);
fp2=atan(sqrt((tan(pi*f3dB))^2/e2))/pi;
ap2=10*log10(1+1/L);
as2=As;
[z,p,k]=ellip(n,ap2,as2,2*fp2);
Bmax2=max(abs(p).^2);
disp(sprintf('βmax值的范围：%f<βmax<%f',Bmax1,Bmax2))
Bmax=input('选择：βmax=');
%求Bi和a1
e=fzero(@Low_Mu_Bmax,(e1+e2)/2);
fp=atan(sqrt((tan(pi*f3dB)^2/e)))/pi;
L=Low_Mu_eL(n,e);
ap=10*log10(1+1/L);
as=10*log10(1+L);
[numD,denD]=ellip(n,ap,as,2*fp);
zplane(numD,denD);       %绘制零极点
[z,p,c]=ellip(n,ap,as,2*fp);
a1=-(1-tan(pi*f3dB))/(1+tan(pi*f3dB));
%对Bi按Ha和Hb进行排序和绘制d/dX图
[Bi,serial_num]=Low_Mu_sort(n,p);
Bi(serial_num)=Bmax; %对Bmax进行替换
Bi(1)=a1;  %把a1归入Bi数组，方便程序编写
Low_Mu_dX(n,a,Bi,serial_num,Fs);

len=0;
while len~=(n-1)/2
    disp('输入量化的顺序，其中Seq=1表示α1,Seq=2,3,...表示β2,β3,...')
    Seq=input('  Seq=');
    len=length(Seq);
end
figure;
[Dp,Ds]=Low_Mul_D(Ap,As);
[Pha_Sub,f]=Low_Mul_Pha_plot(n,a,Bi,Fs,Ds);
Bq=Bi;
for i=1:len
   B=Low_Mul_quan(a,Bq,Fs,Pha_Sub,Ds,f,Seq(i));
   Bq(Seq(i))=B;
   [Pha_Sub,f]=Low_Mul_Pha_plot(n,a,Bq,Fs,Ds);
end
%计算幅频相频响应
[num_result,den_result]=low_Mul_Ma_Ph(n,a,Bi);