Periodogram in MATLAB

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

Here:

• 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]|²

Where:

• 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 

   

  
  
  

   

   

   

   

  Copy the MATLAB Code above from here
  

  % The code is developed by SalimWireless.Com 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; web('https://www.salimwireless.com/search?q=salimwireless.com+power%20spectral%20density', '-browser');

   

  Further Reading

  1.