% Construct a first-derivative estimator filter using LLS methods. % Copyright (c) Larry Trammell, 2017. % This Octave (Matlab) script file is distributed under the % Creative Commons Attribution-ShareAlike 4.0 International % License. In general terms, you can use this freely, with some % restrictions on adapting and redistributing modified versions. % For full details, see: % http://creativecommons.org/licenses/by-sa/4.0/ % % ================================================================== % Adjustable design parameters: % == Bandwidth == nominal end of fit interval, integer %-of-Nyquist BW = 14; % % == Stopband == start of attenuation band, integer %-of-Nyquist SB = 40; % % == Filter size == number of terms in final filter, odd integer Nfilt = 15; % % == Low band fit == weighting factor, positive lweight = 1.0; % % == Stopband rejection == weighting factor, positive sweight = 0.1; % % == Low band flatness == weighting factor, zero or positive flatness = 20.0; % ================================================================== % Derived parameters, do not alter N = Nfilt/2; Mfit = 100; fstep = pi/Mfit; % Reserve constraint and target value matrices mmat = zeros(Mfit,N); rmat = zeros(Mfit,1); % Add constraints for each frequency step for jterm=1:Mfit % Calculate weighting to apply at each discrete frequency weight = 0; if jterm<=BW weight = lweight + flatness/jterm; else if jterm>=SB weight = sweight; end end % Calculate weighted independent values from waveforms for iterm=1:N freq = fstep*iterm; mmat(jterm,iterm) = sin(fstep*iterm*(jterm)) * weight; end % Calculate weighted dependent values from "ideal response" if jterm<=BW rmat(jterm) = fstep*jterm * weight; end end % Solve for filter parameters parms = mmat\rmat; % Reorganize terms in anti-symmetric filter form filt = zeros(1,Nfilt); for term=1:N filt(Nfit+1+term) = parms(term)*(0.5); filt(Nfit+1-term) = parms(term)*(-0.5); end output_precision(7); filt % Reorganize terms in anti-symmetric filter form filt = zeros(1,Nfilt); for term=1:N filt(term) = parms(term)*(-0.5); filt(Nfilt+1-term) = -filt(term); end output_precision(7); filt % Frequency response on linear scale fr = zeros(1,Mfit); ref = zeros(1,Mfit); for jterm=1:Mfit resp = 0.0; for iterm=1:N freq = jterm*fstep; sinval = sin(iterm*freq); resp = resp + sinval*parms(iterm); end if jterm<=20 ref(jterm)=freq; end fr(jterm) = resp; end % Plot frequency response of fit figure(1) plot(1:Mfit, ref(1:Mfit), 'g', 1:Mfit, fr(1:Mfit), 'b') title('Least squares and ideal estimator frequency response') legend('ideal','estimator')