1
$\begingroup$

I was required to program a function in R to estimate power spectral density of EEG data. The function I created is, I think, rather standard, and estimates the PSD of a signal by computing

$$ \frac{2|X[k]|^2}{S \times f_s} $$

where $X[k]$ is the $k$th element in the DFT, $S = \sum_{i=0}^{N-1} w_i^2$ the sum of the squared window function, and $f_s$ the sampling rate. The code, specifically, is:

my_psd <- function(x, sampling_rate, han = TRUE){
    N <- length(x)
    if (han){
        hann <- gsignal::hann(N) # Hanning window
        x <- x * hann
        normalization <- 2/(sum(hann^2) * sampling_rate)
    }else{
        normalization <- 2/(N * sampling_rate)
    }
    ft <- abs(fft(x)) # Take absolute value to remove complex numbers and normalize with 1/N.
    ft <- ft[1:(N/2 + 1)] # Make one sided
    freq <- c(0:(length(ft)-1))*sampling_rate/length(x) # One sided fft
    
    power <- ft^2 * normalization
    return(list(spec = power, freq = freq))
}

My boss gave me an EDF with EEG data and the PSD estimation other researchers obtained in the past with this data. I will call the PSD estimation she sent me $\hat{P}(f)$.

She asked me to run my PSD function and compare my results with $\hat{P}(f)$. But the scale of my results are far from the scale of the PSD she sent me. When $f$ is small (say, in the range $[0, 0.4]$), $\hat{P}(f)$ rounds $\approx 100$. But my function gives me values of $\approx 2000$ in the same range!

$\hat{P}(f)$, the PSD I must compare my own PSD with, comes from a very reputable laboratory and I am pretty sure it must be well computed. And I cannot explain the difference in the scale of my results...

I have tried the obvious things; e.g. taking the square root of my results, and filtering the data in a way which I have reasons to believe the lab that computed $\hat{P}$ may have filtered it. But the difference in scale remains.

What could be behind the scale difference? Are there known factors that affect the scale of the PSD of a signal?

EDIT: I add the plots that were asked for in the comments.

1. For each $5$ second epoch in the EEG, the my_psd function was used to compute the spectrum; then the spectrum was averaged across epochs. This is the same procedure used for $\hat{P}$. My result looks like this:

enter image description here

while $\hat{P}$ looks like this:

enter image description here

2. Using Welch's method. Using Welch's method, with window length $5$ seconds and overlap $0.5$, the PSD looks like this:

enter image description here

Welch's method gives less of a mismatch in terms of scale. But still, power reaches $400$, while in $\hat{P}$ it barely reaches a hundred.

$\endgroup$
7
  • $\begingroup$ Can you show us a plot with your computed PSD and theirs super-imposed? $\endgroup$
    – Jdip
    Commented Apr 8 at 23:52
  • $\begingroup$ The classic way of estimating PSDs is using welch's method $\endgroup$
    – Jdip
    Commented Apr 8 at 23:56
  • $\begingroup$ What is your sample rate? $\endgroup$
    – Baddioes
    Commented Apr 9 at 2:03
  • $\begingroup$ Sample rate is 500Hz. I have used Welch's method and the discrepancies remain. Will update a plot. $\endgroup$
    – lafinur
    Commented Apr 10 at 13:57
  • $\begingroup$ I've uploaded the plots. I didn't superimpose them because the difference in scales made $\hat{P}$ almost impossible to appreciate. $\endgroup$
    – lafinur
    Commented Apr 10 at 14:28

1 Answer 1

2
$\begingroup$

I think the issue here is scaling. Possibly contrary to popular belief, there is not a singular method of scaling that works generically for all PSD estimates, as each PSD type estimates the spectrum differently. Additionally, the scaling you provide is specifically for referencing back to true power as opposed to a sort of "digital" power. Often times, this is simply ignored as it is merely a scale factor difference. I have a feeling that $\hat{P}$ may be in this "digital" power.

Ignoring the $f_{s}$, the scale factor you provided relates specifically to the unmodified periodogram. In the periodogram, the power in the temporal window is $N$. When you start doing averaging, the power in the temporal window goes down as the data length decreases, however, you are not simply adding the windowed spectral estimates, you are averaging them. Therefore, you also have to scale by the number of windows.

Mathematically explaining this, if the "else" conditional is met within your my_psd function, this is a method known as Bartlett's method. Specifically, there is no overlap, and technically no windowing within Bartlett's method, as both of these are introduced in Welch's method. Bartlett's method is a direct averaged periodogram, so if we have

\begin{equation} \hat{\phi}_{j} = \frac{1}{N}|\sum_{k=0}^{N-1}y_{j}[k]e^{-j\omega k}|^{2}\end{equation}

Bartlett's method is estimated as

\begin{equation} \hat{\phi}_{B} = \frac{1}{M}\sum_{j=0}^{M-1}\hat{\phi}_{j} = \frac{1}{MN}\sum_{j=0}^{M-1}|\sum_{k=0}^{N-1}y_{j}[k]e^{-j\omega k}|^{2}\end{equation}

So, you are normalizing by $\frac{1}{\text{(number of samples per window)}\text{(number of windows being averaged)}}$. This will likely get you much closer to $\hat{P}$. If your data is real valued, you may need to keep the 2 in the numerator, see my answer here for more details on that.

If you are going to window the data, then you have to scale by $\frac{1}{MNP}$ where $P$ is the power in the temporal window, ie $\frac{1}{N}\sum_{k=0}^{N-1}|w[k]|^{2}$. Plugging in and simplifying we would get

\begin{equation} \hat{\phi}_{O} = \frac{1}{M\sum_{k=0}^{N-1}|w[k]|^{2}}\sum_{j=0}^{M-1}|\sum_{k=0}^{N-1}w[k]y_{j}[k]e^{-j\omega k}|^{2} \end{equation}

Extending this to Welch's method, if you use 50% overlap you would then substitute $L = 2M$ and get something like

\begin{equation} \hat{\phi}_{W} = \frac{1}{L\sum_{k=0}^{N-1}|w[k]|^{2}}\sum_{j=0}^{L-1}|\sum_{k=0}^{N-1}w[k]y_{j}[k]e^{-j\omega k}|^{2} \end{equation}

Hopefully this helps get you closer to proper scaling. If you have any more questions or concerns, don't hesitate to comment! I also referenced 1 for many of these formulas.

1 Spectral Analysis of Signals by Stoica and Moses.

$\endgroup$
2
  • 1
    $\begingroup$ Superb answer, you hit the nail in the head. The results match now. If we were in MathStackExchange I'd gift you reputation, but I have none on this site. Again, really detailed answer and I'm in fact impressed on how accurately you diagnosed the issue. Thanks! $\endgroup$
    – lafinur
    Commented Apr 11 at 22:02
  • $\begingroup$ Thank you for the kind words! I appreciate that as well as the upvote and accept! $\endgroup$
    – Baddioes
    Commented Apr 11 at 22:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.