For loop for Welch Estimate

I have a problem in creating the for loop to calculate the Welch estimate without using the pwelch function. Can someone help me with this?

'We want to calculate the Welch estimate using a 100s sliding interval over the peridograms. To do this, calculate the first PSD estimate on the peridograms 1 to 10, the second estimate from 2 to 11, the third estimate from 3 to 12 etc.'

f_s = 1000;
M_welch = 10*f_s+1; % hann window length - odd 
w_welch = hann(M_welch); % window function to apply directly to signal
w_welch = w_welch/norm(w_welch); % normalise to unit energy
Lmax = floor(N/M_welch); % maximum number of peridograms 
L2 = 10; % number of peridograms in a 100s interval 
y = x;

Y = reshape(y(1:Lmax*M_welch),M_welch,Lmax)';  % split signal into L segments of length M 
% put each segment of signal into a coloumn and create a matrix 
Y = Y.* repmat(w_welch',Lmax,1);  % multiply by window function - element multiplication 

% phi_sum = zeros(1,2*M-1); % initialize sum for ACF
S_sum_welch = zeros(1,M_welch); % initialize sum for PSD 

f_Hz = (f_s/2)*(0:floor(M_welch/2))/(M_welch/2); % corresponding positive frequency vector (Hz)
Rp = 1;

deltaRatio_welch = zeros(1,Lmax-9);

for i = 1:Lmax-9 % first 100s
    for j = i:i+9 % loop from 1st periodogram to 10th periodogram
        S_sum_welch = S_sum_welch + 2*T*abs((fft(Y(j,:)))).^2; % update PSD % why is 2*T
    end
end

Answer 1

[Rohit Kulkarni]

Hi Yen,

I understand that you want to calculate Welch estimate without using the pwelch function.

Follow these steps to implement Welch estimate:

1. Divide the signal into overlapping segments.
2. Apply a window function to each segment.
3. Compute the Fast Fourier Transform (FFT) of each segment.
4. Average the squared magnitudes of the FFTs.
5. Scale the result by the sampling frequency and the window power.

You can refer to the modified code provided below, which implements all of these steps

Additionally, you can compare the results with the pwelch function and plot them on the same graph for easy comparison.

% set the sampling frequency and the window length
f_s = 1000;
M_welch = 10*f_s+1; % hann window length - odd 

% calculate the PSD using my code
w_welch = hann(M_welch); % window function to apply directly to signal
w_welch = w_welch/norm(w_welch); % normalise to unit energy
Lmax = floor(N/M_welch); % maximum number of peridograms 
L2 = 10; % number of peridograms in a 100s interval 
y = x;
Y = reshape(y(1:Lmax*M_welch),M_welch,Lmax)';  % split signal into L segments of length M 
% put each segment of signal into a coloumn and create a matrix 
Y = Y.* repmat(w_welch',Lmax,1);  % multiply by window function - element multiplication 
S_welch = zeros(Lmax-L2+1,M_welch); % initialize a matrix to store the PSD estimates for each window
T = 1/f_s; % sampling interval (s)

for i = 1:Lmax-L2+1 % loop over the sliding window length
    S_sum_welch = zeros(1,M_welch); % reset the sum for each window
    for j = i:i+L2-1 % loop over the segments in each window
        S_sum_welch = S_sum_welch + 2*T*abs((fft(Y(j,:)))).^2; % update PSD sum
    end
    S_sum_welch = S_sum_welch/L2; % divide by the number of segments to get the average PSD estimate
    S_welch(i,:) = S_sum_welch; % store the PSD estimate in the matrix
end

% calculate the PSD using the pwelch function of MATLAB
[S_pwelch,f_pwelch] = pwelch(x,w_welch,[],[],f_s);

% define the frequency vector for the PSD estimates
f_Hz = (f_s/2)*(0:floor(M_welch/2))/(M_welch/2); % corresponding positive frequency vector (Hz)

% plot the PSD estimates for comparison
figure;
plot(f_pwelch,S_pwelch,'b'); % plot the pwelch estimate in blue
hold on;
plot(f_Hz,S_welch(1,1:M_welch/2+1),'r'); % plot the first window estimate in red
plot(f_Hz,S_welch(end,1:M_welch/2+1),'g'); % plot the last window estimate in green
hold off;
xlabel('Frequency (Hz)');
ylabel('PSD');
legend('pwelch','first window','last window');
title('Comparison of PSD estimates');

Please refer to the documentation of 'pwelch' to know more about it: https://in.mathworks.com/help/signal/ref/pwelch.html

Hope this resolves your query!