Welch's PSD estimate

I'm trying to work out why two different ways of doing a single-sided PSD analysis yields different results. As far as I can tell, these two methods below, are the same, with same FFT parameters, but I guess they can't be, because the results are ~ 40 dB different. Thanks in advance for any advice/help!

%Read audio file
clc;clear
load('xbit.mat');
Fs=96000;
xbit=detrend(xbit);

%FFT inputs
S=-163.287; %calibration
N=Fs; %segment length
r=0.5; %50% overlap
lcut=1; %low frequeny cut-off
hcut=48000;%high frequency cut-off


%% PWelch method
pxx=pwelch(xbit,N,[],Fs); %perform PSD analysis 
%signal, window/segment length, overlap (default = 50%), nfft size

pxx_dB=10*log10(pxx)-S; %convert to dB re 1µPa and calibrate

clearvars -except pxx_dB xbit Fs S N r lcut hcut


%% Manual method
% Divide signal into data segments 
xl = length(xbit);

xbit = single(xbit);                %reduce precision to single for speed
xgrid = buffer(xbit,N,ceil(N*r),'nodelay').';    
                                    %grid whose rows are each (overlapped) 
                                    %   segment for analysis
clear xbit
if xgrid(length(xgrid(:,1)),N) == 0 %remove final segment if not full
    xgrid = xgrid(1:length(xgrid(:,1))-1,:);
end

M = length(xgrid(:,1));             %total number of data segments

% Apply window function 
w = (0.54 - 0.46*cos(2*pi*(1:N)/N)); %Hamming window
alpha = 0.54; %scaling factor
xgrid = xgrid.*repmat(w/alpha,M,1); %multiply segments by window function

% Compute DFT 
X = abs(fft(xgrid.')).'; %calculate DFT of each data segment
clear xgrid

% Compute power spectrum 
P = (X./N).^2; %power spectrum = square of amplitude                                  
clear X

% Compute single-sided power spectrum 
Pss = 2*P(:,2:floor(N/2)+1);        %remove DC (0 Hz) component and 
                                    % frequencies above Nyquist frequency
                                    % Fs/2 (index of Fs/2 = N/2+1), divide
                                    % by noise power bandwidth

% Compute frequencies of DFT bins
f = floor(Fs/2)*linspace(1/(N/2),1,N/2);
                                    %calculate frequencies of DFT bins
flow = find(single(f) >= lcut,1,'first');   %low-frequency cut-off INDEX                                    
fhigh = find(single(f) <= hcut,1,'last');   %high-frequency cut-off INDEX
f = f(flow:fhigh);                  %frequency bins in user-defined range
nf = length(f);                     %number of frequency bins
Pss = Pss(:,flow:fhigh);            


% Compute noise power bandwidth and delta(f)
B = (1/N).*(sum((w/alpha).^2));     %noise power bandwidth
delf = Fs/N;                        %frequency bin width


% Convert to dB and apply calibration
a = 10*log10((1/(delf*B))*Pss./(1^2))-S;  

% Compute time vector
tint = (1-r)*N/Fs;                  %time interval in secs between segments
ttot = M*tint-tint;                 %total duration of file in seconds
t = 0:tint:ttot;                    %time vector in seconds   

% Construct output array
a = double(a);
A = zeros(M+1,nf+1);
A(2:M+1,2:nf+1) = a;
A(1,2:nf+1) = f; A(2:M+1,1) = t;

%Calculate average (frequency-bin wise)
A(2:end,2:end)=10.^(A(2:end,2:end)/10); %convert to linear space
mean_A=mean(A(2:end,2:end));
mean_A=10*log10(mean_A);
mean_A=[0 mean_A];
A=[A(1,:); mean_A];


%% Compare the results of the two methods visually
plot(pxx_dB);
hold on
plot(A(2,2:end));

Answer 1

[Balaji]

Hi Louise,

As per my understanding, you are facing an issue while implementing the single-sided PSD.

The “pwelch” function takes in the following arguments:

pxx = pwelch(x,window,noverlap,nfft)

Assume the default Fs to be 1Hz unless specified specifically by:

[pxx,f] = pwelch(___,fs) 

You will be able to replicate the results you have obtained through the second method that you have used by modifying the 14th line of your code as follows:

pxx=pwelch(xbit,N,[],N,Fs); 

For more information on “pwelch” function please refer to the documentation provided below.

https://in.mathworks.com/help/signal/ref/pwelch.html