Need an example for calculating power spectrum density

Question

Stefan 2014-7-28

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/143654-need-an-example-for-calculating-power-spectrum-density

评论： ROHAN JAIN 2020-5-4

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Greg Dionne 2014-7-29

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/143654-need-an-example-for-calculating-power-spectrum-density#answer_146800

编辑：Greg Dionne 2014-7-29

在 MATLAB Online 中打开

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 个评论
显示 3更早的评论隐藏 3更早的评论

Greg Dionne 2014-8-1

在 MATLAB Online 中打开

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

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

ROHAN JAIN 2020-5-4

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 2014-7-28

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/143654-need-an-example-for-calculating-power-spectrum-density#answer_146723

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