Periodogram in MATLAB

Power Spectral Density Estimation Using the Periodogram

Step 1: Signal Representation

Let the signal be x[n], where:

• n = 0, 1, ..., N-1 (discrete-time indices),
• N is the total number of samples.

Step 2: Compute the Discrete-Time Fourier Transform (DTFT)

The DTFT of x[n] is:

X(f) = ∑ x[n] e-j2πfn

For practical computation, the Discrete Fourier Transform (DFT) is used:

X[k] = ∑ x[n] e-j(2π/N)kn, k = 0, 1, ..., N-1

• k represents discrete frequency bins,
• f_k = k/N * f_s, where f_s is the sampling frequency.

Step 3: Compute Power Spectral Density (PSD)

The periodogram estimates the PSD as:

S_x(f_k) = (1/N) |X[k]|²

• S_x(f_k) represents the power of the signal at frequency f_k.
• The factor 1/N normalizes the power by the signal length.

Step 4: Convert to Decibels (Optional)

For visualization, convert PSD to decibels (dB):

S_xdB(f_k) = 10 log₁₀(S_x(f_k))

Practical Notes

• Frequency Resolution: Depends on the signal duration T = N / f_s. Higher N gives finer frequency resolution.
• Windowing (Optional): Use a window function (e.g., Hamming, Hann) to reduce spectral leakage: x'[n] = x[n] * w[n].
• Frequency Range: PSD spans frequencies:
  • f_k = k/N * f_s for positive frequencies.
  • f_k = -(N-k)/N * f_s for negative frequencies.

MATLAB Code

clc;
clear;
close all;

% Define signal parameters
N = 256; % Number of samples
fs = 1000; % Sampling frequency in Hz
t = (0:N-1)/fs; % Time vector

% Generate a sample signal (sum of two sinusoids)
f1 = 50; % Frequency of first sinusoid in Hz
f2 = 120; % Frequency of second sinusoid in Hz
x = sin(2*pi*f1*t) + 0.5*sin(2*pi*f2*t) + randn(1, N)*0.1; % Signal with noise

% Compute DFT of the signal
X = fft(x, N);

% Calculate the periodogram
Pxx = (1/N) * abs(X).^2;

% Generate frequency vector
f = (0:N-1)*(fs/N);

% Keep only the positive frequencies (up to Nyquist frequency)
Pxx = Pxx(1:N/2+1);
f = f(1:N/2+1);

% Plot the periodogram
figure;
plot(f, 10*log10(Pxx), 'LineWidth', 1.5);
xlabel('Frequency (Hz)');
ylabel('Power/Frequency (dB/Hz)');
title('Periodogram');
grid on;

    

Output

Further Reading

1. Periodogram and Windowed Periododgram in details
2. Correlogram in MATLAB
3. Bartlett Method
4. Welch Method
5. Spectral Estimation Methods (Theory)
6. Comparison of Spectral Estimation Methods in MATLAB
7. How Windowing Affects Your Periodogram