# Filter coefficients for colored noise Voss Algorithm

If I use a generic filter for generating colored noise like pink,brown,white then how do I modify this statement and how do I know what are the coefficients to be used in AR model for different noise.Signal is white gaussian noise and y gives different noise based on the value of a.

• This is not a well-phrased question. If you have MATLAB programming questions, ask them on Stack Overflow. If you have questions on the spectral shape of different types of noise and how to generate them, then I would reword your question to be more specific in that direction. Jan 8 '12 at 15:08
• I have not asked about programming.I have specifically asked about the coefficients of the numerator and denominator required to design the filters.In some answers in Mathworks mathworks.com/matlabcentral/newsreader/view_thread/292182 they mention the value of a=-0.9,i wanted to know on what basis they do that and hence what shall be the values for other colored noises.
– sts
Jan 8 '12 at 16:48
• Even just the conceptual answer for some types of colored noise can be more complex than you might expect (see e.g. dsp.stackexchange.com/questions/322/…). Jan 8 '12 at 19:44
• There must be some values available for a,b for colored noises like the one shown in the link : mathworks.com/matlabcentral/newsreader/view_thread/292182 without going into the mathematical concept else how did the responder suggest that value?
– sts
Jan 8 '12 at 20:10

Here, $a$ is a filter coefficient for a simple single-pole recursive low-pass filter. If you look at the documentation for MATLAB's filter() function, you'll see this corresponds to the difference equation $y(n) = x(n) - a y(n-1)$. Let's analyze the filter in the $z$-domain to get its frequency response (just remembering that $z^{-1}$ corresponds to a single sample delay:

$$y(n) + a y(n-1) = x(n)$$ $$Y(z)*(1 + a z^{-1}) = X(z)$$ $$\frac{Y(z)}{X(z)} = \frac{1}{1 + a z^{-1}} = H(z)$$

That's the $z$-domain transfer function. We can substitute $z = e^{j\omega}$ where $\omega$ is the normalized frequency, where 0 is DC and $\pi$ corresponds to the Nyquist frequency of whatever sampling rate you're at. Now you should realize that if you use the same digital coefficients at different sampling rates, you will get filtered noise that varies depending on the sampling rate. If your sampling rate is $F_s$, then $\omega$ corresponds to the frequency $\dfrac{\omega F_s}{2 \pi}$. Let's find the frequency response:

$$G^2(\omega) = |H(e^{j\omega})|^2 = \left| \frac{1}{1 + a e^{-j\omega}} \right|^2$$

$$G^2(\omega) = \frac{1}{|1 + a(\cos(-\omega) + j \sin(-\omega))|^2}$$

$$G^2(\omega) = \frac{1}{|(1 + a \cos(\omega)) - j \sin(\omega))|^2}$$

$$G^2(\omega) = \frac{1}{(1 + a \cos(\omega))^2 + \sin^2(\omega)}$$

$$G^2(\omega) = \frac{1}{1 + 2 a \cos(\omega) + a^2 \cos^2(\omega) + \sin^2(\omega)}$$

$$G^2(\omega) = \frac{1}{1 + 2 a \cos(\omega) + a^2 \cos^2(\omega) + (1 - \cos^2(\omega))}$$

$$G^2(\omega) = \frac{1}{2 + 2 a \cos(\omega) + (a^2 - 1) \cos^2(\omega)}$$

$$G(\omega) = \frac{1}{\sqrt{2 + 2 a \cos(\omega) + (a^2 - 1) \cos^2(\omega)}}$$

Let's plot the frequency response with $a = -0.9$: WolframAlpha plot. We can also take 20 Log10() to get the frequency response in decibels: WolframAlpha plot .

Now we can find out some things about this filter. For example, its gain at DC is greater than unity! $G(0) = \dfrac{1}{\sqrt{1 + 2a + a^2}}$. When $a = -0.9$, the filter increases the DC gain by a factor of 10, or 20 dB. Often filters are characterized in terms of their -3 dB point, which is where their magnitude decreases by $\sqrt{2}$ from their peak. Let's solve for this point:

$$\frac{1} { \sqrt{2} \sqrt{1 + 2a + a^2} } = \frac{1}{\sqrt{2 + 2 a \cos(\omega) + (a^2 - 1) \cos^2(\omega)}}$$

$$\sqrt{2} \sqrt{1 + 2a + a^2} = \sqrt(2 + 2 a \cos(\omega) + (a^2 - 1) \cos^2(\omega))$$

$$2 + 4a + 2 a^2 = 2 + 2 a \cos(\omega) + (a^2 - 1) \cos^2(\omega)$$

Now treat this as a quadratic in $\cos(\omega)$, let's substitute $x = \cos(\omega)$ and use the quadratic formula:

$$(a^2 - 1) x + 2ax + (-4a - 2 a^2) = 0$$ $$x = \frac{-2a \pm \sqrt( 4a^2 - 4(a^2-1)(-4a - 2 a^2)) }{2 (a^2 - 1)}$$

Now, we want the solution corresponding to positive frequencies, and it turns out that choosing the result with the minus sign gives us that. So, we can take the inverse cosine, and get our -3 dB frequency:

$$\omega = \arccos\left( \frac{-2a - \sqrt{4a^2 - 4(a^2-1)(-4a - 2 a^2)}}{2 (a^2 - 1)} \right)$$

Let's plug it in for $a = -0.9$ and verify:

$w = 0.0958385$

$H(z) = 7.07107 = \dfrac{10}{\sqrt{2}}$

It might be useful to go the other way, to design a filter by choosing a value of a to have some desired -3 dB point:

$$2 + 4a + 2 a^2 = 2 + 2 a \cos(\omega) + (a^2 - 1) \cos^2(\omega)$$

$$2 + 4a + 2 a^2 = 2 + 2 a \cos(\omega) + a^2 cos^2(\omega) - \cos^2(\omega)$$

$$(2 - \cos^2(\omega)) a^2 + (4 - 2 \cos(\omega)) a + \cos^2(w) = 0$$

$$\frac{ -(4 - 2 \cos(\omega)) +/- \sqrt{ (4 - 2 \cos^2(\omega))^2 - 4(2 - \cos^2(\omega))(\cos^2(\omega)) } }{2 (2 - \cos^2(\omega))}$$

And again we can verify by plugging in $\omega = 0.0958385$, and we'll see that we get $a = -0.9$.

So, for example, you could use these formulas to pick values of $a$ at different sampling rates, so that you have the same frequency -3 dB point. You can use the formula for gain at DC, for example, to normalize the filter so its gain at DC is unity. You can take the formula for gain, convert it to decibels, and see that it approaches roughly a -6 dB/oct slope, which is why it's being used to shape noise like this.

• WOW !!This mathematical analysis would take atleast couple of days for me to digest!Thank you for the concept,I really appreciate.So,if you would be kind enough to also tell what should be the values of a (if you remember,since it will help me to double check after the calculations) for different noises like Sine Weiner,red,violet,pink etc.
– sts
Jan 8 '12 at 18:27
• Also,there are these terms "bounded","unbounded" noise/signals. Under what category does White Gaussian noise, Gaussian Colored noise and non Gaussian colored noise fall into?
– sts
Jan 8 '12 at 18:45
• That sort of filter won't generate exactly those types of noise. You'll have to come up with some way to approximate it by choosing a`. I think a good first step would be to normalize the filter's gain to unity at DC, and then to pick some point from the desire noise profile's frequency response, i.e. the 3 dB point, and try to match that. Honestly I don't really know anything about approximating different kinds of noise; sorry I don't have more information to offer. Jan 8 '12 at 19:09
• Hi,I could not follow the last point which u have made >>It might be useful to go the other way, to design a filter by choosing a value of a to have some desired -3 dB point.Which formula to use for this?
– sts
Jan 8 '12 at 19:43
• How is it possible to know which noise has what frequency like in the answer mentioned -3dB.
– sts
Jan 8 '12 at 21:03