MATLAB Code
% The code is developed by SalimWireless.comclc;
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)));
windowEnergy = sum(hannWindow.^2); % Normalization factor
% Segment the Signal into Overlapping Windows
numSegments = floor((N - overlapSamples) / step); % Number of segments
segments = zeros(segmentLength, numSegments);
for i = 1:numSegments
startIdx = (i - 1) * step + 1;
endIdx = startIdx + segmentLength - 1;
segments(:, i) = x(startIdx:endIdx);
end
% Initialize PSD Array
psdSum = zeros(floor(segmentLength/2) + 1, 1);
% Process Each Segment
for i = 1:numSegments
% Apply Hann Window to the Segment
windowedSegment = segments(:, i) .* hannWindow;
% Compute FFT of Windowed Segment
fftSegment = fft(windowedSegment);
% Compute One-Sided PSD
segmentPSD = (1 / (Fs * segmentLength * windowEnergy)) * abs(fftSegment(1:floor(segmentLength/2) + 1)).^2;
% Accumulate PSD
psdSum = psdSum + segmentPSD;
end
% Average PSD Over All Segments
psd = psdSum / numSegments;
% Adjust for One-Sided Spectrum (Excluding DC and Nyquist)
psd(2:end-1) = 2 * psd(2:end-1);
% Generate Frequency Axis
frequencies = (0:floor(segmentLength/2)) * (Fs / segmentLength);
% Plot the PSD
figure;
plot(frequencies, 10*log10(psd)); % Plot PSD in dB/Hz
xlabel('Frequency (Hz)');
ylabel('Power Spectral Density (dB/Hz)');
title('Power Spectral Density (Welch’s Method)');
grid on;
startIdx = (i - 1) * step + 1;
endIdx = startIdx + segmentLength - 1;
segments(:, i) = x(startIdx:endIdx);
end
% Initialize PSD Array
psdSum = zeros(floor(segmentLength/2) + 1, 1);
% Process Each Segment
for i = 1:numSegments
% Apply Hann Window to the Segment
windowedSegment = segments(:, i) .* hannWindow;
% Compute FFT of Windowed Segment
fftSegment = fft(windowedSegment);
% Compute One-Sided PSD
segmentPSD = (1 / (Fs * segmentLength * windowEnergy)) * abs(fftSegment(1:floor(segmentLength/2) + 1)).^2;
% Accumulate PSD
psdSum = psdSum + segmentPSD;
end
% Average PSD Over All Segments
psd = psdSum / numSegments;
% Adjust for One-Sided Spectrum (Excluding DC and Nyquist)
psd(2:end-1) = 2 * psd(2:end-1);
% Generate Frequency Axis
frequencies = (0:floor(segmentLength/2)) * (Fs / segmentLength);
% Plot the PSD
figure;
plot(frequencies, 10*log10(psd)); % Plot PSD in dB/Hz
xlabel('Frequency (Hz)');
ylabel('Power Spectral Density (dB/Hz)');
title('Power Spectral Density (Welch’s Method)');
grid on;
Output
