Generate color noise using Auto-Regressive (AR) model

by

Key focus: Learn how to generate color noise using auto regressive (AR) model. Apply Yule Walker equations for generating power law noises: pink noise, Brownian noise.

Auto-Regressive (AR) model

An uncorrelated Gaussian random sequence x[n] can be transformed into a correlated Gaussian random sequence y[n] using an AR time-series model. If a time series random sequence is assumed to be following an auto-regressive AR(N) model of form,

where x[n] is the uncorrelated Gaussian sequence of zero mean and variance \sigma_x^2, the natural tendency is to estimate the model parameters a_1, a_2,\cdots,a_N. Least Squares method can be applied here to find the model parameters, but the computations become cumbersome as the order N increases. Fortunately, the AR model coefficients can be solved for using Yule Walker equations.

Yule Walker equations relate auto-regressive model parameters a_1, a_2,\cdots,a_N to the auto-correlation R_{yy}[k] of random process y[n]. Finding the model parameters using Yule-Walker equations, is a two step process:

1. Given y[n], estimate auto-correlation of the process R_{yy}[n]. If R_{yy}[n] is already specified as a function, utilize it as it is (see auto-correlation equations for Jakes spectrum or Doppler spectrum in section 11.3.2 in the book).

2. Solve Yule Walker equations to find the model parameters a_k \; \left(k=0,1,...,N\right) and the noise sigma \sigma_x^2.

Yule-Walker equations

Yule-Walker equations can be compactly written as

Yule Walker equation
Equation (2) Yule Walker equation

Written in matrix form, the Yule-Walker equations that comprises of a set of N linear equations and N unknown parameters.

Representing equation (3) in a compact form,

\mathbf{R_{yy}}\mathbf{a}=-\mathbf{r_{yy}} \quad \quad \quad \quad (4)

The AR model parameters \mathbf{a} can be found by solving

\mathbf{a}=-\mathbf{R_{yy}}^{-1}\mathbf{r_{yy}} \quad \quad \quad \quad (5)

After solving for the model parameters a_k, the noise variance \sigma_x^2 can be found by applying the estimated values of a_k in equation (2) by setting l=0. The aryule command (in Matlab and Python’s spectrum package) efficiently solves the Yule-Walker equations using Levinson Algorithm [1][2]. Once the model parameters a_k are obtained, the AR model can be implemented as an \emph{infinte impulse response (IIR)} filter of form

H(z) = \dfrac{1}{1 + \sum \limits_{k=1}^{N} a_k z^{-k}} \quad \quad \quad \quad (6)

Example: 1/f^\alpha power law noise generation

The power law 1/f^\alpha in the power spectrum characterizes the fluctuating observables in many natural systems. Many natural systems exhibit some 1/f^\alpha noise which is a stochastic process with a power spectral density having a power exponent that can take values -2 \geq \alpha \leq 2 . Simply put, 1/f^{\alpha} noise is a colored noise with a power spectral density of 1/f^\alpha over its entire frequency range.

S(f) = \frac{1}{|f|^\alpha} \quad \quad (7)

The 1/f^\alpha noise can be classified into different types based on the value of \alpha.

Violet noise – \alpha = -2, the power spectral density is proportional to f^2.
Blue noise – \alpha = -1, the power spectral density is proportional to f.
White noise – \alpha = 0, the power spectral density is flat across the whole spectrum.
Pink noise – \alpha = 1, the power spectral density is proportional to 1/f, i.e, it decreases by 3 \; dB per octave with increase in frequency.
Brownian noise – \alpha = 2, the power spectral density is proportional to 1/f^2, therefore it decreases by 6 \; dB per octave with increase in frequency.

The power law noise can be generated by sequencing a zero-mean white noise through an auto-regressive (AR) filter of order N:

\sum_{k=0}^{k=N} a_k y[n - k] = x[n] \quad \quad \quad \quad (8)

where, x[n] is a zero-mean white noise process. Referring the AR generation method described in [3], the coefficients of the AR filter can be generated as

\begin{aligned} a_0 &= 1 \\ a_k &= \left(k - 1 - \frac{\alpha}{2}\right) \frac{a_k - 1}{k}\end{aligned} \quad \quad \quad (9)

which can be implemented as an infinite impulse response (IIR) filter using the filter transfer function described in equation (6).

The following script implements this method and the sample results are plotted in the next Figure.

Refer the book for the Matlab code

Simulated color noise samples and their PSD estimates: pink noise (α =1) and Brown noise (α =2)
Figure 1: Simulated color noise samples and their PSD estimates: pink noise (α =1) and Brown noise (α =2)

Rate this article: PoorBelow averageAverageGoodExcellent (4 votes, average: 4.75 out of 5)

References

[1] Gene H. Golub, Charles F. Van Loan, Matrix Computations, ISBN-9780801854149, Johns Hopkins University Press, 1996, p. 143.↗
[2] J. Durbin, The fitting of time series in models, Review of the International Statistical Institute, 28:233-243, 1960.↗
[3] Kasdin, N.J. Discrete Simulation of Colored Noise and Stochastic Processes and 1/f^\alpha Power Law Noise Generation, Proceedings of the IEEE, Vol. 83, No. 5, 1995, pp. 802-827.↗

Rate this article: PoorBelow averageAverageGoodExcellent (4 votes, average: 4.75 out of 5)

Books by the author

Wireless Communication Systems in Matlab
Wireless Communication Systems in Matlab
Second Edition(PDF)
PoorBelow averageAverageGoodExcellent (179 votes, average: 3.62 out of 5)

Digital modulations using Python
Digital Modulations using Python
(PDF ebook)
PoorBelow averageAverageGoodExcellent (134 votes, average: 3.56 out of 5)

digital_modulations_using_matlab_book_cover
Digital Modulations using Matlab
(PDF ebook)
PoorBelow averageAverageGoodExcellent (135 votes, average: 3.63 out of 5)
Hand-picked Best books on Communication Engineering
Best books on Signal Processing

Topics in this chapter

Random Variables - Simulating Probabilistic Systems
● Introduction
● Plotting the estimated PDF
● Univariate random variables
 □ Uniform random variable
 □ Bernoulli random variable
 □ Binomial random variable
 □ Exponential random variable
 □ Poisson process
 □ Gaussian random variable
 □ Chi-squared random variable
 □ Non-central Chi-Squared random variable
 □ Chi distributed random variable
 □ Rayleigh random variable
 □ Ricean random variable
 □ Nakagami-m distributed random variable
● Central limit theorem - a demonstration
● Generating correlated random variables
 □ Generating two sequences of correlated random variables
 □ Generating multiple sequences of correlated random variables using Cholesky decomposition
● Generating correlated Gaussian sequences
 □ Spectral factorization method
 □ Auto-Regressive (AR) model

2 thoughts on “Generate color noise using Auto-Regressive (AR) model”

  1. Interesting to see the link between autocorrelation and the exponent of the PSD.

    I implemented the code in the book on page 99.

    How would I go about characterizing alpa from the PSD?
    If I fit a line in loglog space to the PSD, this line is severely biased by the noise lower end of the PSD.

    Example code:
    log_f = 10log10(F(2:end));
    log_p = 10    log10(Pyy(2:end));
    Const = polyfit(log_f,log_p,1);

    slope = Const(1);
    offset = Const(2);
    Yp = polyval(Const, log_f);

    Do you have any example of a robust estimate of the slope and offset following the example noise generation?

    Reply

    • This method can only be applied for noises that has certain statistical properties. The relationship between alpha, PSD and the AR model cannot be obtained by a simple line-fit model. It has to be estimated from empirical data.

      The noise should exhibit fractal phenomenon (long memory and self-similarity) for this method to work. Following references may help you further.

      [1] Jan Beran, “Statistics for Long-Memory Processes “, ISBN-13: 978-0412049019, Chapman and Hall/CRC; 1 edition (October 1, 1994) URL: https://amzn.to/2w9I0NO
      [2] Stadnitski, Tatjana. “Measuring fractality.” Frontiers in physiology vol. 3 127. 7 May. 2012, doi:10.3389/fphys.2012.00127 URL: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3345945/

      Reply

Post your valuable comments !!!