%% 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 passbandWmt = 0;             %% weight on mag. response in transition bandWms = 0;             %% weight on mag. response in stopbandWph = 0.1;           %% weight on deviation from linear phaseWQ = 0;              %% weight on quality factor costWO = 1;              %% weight on percent overshoot costwp = 2*pi*2000;             %% end of passband frequencyws = 2*pi*4000;             %% start of stopband frequencywmax = wp;           %% max. frequency for filter to exhibit linear phasemaxO = 20;           %% max. percent overshoot of the step responsedp = 5.0119e-04;           %% allowable mag. error in the passbandds = 5.0119e-04;           %% 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 codeglobal 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 constantoptions = 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 vectorweights / numtimes