go2.m

这是一个优化IIR滤波器的源码,能设计出4个极点0个零点的IIR优化滤波器

%% Matlab script to run the optimization
%% Code Generated by Mathematica packages written by Brian L. Evans
%% Copyright 1993 by the Regents of the University of California

%% This program uses the following global variables:
global dp ds maxO wmax Wmp Wms Wmt WO wp Wph WQ ws tpeak1;

%% The following constants have been given default values:
Wmp = 0;             %% weight on mag. response in passband
Wmt = 0;             %% weight on mag. response in transition band
Wms = 1;             %% weight on mag. response in stopband
Wph = .1;           %% weight on deviation from linear phase
WQ = 2;              %% weight on quality factor cost
WO = 2;              %% weight on percent overshoot cost
wp = 20;             %% end of passband frequency
ws = 40;             %% start of stopband frequency
wmax = wp;           %% max. frequency for filter to exhibit linear phase
maxO = 20;           %% max. percent overshoot of the step response
dp = 0.21;           %% allowable mag. error in the passband
ds = 0.31;           %% allowable mag. error in the stopband

%% The program uses one global variable that is updated at each iteration:
%%	tpeak1 the time at which the peak in the step response occurs

%% Debugging code
global numtimes weights;
weights = [0 0 0 0];
numtimes = 0;

%% Initial Guess x0 is a Butterworth filter

%% We can choose the cutoff frequency to be between these two frequencies:
wc1 = wp./(-1. + (1. - 1..*dp).^(-2)).^(1./8);
wc2 = ws./(-1. + ds.^(-2)).^(1./8);
if wc1 > wc2
  error( 'The filter order is too low.' );
end;

%% Choose the middle frequency for wc
wc = (wc1 + wc2)./2;
x0 = wc*[cos((11.*pi)./8), cos((9.*pi)./8), sin((11.*pi)./8), sin((9.*pi)./8)];
tpeak1 = fmin( 'stepresp', 0, max(-(pi./x0(3)), -(pi./x0(4))), [], x0(1:2), x0(3:4),x0(5:4),x0(5:4), -1 );

%[-(pi./x0(3)) -(pi./x0(4))]

if ( stepresp(tpeak1, x0(1:2), x0(3:4), x0(5:4),x0(5:4) , 1) - 1 ) > maxO./100
  fprintf(1, '  overshoot for the initial guess is higher than minimum.\n');
end
%% Do 10 iterations at a time and quit when the first three %%   digits of the value of the cost function remain constant

options = foptions;
options(1) = 1;
options(9) = 1;
options(14) = 100;
errvec = [-100.*eps, -100.*eps, -100.*eps, -100.*eps];
[root, options] = constr('filcost2', x0, options, [], errvec, 'costder2');
%lastcost = filcost2(x0);
%while ( abs( (lastcost - options(8) ) / lastcost) > 1e-3 ),
 % lastcost = options(8);
 % [root, options] = constr('filcost2', root, options, [], errvec, 'costder2');
%end;
disp('The conjugate zero pairs are located at:');
zeroloc = root(5:4) +i*root(5:4);
disp(zeroloc);
disp('The conjugate pole pairs are located at:');
poleloc = root(1:2) +i*root(3:4);
disp(poleloc);

%% Normalized weight vector
weights / numtimes