minimax and choosing weights

Started by MA September 16, 2005
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
"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
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
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

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. &#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;
"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.com
OK.... 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
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
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
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
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. &#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;