How to do PSD measuremnet for the given signal

Hi, I am doing PSD (ms2) analysis for the frequency spectrum attached in the figure. I am able to find PSD but the results are not in PSD(ms2).

Here is the code:

figure
[pxx,f] = pwelch(rwave_decimate_1,[],[],[],fs/N);

%the mean
mean = 0;
for i=1:length(pxx)
  mean = mean + pxx(i);
end
M = mean/length(pxx); 

%the Varaince
var = 0;
for j=1:length(pxx)
    var = var + (pxx(j)-M)^2;
end
V = var/length(pxx); 

plot(f,pxx) 

It plots the PSD as in the attached figure. however, I want PSD in (ms2 or s2) on yaxis as in the figure 2, to make proper analysis on LF/HF region. Also, how can I highlight the LF and HF region with some shadded areas as in the figure2 ?

:- ms2 = mean-variance square

it would be appreciated if someone give solution.

Thanks in the advacnce.

[dpb]

M=mean(pxx);
V=var(pxx);

Don't use mean as a variable that aliases the builtin mean function in MATLAB.

[Sarfaraz Ahmed]

Thank you for correction.

I tried like the way you mentioned, however, it gives only 1 mean and 1 variance value.

What I want is that PSD(s^2/hz) values on y axis, not in by default unit. I mean i want something like in the figure2.

In my case, when I plot(f,pxx) -> it does not plot PSD in s^2/hz, instead it just plot something over frequency.

so, can you please assist how can we plot PSD in s^2/hz ?

Thank you.

Answer 1

[Arthi Sathyamurthi]

Hello Sarfaraz,

The units of the PSD estimate from the pwelch function are in squared magnitude units of the time series data per unit frequency. Thus, if the input data, rwave_decimate is in ms and the sample rate given as fs/N is in Hz, the PSD estimate which in your case pxx, ms^2/Hz.

Further, to highlight the LF and HF region with some shaded areas use the area function in MATLAB. You can also have a look at the answered question here for an example.

[Sarfaraz Ahmed]

Thank you Arthi for such response.

can I know please then how can I scale pxx in ms^2/Hz?

Let say, If I plot like that

figure

[pxx,f] = pwelch(rwave_decimate_1,[],[],[],fs/N);

plot(f,pxx)

would it be in ms^2/Hz?

Thanks !!

[Arthi Sathyamurthi]

If the input data, rwave_decimate is in ms and the sample rate given as fs/N is in Hz, the PSD estimate pxx will be in ms^2/Hz.

[Sarfaraz Ahmed]

Thank you.

Actually, its not the case. rwave_decimate is the signal after passing extrated r_waves signal to low pass filter. I don't think its already in the ms/Hz.