Hi, I am trying to generate noise with an f^-4 PSD, and am having some difficulty getting the spectrum desired. The problem is this: Applying a filter with the transfer function co-efficients b = [1] and a=[1 -1] to white noise creates noise (as expected) with an f^-2 PSD. The rolloff of this filter is 20dB per decade, so it's pretty clear why it works. The theory (as far as I can see) predicts that if I apply the same filter to the data twice, I should get noise with an f^-4 PSD. Alternatively, I should be able to use the double integrator b=[1] and a=[1 -2 1] to achieve the same thing. But, it doesn't work. In MATLAB both y = filter([1], [1 -2 1], randn(1,1e5)) and y = filter([1], [1 -1], filter([1], [1 -1], randn(1, 1e5)) produce noise with a PSD that rolls off at 20dB per decade, not the expected 40. Applying the single integrator (b = [1], a=[1 -1]) produces the expected 20dB rolloff. If I do semilogx(f, 20*log10(abs(fft(y))) the spectrum clearly only rolls off at 20dB per decade. In non-matlab terms, I am applying the same filter to white noise samples twice, which should square the frequency response of the filter, but doesn't seem to in this case. Am I overlooking something obvious, or is this a subtle register length effect? Or something else? Regards Marc
Double integrated white noise
Started by ●February 26, 2008
Reply by ●February 26, 20082008-02-26
Marc Brooker wrote:> Hi, > > I am trying to generate noise with an f^-4 PSD, and am having some > difficulty getting the spectrum desired. The problem is this: > > Am I overlooking something obvious, or is this a subtle register length > effect? Or something else? > > Regards > > MarcLooks like I can answer my own question here. The anomaly comes from the implicit rectangular windowing of the finite number of samples. Applying a suitable window to the data seems to give the expected results. Regards Marc
Reply by ●February 26, 20082008-02-26
Marc, First of all, you have built a digital integrator, which is a unstable system and it has a pole at z=1 (dc freq). For this system the transfer function does not exist. So better approach is to pull the pole slightly inside the unit circle. beta = 15/16 b = 1-beta; % 1 term a = [1 -beta]; % 2 term [hz,f] = freqz(b,a,1024,2); plot(f, 20*log10(abs(hz))); [hz, f] = freqz(conv(b,b), conv(a,a), 1024, Fs); plot(f, 20*log10(abs(hz)); in digital domain, highest frequency possible is Fs/2. and for this reason you may not see 20db/decade rolloff, mainly due to frequency warping effects. for the above example you will get %--------------------------------------------------- f = [1 2 3 4 5 6 7 8 9 10]*Fs/20; a1db = [-14 -20 -23 -25 -27 -28 -28.8 -29 -29.7 -30] % single stage case a2db = [-28 -39 -46 -50 -54 -56 -57.7 -58.7 -59.4 -60] % 2 stage case hence in first case we get 16db/decade and in second case we get 32db/decade. the rolloff is not linear. Regards Bharat Pathak Arithos Designs www.Arithos.com
Reply by ●February 26, 20082008-02-26
>Looks like I can answer my own question here. The anomaly comes from the> implicit rectangular windowing of the finite number of samples. >Applying a suitable window to the data seems to give the expectedresults. Marc, you are feeding random data. hence there is no concept of spectral leakage as spectral leakage will always be there if u window or not. windowing is helpful if you have sinetones at input and then u chop-off some portion of it, and try to compute FFT. Regards Bharat Pathak Arithos Designs www.Arithos.com
Reply by ●February 26, 20082008-02-26
bharat pathak wrote:>> Looks like I can answer my own question here. The anomaly comes from the > >> implicit rectangular windowing of the finite number of samples. >> Applying a suitable window to the data seems to give the expected > results. > > Marc, > > you are feeding random data. hence there is no concept of spectral > leakage as spectral leakage will always be there if u window or not. > > windowing is helpful if you have sinetones at input and then u > chop-off some portion of it, and try to compute FFT. >I think in this case windowing _is_ indicated. Mostly because Marc clearly saw better results when he applied it, but also because of the nature of the sample he's FFT-ing. If you take an FFT of white noise, then spectral leakage is immaterial -- if all samples are independent, then there will be no artificial discontinuity between the 1st sample and the Nth. However, if you take an FFT of colored noise, which should have some correlation between adjacent samples, then the step between the 1st and Nth samples will, indeed, splatter over the spectrum at 20dB/decade. Moreover, the lower the frequency content of the noise, the worse this leakage will be. This gets particularly bad when you try this with data from integrated white noise, because the Nth sample is a random walk away from the 1st; with such high variance your expected step is large. With a double integrator your expected step is _huge_. So, Marc -- window the data, with my blessing. It may even me more valid to remove the mean, the trend, and the x^2 trend from the noise before doing your FFT, then augment the FFT output with an estimate of the effect of the mean, trend and x^2 trend on the spectrum. -- Tim Wescott Wescott Design Services http://www.wescottdesign.com Do you need to implement control loops in software? "Applied Control Theory for Embedded Systems" gives you just what it says. See details at http://www.wescottdesign.com/actfes/actfes.html
Reply by ●February 26, 20082008-02-26
bharat pathak wrote: (context restored) > Marc Brooker wrote: >> Hi, >> >> I am trying to generate noise with an f^-4 PSD, and am having some >> difficulty getting the spectrum desired. The problem is this: >> >> Am I overlooking something obvious, or is this a subtle register >> length effect? Or something else? >> >> Regards >> >> Marc> Marc, > > First of all, you have built a digital integrator, which > is a unstable system and it has a pole at z=1 (dc freq). > For this system the transfer function does not exist. > > So better approach is to <snip> >I must be disagreeable today -- I'm certainly disagreeing enough. Depending on how you constrain your z transform a transfer function with a pole outside just about any region you wish to define on the complex plane "doesn't exist". But then, depending on how you constrain your z transform you can find a region of convergence for just about any pole anywhere. With careful dancing around the math and a creative view of causality you can get just about any z transform that consists of exponents, steps and ramps to converge. But that's not important. What's important is that you can _use_ a transfer function with a pole at z=1 to make accurate predictions, both theoretically and through experimentation. So no matter how creatively you have to state your case to a mathematician, when you're doing engineering work you can just take the thing at face value. The only thing that you _must_ pay attention to is that if you are going to do like Marc did and start using an integrator (or a double integrator) all by itself you have to account for the instability, and the resulting behavior that for most inputs the outputs will go to infinity. Marc did this with windowing, which is probably about as proper as you can get given that the system as a whole is much more amenable to theoretical analysis than to experiment. -- Tim Wescott Wescott Design Services http://www.wescottdesign.com Do you need to implement control loops in software? "Applied Control Theory for Embedded Systems" gives you just what it says. See details at http://www.wescottdesign.com/actfes/actfes.html
Reply by ●February 26, 20082008-02-26
Tim, Good that you are disagreeing. Atleast this way I will understand better. An integrator is an unstable system atleast for z=1. So any small amount of DC at input is going to make it to go into oscillations. Only the period of oscillations will change based on how low is the DC value. How does windowing of gaussian data helps in removingof DC from the incoming signal to avoid oscillations. In order to remove DC completely from the incoming signal, you need a differentiator upfront. i.e. hz = [1 -1] and hence if you look at all CIC systems they are a combination of differentiator and integrators. To remove DC completely you need a zero at z = 1, which essentially a differentiator does. It offers -infinite db attenuation at dc. Regards Bharat Pathak>I must be disagreeable today -- I'm certainly disagreeing enough. > >Depending on how you constrain your z transform a transfer function with>a pole outside just about any region you wish to define on the complex >plane "doesn't exist". > >But then, depending on how you constrain your z transform you can find a>region of convergence for just about any pole anywhere. With careful >dancing around the math and a creative view of causality you can get >just about any z transform that consists of exponents, steps and ramps >to converge. > >But that's not important. What's important is that you can _use_ a >transfer function with a pole at z=1 to make accurate predictions, both >theoretically and through experimentation. So no matter how creatively >you have to state your case to a mathematician, when you're doing >engineering work you can just take the thing at face value. > >The only thing that you _must_ pay attention to is that if you are going>to do like Marc did and start using an integrator (or a double >integrator) all by itself you have to account for the instability, and >the resulting behavior that for most inputs the outputs will go to >infinity. Marc did this with windowing, which is probably about as >proper as you can get given that the system as a whole is much more >amenable to theoretical analysis than to experiment. > >-- > >Tim Wescott >Wescott Design Services >http://www.wescottdesign.com > >Do you need to implement control loops in software? >"Applied Control Theory for Embedded Systems" gives you just what itsays.>See details at http://www.wescottdesign.com/actfes/actfes.html >
Reply by ●February 26, 20082008-02-26
bharat pathak wrote:> Tim, > > Good that you are disagreeing. Atleast this way I will understand > better. > > An integrator is an unstable system atleast for z=1. So any small > amount of DC at input is going to make it to go into oscillations.Zero Hz., if you cpnsider that "oscillation".> Only the period of oscillations will change based on how low is the > DC value.Any DC integrated forever goes to infinity. There is only oscillation _after_ a nonlinearity kicks in, and that's outside the realm of linear analysis.> How does windowing of gaussian data helps in removingof DC from the > incoming signal to avoid oscillations. > > In order to remove DC completely from the incoming signal, you > need a differentiator upfront. i.e. hz = [1 -1] and hence if you > look at all CIC systems they are a combination of differentiator > and integrators. To remove DC completely you need a zero at z = 1, > which essentially a differentiator does. It offers -infinite db > attenuation at dc.Jerry -- Engineering is the art of making what you want from things you can get. �����������������������������������������������������������������������
Reply by ●February 26, 20082008-02-26
On Feb 26, 2:11 am, Marc Brooker <myrealn...@gmail.com> wrote:> Marc Brooker wrote: > > Hi, > > > I am trying to generate noise with an f^-4 PSD, and am having some > > difficulty getting the spectrum desired. The problem is this: > > > Am I overlooking something obvious, or is this a subtle register length > > effect? Or something else? > > > Regards > > > Marc > > Looks like I can answer my own question here. The anomaly comes from the > implicit rectangular windowing of the finite number of samples. > Applying a suitable window to the data seems to give the expected results. > > Regards > > MarcMarc An example of evaluating windows for measuring noise shapes is: Estimation of the In-Band Delta-Sigma Noise Power Based on Windowed Data Nunzi, E. Carbone, P. Petri, D. Dipt. di Ingegneria Elettronica e dell'Informazione, Univ. degli Studi di Perugia IEEE Transactions on Instrumentation and Measurement, Dec. 2006, Volume: 55 , Issue: 6, page(s): 2221 - 2226 To achieve high rolloff rates in window response consider the maximum rolloff windows that Nuttall discusses in: Some windows with very good sidelobe behavior IEEE Trans. Acoust., Speech, Signal Processing, vol. 29, pp. 84 - 91, February 1981 Albert H. Nuttall These are the binomial coefficient windows. Von Hann is the two coefficient version. Steeper slopes come with more coefficients. Dale B. Dalrymple http;//dbdimages.com
Reply by ●February 26, 20082008-02-26
On Feb 26, 6:22 pm, "bharat pathak" <bha...@arithos.com> wrote:> Tim, > > Good that you are disagreeing. Atleast this way I will understand > better. > > An integrator is an unstable system atleast for z=1. So any small > amount of DC at input is going to make it to go into oscillations. > Only the period of oscillations will change based on how low is the > DC value. > > How does windowing of gaussian data helps in removingof DC from the > incoming signal to avoid oscillations. > > In order to remove DC completely from the incoming signal, you > need a differentiator upfront. i.e. hz = [1 -1] and hence if you > look at all CIC systems they are a combination of differentiator > and integrators. To remove DC completely you need a zero at z = 1, > which essentially a differentiator does. It offers -infinite db > attenuation at dc. > > Regards > Bharat Pathak > ...The concern about instability at DC is relevant for a continuous system of infinite extent with a finite DC input.. Marc, however, has chosen a finite signal with roughly zero DC content. It is useful to have people around who know how the maths work. It is also useful to remember when the maths don't apply. Dale B. Dalrymple
