go4.m

这是一个优化IIR滤波器的源码,能设计出8个极点8个零点的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
%% @(#)go4.m	1.5	10/30/01

function foo = go2(fileid,a,b,c,d,e,f,g,h,ii,j,k,l,wmin,fmax,tmin,tmax,npts)

%% 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 = a;             %% weight on mag. response in passband
Wmt = b;             %% weight on mag. response in transition band
Wms = c;             %% weight on mag. response in stopband
Wph = d;           %% weight on deviation from linear phase
WQ = e;              %% weight on quality factor cost
WO = f;              %% weight on percent overshoot cost
wp = g;             %% end of passband frequency
ws = h;             %% start of stopband frequency
wmax = wp;           %% max. frequency for filter to exhibit linear phase
maxO = j;           %% max. percent overshoot of the step response
dp = k;           %% allowable mag. error in the passband
ds = l;           %% 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;

x0=[-8107.2614 -5483.3416 -2751.9711 -802.3032 -3836.7800 -9646.0115 -12507.3953 -13554.1190  -9.0949e-13 -9.0949e-13 -9.0949e-13 -9.0949e-13 -17813.1598 -76471.8605 -28205.4248 -20117.5622];

tpeak1 = fmin( 'stepresp', 0, max([-(pi./x0(5)), -(pi./x0(6)), -(pi./x0(7)), -(pi./x0(8))]), [], x0(1:4), x0(5:8),x0(9:12),x0(13:16), -1 );
if ( stepresp(tpeak1, x0(1:4), x0(5:8), x0(9:12),x0(13:16) , 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) = 2;
options(9) = 1;
options(14) = 100;
errvec = [-100.*eps, -100.*eps, -100.*eps, -100.*eps, -100.*eps, -100.*eps, -100.*eps, -100.*eps, -100.*eps, -100.*eps, -100.*eps, -100.*eps, -100.*eps, -100.*eps, -100.*eps, -100.*eps];
[root, options] = constr('filcost4', x0, options, [], errvec, 'costder4');
lastcost = filcost4(x0);
while ( abs( (lastcost - options(8) ) / lastcost) > 1e-3 ),
  lastcost = options(8);
  [root, options] = constr('filcost4', root, options, [], errvec, 'costder4');
end;
disp('The conjugate zero pairs are located at:');
zeroloc = root(9:12) +i*root(13:16);
disp(zeroloc);
disp('The conjugate pole pairs are located at:');
poleloc = root(1:4) +i*root(5:8);
disp(poleloc);

%% Normalized weight vector
weights / numtimes
mag = magphase(fileid,poleloc(1:1),poleloc(2:2),poleloc(3:3),poleloc(4:4),zeroloc(1:1),zeroloc(2:2),zeroloc(3:3),zeroloc(4:4),fmax,wmin,npts)
step = stepresp([tmin :((tmax-tmin)/npts):tmax], [real(poleloc(1:1)) real(poleloc(2:2)) real(poleloc(3:3)) real(poleloc(4:4))],[imag(poleloc(1:1)) imag(poleloc(2:2)) imag(poleloc(3:3)) imag(poleloc(4:4))], [real(zeroloc(1:1)) real(zeroloc(2:2)) real(zeroloc(3:3)) real(zeroloc(4:4))],[imag(zeroloc(1:1)) imag(zeroloc(2:2)) imag(zeroloc(3:3)) imag(zeroloc(4:4))],1)
fid = fopen(fileid,'w');
for i = 1 : length(mag)
   fprintf(fid,'%f\t%f\t%f\n',mag(1,i),mag(2,i),step(i));
end
fprintf(fid,'%f\n%f\n%f\n%f\n%f\n%f\n%f\n%f\n%f\n%f\n%f\n%f\n%f\n%f\n%f\n%f',real(poleloc(1:1)),imag(poleloc(1:1)),real(poleloc(2:2)),imag(poleloc(2:2)),real(poleloc(3:3)),imag(poleloc(3:3)),real(poleloc(4:4)),imag(poleloc(4:4)),real(zeroloc(1:1)),imag(zeroloc(1:1)),real(zeroloc(2:2)),imag(zeroloc(2:2)),real(zeroloc(3:3)),imag(zeroloc(3:3)),real(zeroloc(4:4)), imag(zeroloc(4:4)))
fclose(fid);
foo = 4;

end