Welch Method for Spectral Estimation in MATLAB

Welch’s Method

An approach for estimating the power spectral density by averaging periodograms to lower variance.

The signal is split into overlapping sections.
A window function is applied to each section.
The FFT is computed for each windowed section to obtain its power spectrum.
The power spectra are then averaged across all sections.

Copy the MATLAB Code from here

% The code is developed by SalimWireless.com

clc;
clear;
close all;

% Input Signal and Parameters

Fs = 1000;                 % Sampling frequency
t = 0:1/Fs:0.3;            % Time vector
x = cos(2*pi*200*t) + randn(size(t)); % Signal: 200 Hz cosine + noise

% Welch's Method Parameters

segmentLength = 256;       % Length of each segment
overlapFraction = 0.5;     % Fractional overlap (50%)
overlapSamples = floor(segmentLength * overlapFraction); % Overlap in samples
step = segmentLength - overlapSamples; % Step size
N = length(x);             % Length of the input signal

% Define Hann Window

hannWindow = 0.5 * (1 - cos(2 * pi * (0:segmentLength-1)' / (segmentLength - 1)));

% Segment the Signal into Overlapping Windows

numSegments = floor((N - segmentLength) / step) + 1;
segments = zeros(segmentLength, numSegments);

for i = 1:numSegments
    startIdx = (i - 1) * step + 1;
    endIdx = startIdx + segmentLength - 1;
    segments(:, i) = x(startIdx:endIdx);
end

% Initialize Power Spectrum Array

powerSum = zeros(floor(segmentLength/2) + 1, 1);

% Process Each Segment

for i = 1:numSegments

% Remove DC (mean)
    segment = segments(:, i) - mean(segments(:, i));

% Apply Hann Window
    windowedSegment = segment .* hannWindow;

% Compute FFT
    fftSegment = fft(windowedSegment);

% Compute One-Sided UNNORMALIZED Power Spectrum
    segmentPower = abs(fftSegment(1:floor(segmentLength/2) + 1)).^2;

% Accumulate
    powerSum = powerSum + segmentPower;
end

% Average Over Segments

powerSpectrum = powerSum / numSegments;

% One-Sided Adjustment (Exclude DC & Nyquist)

powerSpectrum(2:end-1) = 2 * powerSpectrum(2:end-1);

% Frequency Axis

frequencies = (0:floor(segmentLength/2)) * (Fs / segmentLength);

% Plot Unnormalized Spectrum

figure;
plot(frequencies, 10*log10(powerSpectrum));
xlabel('Frequency (Hz)');
ylabel('Power (dB)');
title('Welch Spectrum');
grid on;

web('https://www.salimwireless.com/search?q=salimwireless.com+welch%20spectral%20estimation', '-browser');

Search This Blog

Menu