2
$\begingroup$

I have written this function in R to compute the power spectrum density (PSD) of a signal:

my_psd <- function(x, sampling_rate, han = TRUE){
    N <- length(x)
    if (han){
        hann <- gsignal::hann(N) # Hanning window
        x <- x * hann
    }
    ft <- abs(fft(x)) # Take absolute value of fft.
    ft <- ft[1:(N/2 + 1)] # Make one sided
    freq <- c(0:(length(ft)-1))*sampling_rate/length(x) # One sided fft
    
    power <- ( 1/( N * sampling_rate )) * ft^2
    power <- 2 * power[2:( length(ft) - 1 )]
    power <- append(power, ft[1], after=0)
    power <- append(power, ft[length(ft)])
    return(list(spec = power, freq = freq))
}

As far as my understanding goes, it is correct. I computes the absolute value of the FFT of a signal (after optionally applying a Hanning window), which is ft. Then disregards half of the values in ft (one-sided). It computes the frequencies in freq. To compute the power spectrum density, it calculates

$$\frac{1}{N \cdot fs} |H(f)|^2$$

where $H(f)$ is the fft, where $N$ is the length of the signal. It then multiplies by two all values in the resulting vector, except the first ($0 Hz$) and the last (the one corresponding to the Nyquist frequency). It returns a named list where spec holds the spectrum and freq the frequencies.

I simulated a simple signal as follows:

# Set seed for reproducibility
set.seed(123)

# Parameters
sampling_rate <- 500 # Sampling rate (samples per second)
duration <- 30           # Duration of the signal in seconds
frequencies <- c(10, 30, 50)  # Frequencies of the sine waves (in Hz)
amplitudes <- c(5, 10, 15)   # Amplitudes of the sine waves

# Generate time vector
t <- seq(0, duration, by = 1/sampling_rate)
n <- length(t)

# Initialize signal as zeros
signal <- rep(0, n)

# Create a signal as a mixture of sine waves
for (i in seq_along(frequencies)) {
  frequency <- frequencies[i]
  amplitude <- amplitudes[i]
  
  # Generate the sine wave component
  sine_wave <- amplitude * sin(2 * pi * frequency * t)
  
  # Add the sine wave component to the signal
  signal <- signal + sine_wave
}

The signal is simply the sum of three sine waves with amplitudes $5, 10, 15$ and frequencies $10, 30, 50$ respectively.

I am comparing the difference between the power spectrum with and without the Hanning window. I am finding that the Hanning window exacerbates the most negative powers. See the image below, where the left plot is with Hanning window and the second plot without Hanning window:

With Hanning left / Without Hanning right

(Note. The spectrums in the plot were converted to decibels) Why is this huge difference occurring? The most negative values with a Hanning window reach $-200$, whereas without a Hanning window power barely reaches $-100$. This is a difference of a factor of $2$, which I don't think is negligible.

I am quite new to FFT and I'm doing this to learn, so excuse me if I'm missing something obvious.

$\endgroup$

1 Answer 1

2
$\begingroup$

When analyzing pure tones, you’d want to calculate the power spectrum, not the PSD. I’ll answer your question first, and give you the correct power spectrum computation second.

Denote by $|X|^2$ the Squared Magnitude Spectrum and $w$ the window applied.

  • The correctly scaled PSD is: $$\frac{2|X|^2}{f_s \times S_2} \quad \texttt{with} \quad S_2 = \sum_{i = 0}^{N - 1} w_i^2$$

For a rectangular window (i.e. no windowing applied), $S_2 = N$.

For a Hann window, $S_2 \approx 0.37N$, hence the scaling issue you’re seeing.

  • The Power Spectrum is: $$\frac{2|X|^2}{(S_1)^2} \quad \texttt{with} \quad S_1 = \sum_{i = 0}^{N - 1} w_i$$

For more detail, please refer to this

$\endgroup$
7
  • $\begingroup$ I suppose I am a bit confused about the difference between power spectrum and PSD, and why the latter should not be used to "pure tones" (by which I gather you mean the absolute resulting values). $\endgroup$
    – lafinur
    Commented Apr 6 at 23:10
  • $\begingroup$ Btw is $X$ the one-sided spectrum (and hence the factor of $2$ in the scaled PSD)? $\endgroup$
    – lafinur
    Commented Apr 6 at 23:16
  • $\begingroup$ Did you read the answer I linked you to? It has to do with the distribution of power across frequency mainly, but that’s a whole different question than the one you’ve asked. Maybe this can help you make the distinction. $\endgroup$
    – Jdip
    Commented Apr 6 at 23:19
  • $\begingroup$ Yes $X$ is one-sided ! $\endgroup$
    – Jdip
    Commented Apr 6 at 23:20
  • $\begingroup$ Thanks, I'm giving a careful read to the answer you linked. Excellent material. I really appreciate the help! $\endgroup$
    – lafinur
    Commented Apr 6 at 23:20

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.