Need an example for calculating power spectrum density

27 views (last 30 days)
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.

Accepted Answer

Greg Dionne
Greg Dionne on 29 Jul 2014
Edited: Greg Dionne on 29 Jul 2014
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
  5 Comments
Greg Dionne
Greg Dionne on 1 Aug 2014
  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
ROHAN JAIN on 4 May 2020
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

Sign in to comment.

More Answers (1)

Greg Dionne
Greg Dionne on 28 Jul 2014
Hi Stefan,
There are numerous programs on MATLAB Central's File Exchange that can extract various features from ECG waveforms.
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
  1 Comment
Stefan
Stefan on 29 Jul 2014
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.

Sign in to comment.

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!