Need an example for calculating power spectrum density

Hello,

I want to plot a Power Spectrum Density(having units s^2/Hz)plot against frequency(Hz) as shown in this link PSD and want to calculate the variables PeakFreq,VLFpower,LFpower,UFpower,VLF,LF,UF as shown in link.

Can someone please give a demo with some example values of

1)how to plot the PSD as shown in the link

2)how to calculate those variables.

Thankyou.

Answer 1

[Greg Dionne]

Hi Stefan,

Thanks for the data set. Here's an idea to get you started.

% RRintervals is a list of the intervals between successive R waves

% approximate the time relative to the first interval
t = cumsum(RRintervals);

% make a uniform grid inside the full interval 
% I arbitrarily picked a grid at 1 second intervals... you may find something else more useful.
tGrid = 1:55;

% interpolate to get the beats per second on a uniform grid
bps = interp1(t,RRintervals,tGrid,'spline');

% compute the PSD
% units of Pxx are squared seconds/Hz.
% sample rate is 1 Hz.
[Pxx,F] = periodogram(bps,[],numel(bps),1);

% Zero out the DC bin so it fits nicely on screen
Pxx(1) = 0;

% make a pretty plot
plot(F,Pxx);
xlabel('Hertz');
ylabel('s^2 / Hz');

% compute the power in the various bands...
% note that I removed the DC bin earlier, so VLF is somewhat suspect...
vlf = bandpower(Pxx,F,[0 0.04],'psd') % units of sec^2
lf = bandpower(Pxx,F,[0.04 0.15],'psd') % units of sec^2
hf = bandpower(Pxx,F,[0.15 0.4],'psd') % units of sec^2

totPower = bandpower(Pxx,F,'psd') % units of sec^2
% you can then take the ratio of lf, hf, etc. to totPower * 100 to get the percentages etc.

Hope this helps.

-G

[Stefan]

Hi Greg Dionne,

Thanks for providing such a good implementation which is the needed. But I have got some doubts as

1)Why interp1 is used(or what it does) and why the method ‘spline’ is used

2)Also can please explain why the sample rate of 1Hz is selected

Can you help me in understanding these.

[Greg Dionne]

Hi Stefan,

Good questions.

I used interp1 in conjunction with cumsum as an expedient way to compute beats per second on a uniform grid so I could pass it periodogram.

Spline has a smoother way of connecting points so it introduces less high frequency noise. It probably doesn't matter too much, since you're looking at low frequencies anyway.

As for the 1 Hz sample rate, I saw that you had about a minute's worth of data, and just around 60 samples within it. Was "easy" to think in terms of seconds. Don't feel like you have to use that number. You can always increase the sample rate a bit more if you're nervous about aliasing effects.

Hope this helps.

-G

[Stefan]

Thanks Greg, Can you please help me otu wiht these as well.

1)how to get the units in msec^2

2)How to get the frequency corresponding to maximum power in a each particular band(for example how to get the peak freq in the VLF band of 0 to 0.04Hz)

Thanks.

[Greg Dionne]

1. To convert X s^2 to X ms^2, multiply by 1e6.
2. You could try something like:

iRange = find(0 <= F & F <= 0.04);
[pMax, iMax] = max(Pxx(iRange))
fMax = F(iRange(iMax))

[ROHAN JAIN]

Hi, I have a general doubt regarding the previous question. For very low frequencies (below 2Hz) the power of the psd is by default very high because of the presence of 1/f inherent noise, and as a result of which many interpretations go wrong. So, the high amplitudes or power values at such low frequencies is not due to the signal but noise actually.

Do you have any suggestions for the removal of such noise from EEG signal in an effective way ?

Many Thanks,

Best Wishes,

RJ

Answer 2

[Greg Dionne]

Hi Stefan,

There are numerous programs on MATLAB Central's File Exchange that can extract various features from ECG waveforms.

http://www.mathworks.com/matlabcentral/fileexchange/index?utf8=%E2%9C%93&term=ECG

If you just need something simple and have a recent copy of the Signal Processing Toolbox, you can use PERIODOGRAM without output arguments to plot a PSD. You can also use BANDPOWER to obtain the bandpower between two frequencies of a uniformly sampled signal or a PSD.

Hope this helps.

-Greg

[Stefan]

Hi,

I tried both PSD and periodogram for the RRinterval(as attached)I have calculated.

Can I know

1) how to get the plot with PSD units as 's^2/Hz'

2) how to calculate the variables shown in the below image with respective units from PSD.

Thanks.