Good evening (here it is, at least)! I'm having a little trouble with the design of a linear phase FIR. I'm using remez (minimax/equiripple) in matlab. My specification is; +-A dB in the passband and -B dB in the stopband, band edges at \omega_{p} and \omega_{s}. Now I need to choose my weighting so that I get the right amount of ripple in each band, but I just can't get it right!! Here is what I do: I'm using an equiripple method so the ripple in the pass- & stopbands must be equal in magnitude when the algorithm has done its thing : max|passband error| = max|stopband error| (|...|=abs(...)) Now the weights come in: 1*max|passband error| = vs*max|stopband error| taking 20*log10(...) on both sides, and inserting values from the spec: -A dB = 20*log10(vs) - B dB finally: vs = 10^((A-B)/20) This doesn't work as intended, remez(N,[0 wp ws 1][1 1 0 0][1 vs]) does not give me the filter I wanted. What am I doing wrong? Anyone? This message was sent using the Comp.DSP web interface on www.DSPRelated.com
minimax and choosing weights
Started by ●September 16, 2005
Reply by ●September 16, 20052005-09-16
"MA" <oooncs@yahoo.se> wrote in message news:Yc-dneEOW71_s7beRVn-sQ@giganews.com...> Good evening (here it is, at least)! I'm having a little trouble with the > design of a linear phase FIR. I'm using remez (minimax/equiripple) in > matlab. My specification is; +-A dB in the passband and -B dB in the > stopband, band edges at \omega_{p} and \omega_{s}. Now I need to choose my > weighting so that I get the right amount of ripple in each band, but I > just > can't get it right!! Here is what I do: > > I'm using an equiripple method so the ripple in the pass- & stopbands must > be equal in magnitude when the algorithm has done its thing : > > max|passband error| = max|stopband error| > > (|...|=abs(...)) Now the weights come in: > > 1*max|passband error| = vs*max|stopband error| > > taking 20*log10(...) on both sides, and inserting values from the spec: > > -A dB = 20*log10(vs) - B dB > > finally: vs = 10^((A-B)/20) > > > This doesn't work as intended, remez(N,[0 wp ws 1][1 1 0 0][1 vs]) does > not give me the filter I wanted. What am I doing wrong? Anyone? >Most probably the weights apply to the absolute error and not the error in dB. Example: Passband ripple .5 dB (assuming passband gain is specified to be 1.0) .5dB is 1.00925 absolute or 0.00925 peak error. Stopband ripple -40dB which is 0.01 peak error. In this case, the passband error is least by 0.01/0.00925 = 1.01108 So, weigh the passband by 1.01108 and the stopband by 1.0. I hope this helps. Fred
Reply by ●September 16, 20052005-09-16
Thanks for replying Fred! I still don't get it Did you really get the numbers right there? 20*log10(1.00925) = 0.08dB If we're looking at a +-1dB oscillation around 1, we have A_minimum=0.8913 and A_maximum=1.1220. The midpoint in between is at (A_maximum+A_minimum)/2=1.0066. The maximum pass-band error must be A_maximum - 1 then? I'm getting more and more confused... Those dB's... This message was sent using the Comp.DSP web interface on www.DSPRelated.com
Reply by ●September 16, 20052005-09-16
You only really have control over the ratio of the ripples (peak-peak, not in dB) using the weights. This is controlled by the ratio of the weights, so if you increase all of the weights by a single factor the filter shouldn't change. Therefor if you set one weight to a fixed number, and vary the other, keeping the other parameters fixed, you will see all the possible responses you can get. But ... keep in mind that you may not be able to meet the spec simply by varying the weights.... like if N is the wrong size (strongly tied to wp and ws) you may never meet spec. dirk
Reply by ●September 16, 20052005-09-16
MA wrote:> Thanks for replying Fred! I still don't get it > > Did you really get the numbers right there? 20*log10(1.00925) = 0.08dB > > If we're looking at a +-1dB oscillation around 1, we have A_minimum=0.8913 > and A_maximum=1.1220. The midpoint in between is at > (A_maximum+A_minimum)/2=1.0066. The maximum pass-band error must be > A_maximum - 1 then? > > I'm getting more and more confused... Those dB's...You will probably benefit from reviewing the basic definitions you need to deal with. Ripple:(cyclic) deviation from nominal. Often measured as a fraction, in percent, or in decibels. Decibel: a ratio of powers: 10*log(P/P_std). Because power in a constant impedance varies with the square of voltage or current, that's equivalent to 10*log(v^2/P_std^2) = 20*log(V/V_std). I hope that clears up some of your confusion. Jerry -- Engineering is the art of making what you want from things you can get. �����������������������������������������������������������������������
Reply by ●September 16, 20052005-09-16
"MA" <oooncs@yahoo.se> wrote in message news:Ofidnf1xgOVyqrbeRVn-vA@giganews.com...> Thanks for replying Fred! I still don't get it > > Did you really get the numbers right there? 20*log10(1.00925) = 0.08dB > > If we're looking at a +-1dB oscillation around 1, we have A_minimum=0.8913 > and A_maximum=1.1220. The midpoint in between is at > (A_maximum+A_minimum)/2=1.0066. The maximum pass-band error must be > A_maximum - 1 then? > > I'm getting more and more confused... Those dB's... > > This message was sent using the Comp.DSP web interface on > www.DSPRelated.comOK.... The ripple in the passband is the *deviation* from the desired gain of 1.0000 The ripple in the stopband is the *deviation* from the desired gain of 0.00000 However, zero is not a good reference for dB (which is always a ratio of two gains) so stopband ripple is typically measured with respect to the desired passband gain (usually 1.0) and, The Remez algorithm internally deals with absolute values (weighted) and not dB. Only the results are stated in dB. So the example I gave said a passband ripple of 0.5 dB. If the measure is made absolute instead of in dB, the peak is 10^(0.5/20) = 1.00925 which means an absolute deviation of 0.00925. To get back to dB, you have to take the ratio of 1.00925:1 20*log10(1.00925) = 0.5dB I guess it's the passband that causes the confusion. A variation of 0.00925 relative to 1.000 is -40dB. A gain of 1.0925 relative to 1.000 is 0.5dB and it's this one we normally use for passband specifications. In the stopband the ripple given -40dB. 10^(-40/20) = 0.01 the absolute value of the peaks in the stopband. Using the ratio 0.01 : 1.000 gives: 20*log10(.010) = -40dB Fred
Reply by ●September 16, 20052005-09-16
If you're calculating the deviations incorrectly then how are you calculating the order of the filter. I believe this calculation needs the ripples expressed in terms deviation (absolute value). So you have the filter length incorrect - resulting in a filter that doesn't meet the spec. The following bit of code shows how to convert the ripples in dB to absolute scale. The weight vector is simply the ratio of the deviations. rp = 3; % Passband ripple rs = 40; % Stopband ripple fs = 2000; % Sampling frequency f = [500 600]; % Cutoff frequencies a = [1 0]; % Desired amplitudes % Compute deviations dev = [(10^(rp/20)-1)/(10^(rp/20)+1) 10^(-rs/20)]; [n,fo,ao,w] = remezord(f,a,dev,fs); b = remez(n,fo,ao,w); Hope that helps. Cheers, David MA wrote:> Good evening (here it is, at least)! I'm having a little trouble with the > design of a linear phase FIR. I'm using remez (minimax/equiripple) in > matlab. My specification is; +-A dB in the passband and -B dB in the > stopband, band edges at \omega_{p} and \omega_{s}. Now I need to choose my > weighting so that I get the right amount of ripple in each band, but I just > can't get it right!! Here is what I do: > > I'm using an equiripple method so the ripple in the pass- & stopbands must > be equal in magnitude when the algorithm has done its thing : > > max|passband error| = max|stopband error| > > (|...|=abs(...)) Now the weights come in: > > 1*max|passband error| = vs*max|stopband error| > > taking 20*log10(...) on both sides, and inserting values from the spec: > > -A dB = 20*log10(vs) - B dB > > finally: vs = 10^((A-B)/20) > > > This doesn't work as intended, remez(N,[0 wp ws 1][1 1 0 0][1 vs]) does > not give me the filter I wanted. What am I doing wrong? Anyone? > > > This message was sent using the Comp.DSP web interface on > www.DSPRelated.com
Reply by ●September 17, 20052005-09-17
Hi all! Thanks for helping me! Gotta love this forum and its users:) I think I got it right now... First: passband ripple between 1 and -1 dB magnitude means in linear scale magnitude between 1.1220 and 0.8913. The average magnitude in the passband is then 1.0066. It is important to distinguish this case with passband ripple given in dB from the case with symmetric deviation from unit magnitude. Consequently we have to adjust the "desired magnitude" to match the passband average: passbandAverage = (passbandmax+passbandmin)/2 passbandError = (passbandmax-passbandmin)/2 Next, the amplitude corresponding to 40dB below unity lies slightly more than 40 dB below the new value of the "desired magnitude" in the passband, so we have to compute a new deviation from 0 in the stopband: stopbandmax_dB2 = 20*log10(stopbandmax/passbandAverage) stopbandmax2 = 10^(stopbandmax_dB2/20) We can now send this info to our friend remezord: remezord([5 6],[passbandAverage 0],[ passbandError stopbandmax2 ],20) Or we can choose the weights ourselves: stopbandweight = passbandError/stopbandmax weight vector = [1 stopbandweight] Here below comes the code I used to sort out this, you can copy-paste it to see a nice plot: ------------------------------------------------------------------------ clear all % ------- spec -------- passbandmax_dB = 4; passbandmin_dB = -4; stopbandmax_dB = -40; % --------------------- passbandmax = 10^(passbandmax_dB/20) passbandmin = 10^(passbandmin_dB/20) stopbandmax = 10^(stopbandmax_dB/20) passbandPosDev = passbandmax-1 passbandNegDev = 1-passbandmin passbandMaxDev = max(passbandPosDev,passbandNegDev) passbandAverage = (passbandmax+passbandmin)/2 passbandError = (passbandmax-passbandmin)/2 stopbandmax_dB2 = 20*log10(stopbandmax/passbandAverage) stopbandmax2 = 10^(stopbandmax_dB2/20) % deviation = [ passbandMaxDev stopbandmax ]; deviation2 = [ passbandError stopbandmax2 ]; [N,F,A,W] = remezord([5 6],[1 0],deviation,20); [N2,F2,A2,W2] = remezord([5 6],[passbandAverage 0],deviation2,20) A2 = [passbandAverage passbandAverage 0 0] W2 = [1 passbandError/stopbandmax] b = remez(N,F,A,W); b2 = remez(N2,F2,A2,W2); H = freqz(b,1,2048); H2 = freqz(b2,1,2048); figure plot([0:2047]/2048,max(20*log10(abs(H)),-70),'k') hold on plot([0:2047]/2048,max(20*log10(abs(H2)),-70),'k:') hold on plot([0:2047]/2048,20*log10(passbandmax*ones(1,2048)),'k') hold on plot([0:2047]/2048,20*log10(passbandmin*ones(1,2048)),'k') hold on plot( [0:2047]/2048,20*log10(stopbandmax*ones(1,2048)),'k') title('Magnitude |H(e^{j \omega})| and |H_{2}(e^{j \omega})|') xlabel('Normalized frequency [units of \pi]') ylabel('Magnitude [dB]') legend('|H(e^{j \omega})|','|H_{2}(e^{j \omega})|') This message was sent using the Comp.DSP web interface on www.DSPRelated.com
Reply by ●September 17, 20052005-09-17
Hmmm... remezord seems to make an optimistic estimate of the required order. I get the same result using fdatool, I must add 2 to the order to be within specified limits. This message was sent using the Comp.DSP web interface on www.DSPRelated.com
Reply by ●September 17, 20052005-09-17
MA wrote:> Hi all! Thanks for helping me! Gotta love this forum and its users:) > > I think I got it right now... > > First: passband ripple between 1 and -1 dB magnitude means in linear scale > magnitude between 1.1220 and 0.8913. The average magnitude in the passband > is then 1.0066. It is important to distinguish this case with passband > ripple given in dB from the case with symmetric deviation from unit > magnitude. Consequently we have to adjust the "desired magnitude" to match > the passband average: > > passbandAverage = (passbandmax+passbandmin)/2 > passbandError = (passbandmax-passbandmin)/2 > > Next, the amplitude corresponding to 40dB below unity lies slightly more > than 40 dB below the new value of the "desired magnitude" in the passband, > so we have to compute a new deviation from 0 in the stopband: > > stopbandmax_dB2 = 20*log10(stopbandmax/passbandAverage) > stopbandmax2 = 10^(stopbandmax_dB2/20) > > We can now send this info to our friend remezord: > > remezord([5 6],[passbandAverage 0],[ passbandError stopbandmax2 ],20) > > > Or we can choose the weights ourselves: > > stopbandweight = passbandError/stopbandmax > weight vector = [1 stopbandweight] > > > > Here below comes the code I used to sort out this, you can copy-paste it > to see a nice plot: > > ------------------------------------------------------------------------ > > > > clear all > > % ------- spec -------- > passbandmax_dB = 4; > passbandmin_dB = -4; > stopbandmax_dB = -40; > % --------------------- > passbandmax = 10^(passbandmax_dB/20) > passbandmin = 10^(passbandmin_dB/20) > stopbandmax = 10^(stopbandmax_dB/20) > passbandPosDev = passbandmax-1 > passbandNegDev = 1-passbandmin > passbandMaxDev = max(passbandPosDev,passbandNegDev) > > passbandAverage = (passbandmax+passbandmin)/2 > passbandError = (passbandmax-passbandmin)/2 > > stopbandmax_dB2 = 20*log10(stopbandmax/passbandAverage) > stopbandmax2 = 10^(stopbandmax_dB2/20) % > > deviation = [ passbandMaxDev stopbandmax ]; > deviation2 = [ passbandError stopbandmax2 ]; > > [N,F,A,W] = remezord([5 6],[1 0],deviation,20); > [N2,F2,A2,W2] = remezord([5 6],[passbandAverage 0],deviation2,20) > > A2 = [passbandAverage passbandAverage 0 0] > W2 = [1 passbandError/stopbandmax] > > b = remez(N,F,A,W); > b2 = remez(N2,F2,A2,W2); > > H = freqz(b,1,2048); > H2 = freqz(b2,1,2048); > > figure > plot([0:2047]/2048,max(20*log10(abs(H)),-70),'k') > hold on > plot([0:2047]/2048,max(20*log10(abs(H2)),-70),'k:') > hold on > plot([0:2047]/2048,20*log10(passbandmax*ones(1,2048)),'k') > hold on > plot([0:2047]/2048,20*log10(passbandmin*ones(1,2048)),'k') > hold on > plot( [0:2047]/2048,20*log10(stopbandmax*ones(1,2048)),'k') > title('Magnitude |H(e^{j \omega})| and |H_{2}(e^{j \omega})|') > xlabel('Normalized frequency [units of \pi]') > ylabel('Magnitude [dB]') > legend('|H(e^{j \omega})|','|H_{2}(e^{j \omega})|')I don't think that I've ever written this before (except for issues of significance), but I think you're trying to be too exact. 30 +/- 10% is the same as 30 +/- 3, and the arithmetic mean is exactly right. Decibels are a ratio, hence multiplicative, so the geometric mean is exact. For small multipliers -- much less than one -- the geometric and arithmetic means are so close that the distinction is ignored. (The derivations of many useful formulas in mechanical and electrical engineering depend on ignoring it.) Consider your computed correction: an average of 1.0066 instead of 1. It is a rare application that is affected by a gain error of 2/3%. We avoid such designs because they lack robustness in use. There's a second issue. Filter specifications aren't intended to be met exactly. You want a stop-band attenuation of _at_least_ 40 dB. Your design might be faulted as overkill if it achieved 50 dB, but 42 or 43 is fine. You can't build a filter with a half tap, and a little margin is good. Two conclusions come from considering the effects of finite-precision arithmetic: exactness beyond a certain point is useless and misleading, and margins of safety are as necessary in calculated structures as they are in bridges. Jerry -- Engineering is the art of making what you want from things you can get. �����������������������������������������������������������������������