In deterministic analysis, we transform a known signal into frequency domain. In stochastic analysis, the signal itself is random, so we characterize its ensemble second-order structure through \( R_{xx} \) and \( S_{xx} \).

\[ R_{xx}(\theta)=\mathbb{E}\!\left[(x(t)-\mu_x)(x(t+\theta)-\mu_x)\right], \qquad R_{xx}[m]=\mathbb{E}\!\left[(x[n]-\mu_x)(x[n+m]-\mu_x)\right] \]

For WSS processes these functions depend only on lag. At zero lag they equal the variance.

flowchart TD
  A["Random process x(t) or x[n]"] --> B["Check WSS assumptions and remove mean"]
  B --> C["Compute autocorrelation Rxx(lag)"]
  C --> D["Transform to PSD Sxx(freq)"]
  D --> E["Interpret bandwidth and dominant content"]
  E --> F["Build white / colored noise models"]
        

Even symmetry (real WSS):

\[ R_{xx}(-\theta)=R_{xx}(\theta), \qquad R_{xx}[-m]=R_{xx}[m] \]

Proof: For continuous time,

\[ R_{xx}(-\theta)=\mathbb{E}[x(t)x(t-\theta)] =\mathbb{E}[x(t+\theta)x(t)] =\mathbb{E}[x(t)x(t+\theta)] =R_{xx}(\theta) \]

Positive semidefinite property: for any real coefficients \( a_1,\dots,a_p \) and lags \( \theta_1,\dots,\theta_p \),

\[ \sum_{i=1}^p\sum_{j=1}^p a_i a_j R_{xx}(\theta_i-\theta_j) \ge 0 \]

Proof: Define \( z=\sum_i a_i x(t+\theta_i) \). Then \( \mathbb{E}[z^2] \ge 0 \), and expansion gives the claim.

\[ \mathbb{E}[z^2] =\sum_{i=1}^p\sum_{j=1}^p a_i a_j \mathbb{E}[x(t+\theta_i)x(t+\theta_j)] =\sum_{i=1}^p\sum_{j=1}^p a_i a_j R_{xx}(\theta_i-\theta_j) \]

The PSD is the Fourier transform of the autocorrelation. Continuous-time:

\[ S_{xx}(\omega)=\int_{-\infty}^{\infty} R_{xx}(\theta)e^{-j\omega\theta}\,d\theta, \qquad R_{xx}(\theta)=\frac{1}{2\pi}\int_{-\infty}^{\infty} S_{xx}(\omega)e^{j\omega\theta}\,d\omega \]

Discrete-time:

\[ S_{xx}(e^{j\omega})=\sum_{m=-\infty}^{\infty} R_{xx}[m]e^{-j\omega m}, \qquad R_{xx}[m]=\frac{1}{2\pi}\int_{-\pi}^{\pi} S_{xx}(e^{j\omega})e^{j\omega m}\,d\omega \]

A direct consequence is \( S_{xx}(e^{j\omega})\ge 0 \), and for zero-mean processes, the total spectral area gives variance:

\[ R_{xx}[0]=\operatorname{Var}(x[n])=\frac{1}{2\pi}\int_{-\pi}^{\pi} S_{xx}(e^{j\omega})\,d\omega \]

Ideal continuous-time white noise is a generalized process:

\[ R_{ww}(\theta)=\frac{N_0}{2}\delta(\theta), \qquad S_{ww}(\omega)=\frac{N_0}{2} \]

It has flat PSD on all frequencies, which makes it useful analytically but idealized physically.

Discrete-time white noise with variance \( \sigma_w^2 \):

\[ \mathbb{E}[w[n]]=0, \qquad R_{ww}[m]=\sigma_w^2\delta[m], \qquad S_{ww}(e^{j\omega})=\sigma_w^2 \]

Band-limited white noise (practical simulation model):

\[ S_{ww}(\omega)= \begin{cases} S_0, & |\omega| \le \Omega_c \\ 0, & |\omega| > \Omega_c \end{cases} \]

\[ R_{ww}(\theta)=\frac{1}{2\pi}\int_{-\Omega_c}^{\Omega_c} S_0 e^{j\omega\theta}\,d\omega =\frac{S_0}{\pi}\frac{\sin(\Omega_c\theta)}{\theta} \]

For a finite record \( x[0],\dots,x[N-1] \), the ensemble moments are replaced by estimators.

\[ \hat{R}_{xx}^{(b)}[m]=\frac{1}{N}\sum_{n=0}^{N-1-m} x[n]x[n+m], \qquad \hat{R}_{xx}^{(u)}[m]=\frac{1}{N-m}\sum_{n=0}^{N-1-m} x[n]x[n+m] \]

The periodogram estimator is

\[ I_N(e^{j\omega})=\frac{1}{N}\left|\sum_{n=0}^{N-1}x[n]e^{-j\omega n}\right|^2 \]

In practice, Welch/Bartlett averaging or lag-window smoothing is used to reduce variance.

\[ \hat{S}_{xx}(e^{j\omega})=\sum_{m=-M}^{M} w[m]\hat{R}_{xx}[m]e^{-j\omega m} \]

flowchart TD
  A["Finite record x[n]"] --> B["Subtract mean"]
  B --> C["Estimate Rxx[m]"]
  C --> D["Window / truncate lags"]
  D --> E["FFT to estimate PSD"]
  B --> F["Alternative: periodogram / Welch"]
  E --> G["Variance and bandwidth interpretation"]
  F --> G
        

Let \( x[n]=w[n]+b\,w[n-1] \) with white input \( R_{ww}[m]=\sigma_w^2\delta[m] \). Then

\[ R_{xx}[m]=\sigma_w^2\Big((1+b^2)\delta[m]+b\delta[m-1]+b\delta[m+1]\Big) \]

\[ S_{xx}(e^{j\omega})=\sigma_w^2(1+b^2+2b\cos\omega) =\sigma_w^2\left|1+be^{-j\omega}\right|^2 \]

This example is simple enough for full derivation and rich enough to illustrate correlation-to-PSD mapping.

Chapter17_Lesson2.py

# Chapter17_Lesson2.py
# Autocorrelation, Power Spectral Density, and White Noise Models

import numpy as np
import matplotlib.pyplot as plt

def autocorr_biased(x, M):
    x = np.asarray(x, float) - np.mean(x)
    N = len(x)
    return np.array([np.dot(x[:N-m], x[m:]) / N for m in range(M + 1)])

rng = np.random.default_rng(17)
N, sigma, b = 2048, 2.0, 0.6
w = sigma * rng.standard_normal(N)                 # white noise
x = np.r_[w[0], w[1:] + b * w[:-1]]               # colored: x[n]=w[n]+b w[n-1]

M = 80
Rw = autocorr_biased(w, M)
Rx = autocorr_biased(x, M)

# PSD from autocorrelation (Wiener-Khinchin route)
Re = np.r_[Rx, Rx[-2:0:-1]]                       # even extension
Sx = np.real(np.fft.fft(Re))
omega = 2 * np.pi * np.arange(len(Sx)) / len(Sx)

# Theory for x[n]=w[n]+b w[n-1]
Sx_th = sigma**2 * (1 + b**2 + 2 * b * np.cos(omega))

print("R_w[0] (variance estimate) =", Rw[0])
print("R_x[0] (variance estimate) =", Rx[0])
print("Theoretical var(x)        =", sigma**2 * (1 + b**2))

plt.figure()
plt.plot(Rw, label="White")
plt.plot(Rx, label="Colored")
plt.xlabel("Lag m"); plt.ylabel("Autocorrelation"); plt.grid(True); plt.legend()

plt.figure()
half = len(omega) // 2
plt.plot(omega[:half], Sx[:half], label="Estimated PSD")
plt.plot(omega[:half], Sx_th[:half], label="Theory")
plt.xlabel("rad/sample"); plt.ylabel("PSD"); plt.grid(True); plt.legend()
plt.show()

Chapter17_Lesson2.cpp

// Chapter17_Lesson2.cpp
// Autocorrelation and PSD (compact C++17 example)

#include <cmath>
#include <complex>
#include <iostream>
#include <numeric>
#include <random>
#include <vector>

std::vector<double> autocorrBiased(const std::vector<double>& x, int M) {
    int N = (int)x.size();
    double mu = std::accumulate(x.begin(), x.end(), 0.0) / N;
    std::vector<double> xc(N), R(M + 1, 0.0);
    for (int i = 0; i < N; ++i) xc[i] = x[i] - mu;
    for (int m = 0; m <= M; ++m) {
        for (int n = 0; n < N - m; ++n) R[m] += xc[n] * xc[n + m];
        R[m] /= N;
    }
    return R;
}

std::vector<double> evenExtend(const std::vector<double>& R) {
    int M = (int)R.size() - 1;
    std::vector<double> e;
    for (int i = 0; i <= M; ++i) e.push_back(R[i]);
    for (int i = M - 1; i >= 1; --i) e.push_back(R[i]);
    return e;
}

std::vector<double> dftReal(const std::vector<double>& x) {
    int N = (int)x.size();
    std::vector<double> S(N, 0.0);
    for (int k = 0; k < N; ++k) {
        std::complex<double> Xk(0.0, 0.0);
        for (int n = 0; n < N; ++n) {
            double a = -2.0 * M_PI * k * n / N;
            Xk += x[n] * std::complex<double>(std::cos(a), std::sin(a));
        }
        S[k] = Xk.real(); // real because autocorrelation extension is even
    }
    return S;
}

int main() {
    const int N = 512, M = 20;
    const double sigma = 1.5, b = 0.6;
    std::mt19937 gen(17);
    std::normal_distribution<double> gauss(0.0, sigma);

    std::vector<double> w(N), x(N);
    for (int n = 0; n < N; ++n) w[n] = gauss(gen);
    x[0] = w[0];
    for (int n = 1; n < N; ++n) x[n] = w[n] + b * w[n - 1];

    auto Rw = autocorrBiased(w, M);
    auto Rx = autocorrBiased(x, M);
    auto Sx = dftReal(evenExtend(Rx));

    std::cout << "R_w[0] = " << Rw[0] << "\\n";
    std::cout << "R_x[0] = " << Rx[0] << "\\n";
    std::cout << "Theo R_x[0] = " << sigma * sigma * (1.0 + b * b) << "\\n";
    std::cout << "First 8 PSD bins from autocorrelation-DFT:\\n";
    for (int k = 0; k < 8; ++k) std::cout << "k=" << k << "  " << Sx[k] << "\\n";
    return 0;
}

Chapter17_Lesson2.java

// Chapter17_Lesson2.java
// Autocorrelation and PSD (compact Java example)

import java.util.Random;

public class Chapter17_Lesson2 {
    static double[] autocorrBiased(double[] x, int M) {
        int N = x.length;
        double mu = 0.0;
        for (double v : x) mu += v;
        mu /= N;

        double[] xc = new double[N];
        for (int i = 0; i < N; i++) xc[i] = x[i] - mu;

        double[] R = new double[M + 1];
        for (int m = 0; m <= M; m++) {
            for (int n = 0; n < N - m; n++) R[m] += xc[n] * xc[n + m];
            R[m] /= N;
        }
        return R;
    }

    static double[] evenExtend(double[] R) {
        int M = R.length - 1;
        double[] e = new double[2 * M];
        for (int i = 0; i <= M; i++) e[i] = R[i];
        for (int i = M - 1; i >= 1; i--) e[2 * M - i] = R[i];
        return e;
    }

    static double[] dftReal(double[] x) {
        int N = x.length;
        double[] S = new double[N];
        for (int k = 0; k < N; k++) {
            double re = 0.0, im = 0.0;
            for (int n = 0; n < N; n++) {
                double a = -2.0 * Math.PI * k * n / N;
                re += x[n] * Math.cos(a);
                im += x[n] * Math.sin(a);
            }
            S[k] = re; // autocorrelation extension is even, so imaginary part is near zero
        }
        return S;
    }

    public static void main(String[] args) {
        int N = 512, M = 20;
        double sigma = 1.2, b = 0.6;
        Random rng = new Random(17);

        double[] w = new double[N];
        for (int n = 0; n < N; n++) w[n] = sigma * rng.nextGaussian();

        double[] x = new double[N];
        x[0] = w[0];
        for (int n = 1; n < N; n++) x[n] = w[n] + b * w[n - 1];

        double[] Rw = autocorrBiased(w, M);
        double[] Rx = autocorrBiased(x, M);
        double[] Sx = dftReal(evenExtend(Rx));

        System.out.println("R_w[0] = " + Rw[0]);
        System.out.println("R_x[0] = " + Rx[0]);
        System.out.println("Theo R_x[0] = " + (sigma * sigma * (1 + b * b)));
        for (int k = 0; k < 8; k++) {
            System.out.println("PSD bin " + k + " = " + Sx[k]);
        }
    }
}

Chapter17_Lesson2.m

% Chapter17_Lesson2.m
% Autocorrelation, PSD, and White Noise Models (MATLAB / Simulink)

clear; clc; close all; rng(17);

N = 2048; sigma = 2.0; b = 0.6;
w = sigma * randn(N,1);             % white noise
x = filter([1 b], 1, w);            % x[n] = w[n] + b w[n-1]

M = 80;
RwFull = xcorr(w-mean(w), M, 'biased');
RxFull = xcorr(x-mean(x), M, 'biased');
Rw = RwFull(M+1:end);
Rx = RxFull(M+1:end);

Re = [Rx; Rx(end-1:-1:2)];
Sx = real(fft(Re));
omega = 2*pi*(0:length(Sx)-1)'/length(Sx);
SxTheory = sigma^2 * (1 + b^2 + 2*b*cos(omega));

fprintf('R_w(0) = %.6f\n', Rw(1));
fprintf('R_x(0) = %.6f\n', Rx(1));
fprintf('Theo R_x(0) = %.6f\n', sigma^2*(1+b^2));

figure; plot(0:M, Rw, 0:M, Rx); grid on;
xlabel('Lag m'); ylabel('Autocorrelation'); legend('White','Colored');

figure; half = 1:floor(length(Sx)/2);
plot(omega(half), Sx(half), omega(half), SxTheory(half), '--'); grid on;
xlabel('rad/sample'); ylabel('PSD'); legend('Estimated','Theory');

% Simulink setup:
% 1) Band-Limited White Noise block
% 2) Discrete Transfer Fcn with numerator [1 b], denominator [1]
% 3) Spectrum Analyzer + Scope

Chapter17_Lesson2.nb

(* Chapter17_Lesson2.nb *)
(* Wolfram Language script for autocorrelation and PSD *)

SeedRandom[17];
n = 2048; sigma = 2.0; b = 0.6;
w = RandomVariate[NormalDistribution[0, sigma], n];
x = Join[{w[[1]]}, Rest[w] + b Most[w]];

autocorrBiased[data_, mmax_] := Module[{z = N[data - Mean[data]], Nn = Length[data]},
  Table[Total[Take[z, Nn - m] Take[z, {1 + m, Nn}]]/Nn, {m, 0, mmax}]
];

rw = autocorrBiased[w, 80];
rx = autocorrBiased[x, 80];

re = Join[rx, Reverse[rx[[2 ;; -2]]]];
sx = Re[Fourier[re, FourierParameters -> {1, -1}]];
omega = Table[2 Pi (k - 1)/Length[sx], {k, Length[sx]}];
sxTheory = sigma^2 (1 + b^2 + 2 b Cos[omega]);

Print["R_w(0) = ", rw[[1]]];
Print["R_x(0) = ", rx[[1]]];
Print["Theo R_x(0) = ", sigma^2 (1 + b^2)];

ListLinePlot[{rw, rx}, PlotLegends -> {"White", "Colored"}, AxesLabel -> {"m", "R"}]
ListLinePlot[{Take[sx, Length[sx]/2], Take[sxTheory, Length[sxTheory]/2]},
  PlotLegends -> {"Estimated", "Theory"}, AxesLabel -> {"bin", "PSD"}]

Problem 1: Prove that for a real WSS process, \( R_{xx}[-m]=R_{xx}[m] \).

Solution: Use stationarity and commutativity of multiplication:

\[ R_{xx}[-m]=\mathbb{E}[x[n]x[n-m]] =\mathbb{E}[x[n-m]x[n]] =R_{xx}[m] \]

Problem 2: Show that \( \mathbb{E}[\hat{R}_{xx}^{(b)}[m]] = \frac{N-m}{N}R_{xx}[m] \).

Solution: Take expectation term-by-term:

\[ \mathbb{E}[\hat{R}_{xx}^{(b)}[m]] =\frac{1}{N}\sum_{n=0}^{N-1-m}\mathbb{E}[x[n]x[n+m]] =\frac{N-m}{N}R_{xx}[m] \]

Problem 3: For \( x[n]=w[n]+b\,w[n-1] \), compute \( R_{xx}[0] \) and \( R_{xx}[1] \).

Solution:

\[ R_{xx}[0]=(1+b^2)\sigma_w^2, \qquad R_{xx}[1]=b\sigma_w^2 \]

and by evenness \( R_{xx}[-1]=b\sigma_w^2 \).

Problem 4: If a zero-mean process has band-limited PSD \( S_{xx}(\omega)=S_0 \) on \( |\omega|\le\Omega_c \), compute the variance.

Solution:

\[ \operatorname{Var}(x(t)) =\frac{1}{2\pi}\int_{-\Omega_c}^{\Omega_c} S_0\,d\omega =\frac{S_0\Omega_c}{\pi} \]

Problem 5: Derive the PSD for \( R_{xx}[m]=\sigma_x^2 a^{|m|} \) with \( |a| < 1 \).

Solution:

\[ S_{xx}(e^{j\omega}) =\sigma_x^2\left(1+2\sum_{m=1}^{\infty}a^m\cos(m\omega)\right) =\sigma_x^2\frac{1-a^2}{1-2a\cos\omega+a^2} \]

Autocorrelation and PSD are the central second-order tools for stochastic system dynamics. The Wiener–Khinchin theorem links them exactly, white-noise models provide canonical inputs, and finite-record estimators make the theory operational for simulation and laboratory data.