Power spectra

A couple of places calculate a power spectrum. The process is as follows.

Demean

ynorm = y-mean(y);

Calculate the power spectral density with “periodogram”

fs = 1/TR;
[pxx,f] = periodogram(ynorm,[],[],fs);

Remove meaningless low frequencies

% Remove meaningless low frequencies.
minf = 1/(t(end)/2);
indCutoff = find( f < minf , 1, 'last' );

pxx(1:indCutoff) = [];
f(1:indCutoff) = [];

Calculate total power

p = bandpower(pxx,f,'psd') returns the average power computed by integrating the power spectral density (PSD) estimate, pxx. The integral is approximated by the rectangle method. The input, f, is a vector of frequencies corresponding to the PSD estimates in pxx. The 'psd' option indicates that the input is a PSD estimate and not time series data.

For the input timeseries y: