How to perform a power spectral density (PSD) estimation of heart rate variability correctly?

I am currently trying to perform a PSD of heart rate variabiliy. I am using matlabs pwelch periodogram but apperently I am doing something wrong. To check the correctness of my results I am using a HRV analysis software called Kubios HRV. I hope someone knows what I am doing wrong.

The datainput is an Array of RR (had been filtered for artefacts before). I am detrending and resampling the data before I perform the pwelch function. I use the same Array in Kubios. There are no further artefact corrections applied in Kubios.

Here is my current matlab code:

inputArray = testData;
xValues = arrayfun(@(x) (sum(inputArray(1:x))/1000), 1:length(inputArray));

% detrending
detrendArray = detrend(inputArray);
trend = inputArray - detrendArray;

% resampling
Fs = 1/((sum(inputArray)/1000)/length(inputArray));
[resampled,xValuesNew] = resample(detrendArray,xValues,Fs);

% pwelch
pwelch(resampled,[],[],[],Fs);

In the attached images you can see the plotting result of my matlab pwech periodogram(blue line) as well as the correct periodogram from Kubios (grey, orange, blue). I am also attaching the test data that I used as a .csv file.

[Utsav Singh]

Hi, were you able to solve this probelm. I am currently facing a roadblock too as I am trying to do spectral analysis of HRV on similar data type as yours. Any help would be appreciated!