A signal processing exercise

Problem:

We are working on the problem 3 of this picture.

Solution:

For part (a), generating a realization of the random process is simple enough.

Let's do a simple review on white noise. A noise is white if its power spectrum is flat. The power spectrum is the Fourier transform of the autocorrelation function. Note that the term 'autocorrelation function' make sense only when we are talking about wide sense stationary signal.

The fact that the power spectrum is flat implies the autocorrelation function is a delta function, which in turns simply implies the samples are uncorrelated.

Therefore, in order to generate the Gaussian white noise, we simply generate a sequence of uncorrelated Gaussian random variables.

To generate $ x[n] $, we simply generate $ w[n] $ and then run the recursion.

For part (b), we will use the Yule Walker's equation.

The key to the derivation of the Yule Walker's equation is to massage the defining equation to use autocorrelation function, how, just multiply $ x[k] $ on both side and take expectation.

For Matlab, this is easily done with a single call to aryule function. Loved Matlab for that.

For part (c), we interpret the autoregressive random process as an Infinite Impulse Response filter on Gaussian white noise. There is a theorem stating that the power spectrum of the filtered signal is simply the power spectrum of the origina lsignal multiply by the square of the transfer function.

Therefore to determine the power spectrum we simply determine the transfer function.

Using z-transform, we can state the following equation:

$ X(z) = a_1z^{-1}X(z) + a_2z^{-1}X(z) + a_3z^{-3}X(z) + W(z) $

Rearranging, we get $ H(z) = \frac{X(z)}{W(z)} = \frac{1}{1 - a_1z^{-1} - a_2z^{-1} - a_3z^{-3}} $

To obtain the power spectrum, remember the input power spectrum is flat (in fact, identically equals to 1), so we simply ignore that as multiply 1 gives nothing, and take the square of the expression above and evaluate on the unit circle to give coefficients.

For part (d), Bartlett method simply chop the signal into equal non-overlapping windows and then compute the Fourier transform for each window, average them. The point of doing so is to make sure the statistics converge.

Without further ado, this is all the Matlab source code:

function problem3()
  clear; clc;
  
  % part a) generate the samples for the given AR(3) process
  length = 1000;
  old1 = 0;
  old2 = 0;
  old3 = 0;
  x = zeros(1, length);
  w = randn(1, length);
  a1 = 0.8;
  a2 = -0.5;
  a3 = 0.1;
  b0 = 1;
  for i = 1:length
    x(i) = a1 * old1 + a2 * old2 + a3 * old3 + b0 * w(i);
    old3 = old2;
    old2 = old1;
    old1 = x(i);
  end

  % part b) estimate the parameter of the AR process using yule walker
  % equations
  params = aryule(x, 3);
  a1 = params(1)
  a2 = params(2)
  a3 = params(3)

  % part c) plot the power spectral density based on the estimated
  % parameters  
  samples = 0:99;
  z = exp(-2 * pi * sqrt(-1) * samples / 99);
  
  % the process is interpreted as an IIR filter acted on a white noise
  % this the transfer function of the IIR filter  
  h = 1 ./ (1 - a1 ./ z -  a2 ./z ./z - a3 ./z ./z ./z);
  % the power spectral density is given by the square of the transfer 
  % function times to power spectral density of the white noise, which is 
  % identically 1
  p = h .* conj(h);
  
  % estimate the power spectral density based on the Bartlett's method
  num_segments = 10;
  window_size = length / num_segments;
  periodogram = zeros(1, window_size);  
  for i = 0:(num_segments - 1)
    window = x(1, window_size * i + 1: window_size * i + window_size);
    local_fft = fft(window);
    local_fft_magntitude = local_fft .* conj(local_fft);
    periodogram = periodogram + local_fft_magntitude / window_size / num_segments;
  end
  
  figure
  hold on
  plot(p);
  plot(periodogram)