# What factors could affect the scale of PSD estimation?

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?

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:

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

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

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.

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

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

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

Bartlett's method is estimated as

$$$$\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}$$$$

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

$$$$\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}$$$$

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

$$$$\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}$$$$

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.

• 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! Commented Apr 11 at 22:02
• Thank you for the kind words! I appreciate that as well as the upvote and accept! Commented Apr 11 at 22:30