How to obtain PSD of time series data using pwelch()

Hi, I use pwelch() to calculate PSD of one set time-series data with 40,000 samples. The data is simulation result, and the time step is 0.01, so there are 100 samples per unit time, i.e. the sample rate fs = 100.

In the file of pwelch(), there are various commands, and I am not sure which one is more accurate. So I try different methods. It seems that method 1 and 2 have different value range of frequency, while method 3 looks very different with method 1 and 2. But all of these are not expected results, which should be something like:

or

Can anyone tell me is there any error in my code or the way I draw the PSD is incorrect? Thanks!!

The code is shown below:

load t;
load signal;
fs = 100;

% method 1
[pxx1, freq1] = pwelch(signal);
figure; plot(freq1, 10*log10(pxx1)); xlim([0 0.1]);

% method 2
[pxx2, freq2] = pwelch(signal,[],[],[],fs);
figure; plot(freq2, 10*log10(pxx2));  xlim([0 3]);

% method 3
[pxx3, freq3] = pwelch(signal,hanning(512),50,512,fs); 
figure; plot(freq3, 10*log10(pxx3)); xlim([0 3]);

Answer 1

[Mathieu NOE]

hello

if you want the finest frequency resolution, the fft length must be equal to your signal length. You could even virtually increase the resolution with zero padding the signal to further increase the fft length , but that comes with some drawbacks (sometimes)

Zero-Padding of FFTs (mechanicalvibration.com)

try this code : (it comes with additionnal features like filetring the signal , you can keep it or remove it)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

filename = 'signal.mat'; 
data = load(filename);
signal = (data.signal);
Fs = 100; % sampling rate is 100 Hz
dt = 1/Fs;

[samples,channels] = size(signal);

% time vector 
time = (0:samples-1)*dt;

%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 50;
if decim>1
    for ck = 1:channels
    newsignal(:,ck) = decimate(signal(:,ck),decim);
    Fs = Fs/decim;
    end
   signal = newsignal;
end
samples = length(signal);
time = (0:samples-1)*1/Fs;

%% band pass  filter section %%%%%%
    f_low = 0.05;
    f_high = 0.5;
    N_bpf = 2;           % filter order
    [b,a] = butter(N_bpf,2/Fs*[f_low f_high]);
    signal_filtered = filtfilt(b,a,signal);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = samples;    % for maximal resolution take NFFT = signal length.
OVERLAP = 0.75; % overlap will be used only if NFFT is shorter than signal length


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : time domain plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1),
subplot(211),plot(time,signal,'b');grid on
title(['Time plot  / Fs = ' num2str(Fs) ' Hz / raw data ']);
xlabel('Time (s)');ylabel('Amplitude');
subplot(212),plot(time,signal_filtered,'r');grid on
title(['Time plot  / Fs = ' num2str(Fs) ' Hz / filtered data ']);
xlabel('Time (s)');ylabel('Amplitude');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[freq, spectrum_raw] = myfft_peak(signal,Fs,NFFT,OVERLAP);
[freq, spectrum_filtered] = myfft_peak(signal_filtered,Fs,NFFT,OVERLAP);

figure(2),plot(freq,20*log10(spectrum_raw),'b',freq,20*log10(spectrum_filtered),'r');grid on
df = freq(2)-freq(1); % frequency resolution 
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz ']);
xlabel('Frequency (Hz)');ylabel('Amplitude (dB)');
legend('raw','filtered');

function  [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
% FFT peak spectrum of signal  (example sinus amplitude 1   = 0 dB after fft).
% Linear averaging
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer percentage of overlap % (between 0 and 0.95)

[samples,channels] = size(signal);

% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,channels);
    s_tmp((1:samples),:) = signal;
    signal = s_tmp;
    samples = nfft;
end

% window : hanning
window = hanning(nfft);
window = window(:);

%    compute fft with overlap 
 offset = fix((1-Overlap)*nfft);
 spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows

    % main loop
    fft_spectrum = 0;
    for i=1:spectnum
        start = (i-1)*offset;
        sw = signal((1+start):(start+nfft),:).*(window*ones(1,channels));
        fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft);     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end
    fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling

% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_spectrum = fft_spectrum(select,:);
freq_vector = (select - 1)*Fs/nfft;
end