I'm trying to calculate the power spectral density (PSD) to analyze signal power. To do so, I wrote a MATLAB code with the great help of others and below is my code.

function [psd, freq] = mypsd(data, fs)
    N = length(data);
    f_delta = fs/floor(N/2)/2; % Modify N to even number
    freq = (0:f_delta:fs/2)';

    xdft = fft(data);
    xdft = abs(xdft/N);
    xft = xdft(1:floor(N/2)+1);
    xft(2:end-1) = 2*xft(2:end-1);
    psdx = xft.^2 / (2*f_delta);
    psd = 10*log10(psdx);

To confirm the code, I compared the result with the one from the pwelch function

my_result = mypsd(data, fs);
mat_result = pwelch(data, rectwin(length(data)), 0, length(data), fs);

In this case above, two results were same.

After that, I tried to segment the data with overlapping 50% of each segment like below, where nfft is 256, n_overlap is nfft/2 = 128.

my_psd_total = zeros(nfft/2+1, 1);

for iter = 1:num_segment
    data_win = data((iter-1)*(nfft - n_overlap)+1:(iter-1)*(nfft - n_overlap)+nfft);
    [psd_win, freq_win] = mypsd(data_win.*hann(nfft), fs);
    my_psd_total = my_psd_total + psd_win;

my_psd_total = my_psd_total / num_segment

However, the result is different with the result of pwelch(data, hann(nfft), n_overlap, nfft, fs)

As I know, elementwise multiplying of window with segment in time domain and average the result is the proper method of using window functions like hanning or hamming et al.

Is there anyone can explain why my result is different with built-in pwelch function?


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.