% This script constructs a second-derivative filter design using % mini-max methods. It requires the glpk function provided by the % Octave package (a wrapper for the GNU GLPK). % ENTER THE DESIGN PARAMETERS HERE: % Passband width (typicaly 0.05 to 0.15) % Transition width (typically about twice the passband) % Number of unique filter terms to compute (besides frequency 0) % Error sensitivity in passband (typically 50 to 500) passb = 0.08 transition = 0.14 nfreqs = 7 psens = 150.0 % Fixed implementation parameters ncycle = 800; zerosens = 1e-5; % Derived internal parameters nsteps = ncycle/2 ; delfreq = 2*pi/ncycle; % Accounting for matrix rows used to span passband ipstart = 1 ipend = ceil(passb*ncycle+1) nfit = ipend-ipstart+1 % Accounting for matrix rows used to span stopband isstart = nfit + 1 + floor(transition*ncycle) nstop = nsteps-isstart+1 nterms = nfit+nstop % Matrix form: | -B B m | * V <= | -b | % | B -B m | | b | % % Cost: C = | 0 0 1 | % Vector: V = | +x -x tol |' % Rows in 'tableau' matrix matrows = 2*(nterms); % Columns in 'tableau' matrix, frequencies plus 'tolerance' column matcols = 2*(nfreqs+1)+1; % Construct the desired frequency response level for a % derivative filter. Parabolic values in fit band, 0 elsewhere. % Enforced 0 at frequency 0. bterms = zeros(nterms,1); for iterm=1:nfit bterms(iterm) = -(delfreq*(iterm-1))^2; end % Construct basis terms by sampling cosine functions. % Omit transition band terms. cosbasis = zeros(nterms,nfreqs+1); for kterm=1:nfreqs+1 for iterm=1:nfit cosbasis(iterm,kterm) = cos((kterm-1)*(iterm-1)*delfreq); end for iterm=1:nstop cosbasis(nfit+iterm,kterm) = cos((kterm-1)*(isstart+iterm-2)*delfreq); end end % Construct passband/stopband weights for tolerance column mcol = [ones(nfit,1); ones(nstop,1)*psens]; mcol(1) = zerosens; %mcol % We now have enough information to construct the "tableau" matrix tableau = zeros(matrows,matcols); % -- sine function basis columns, two copies of sine basis tableau(1:nterms,1:nfreqs+1) = -cosbasis; tableau(nterms+1:2*nterms,1:nfreqs+1) = cosbasis; tableau(1:nterms,nfreqs+2:2*nfreqs+2) = cosbasis; tableau(nterms+1:2*nterms,nfreqs+2:2*nfreqs+2) = -cosbasis; % -- tolerance bound column tableau(1:nterms,matcols) = -mcol; tableau(nterms+1:2*nterms,matcols) = -mcol; %tableau % -- target values vector vb = [ -bterms; bterms ]; %vb % -- cost vector vcost = zeros(1,matcols); vcost(matcols) = 1.0; %disp('Initial tableau:') %tableau % Required for sign-restricted form of glpk solver lbmat = zeros(matcols,1); ubmat = ones(matcols,1)*1e12; clear ctype; for kterm=1:matrows ctype(kterm)="U"; end % Calculate the min-max optimal fit over the available data set [xopt, fmin, status, extra] = glpk(vcost, tableau, vb, lbmat,ubmat,ctype); disp('solution status (look for 180):'),disp(status); %disp('xopt solution terms: '),disp(xopt'); disp('stopband ripple: '),disp(xopt(matcols)*psens); disp('passband error: '),disp(xopt(matcols)); % Reorganize solution terms as a filter cfvect = zeros(1,2*nfreqs+1); center = nfreqs+1; cfvect(center) = (xopt(1)-xopt(nfreqs+2)); for kterm=1:nfreqs termval = (xopt(kterm+1)-xopt(nfreqs+2+kterm)) * 0.5; cfvect(center+kterm) = termval; cfvect(center-kterm) = termval; end disp('DC error: '),disp(sum(cfvect)); [mags,phase] = freqz (cfvect, [1], nsteps+1); figure(1) plot([0:nsteps-1]*(0.5/nsteps), abs(mags(1:nsteps)),'b') title('Filter frequency response characteristic'); axis([0.0, 0.5, 0.0, 1.0]) disp('filter: '),disp(cfvect);