Welch's PSD estimate

10 次查看（过去 30 天）
Louise Wilson
Louise Wilson 2023-8-11
回答： Balaji 2023-8-23
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));

请先登录，再回答此问题。

回答（1 个）

Balaji
Balaji 2023-8-23
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

类别

在 Help Center 和 File Exchange 中查找有关 Parametric Spectral Estimation 的更多信息

标签

产品


版本

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by