4

I need to estimate the power spectral density of some signal and use the welch algorithm as provided by scipy.signal for this aim. Even after hours of research I couldn't find out, what exactly is the difference in the output when I set the kwarg 'scaling' to 'spectrum' in comparison to the default 'density' setting.

I thought the 'spectrum' might be the 'density' multiplied by frequency difference between two adjacent frequencies returned by welch, which would indeed turn a power density into a power. This is true only up to a factor of 1.5 which I really struggle to understand where it comes from. Any ideas?

import numpy as np
from scipy.signal import welch
dt = 1/4
tt = np.arange(0,1000,dt)
w = 1/5
data = 2 * np.sin(2*np.pi * w * tt) + np.random.normal(size = len(tt), scale = 1)

ff_welch, spectrum = welch(data, nperseg = 250, fs = 1/dt, scaling = 'spectrum')
ff_welch, density = welch(data, nperseg = 250, fs = 1/dt, scaling = 'density')
df_welch = (ff_welch[1]-ff_welch[0])
test =  spectrum / (density * df_welch)

the variable test will always be 1.5...

Olmo
  • 325
  • 4
  • 13
  • It looks like the scaling is a bit more complicated than just a factor of frequency resolution; I'm afraid I don't understand it (hence this isn't an answer!), but see https://holometer.fnal.gov/GH_FFT.pdf section 9, found via https://github.com/scipy/scipy/issues/8368. The 1.5 factor you found is a function of the window used; with window 'boxcar', the scale factor is 1. – Rory Yorke Nov 19 '21 at 08:25
  • I'm surprised it's 1.5, not sqrt(2), but that makes a lot of sense now... – Mad Physicist Jul 14 '22 at 03:45

1 Answers1

0

Density is Watts per Hertz whereas Spectrum is the Watts per bin. If you convert to dB, you just get the same plot shifted up or down.

In your case (at a sampling rate of 4 Hz with 126 bins in the output of the Welch function), you get ~0.0159 Hz per bin. That's less then a Hz so intuitively, it makes sense that the spectrum values are smaller. Typically, you could get from spectrum to density by doing 1/0.0159 = 63 and use 63 as a scaling factor. Windowing makes it trickier, but you can do it with the following:

  • Scaling factor for spectrum is: sum(ff_welch)^2
  • Scaling factor for density is: sum(ff_welch)**2 * fs/NFFT * sum(ff_welch) / sum(ff_welch**2)

Here's some code (apologies: I'm on an old Python 2 box) that shows you the ratio between the Welch output for spectrum and for density. If you do an FFT with your own windowing and shifting then use the scaling factors, you can get the same values as the welch function for spectrum. I've never gotten the density one exact so it varies a little bit past the third decimal place.

#!/usr/bin/env python
from __future__ import print_function

import numpy as np
from scipy.signal import welch
import matplotlib
import matplotlib.pyplot as plt

#Added '4.0'. This or 'float(1)/4' works. Python 2 doesn't know to make this a float.
dt = 1/4.0
tt = np.arange(0,1000,dt)
#Added '5.0'. This or 'float(1)/5' works. Python 2 doesn't know to make this a float.
w = 1/5.0
data = 2 * np.sin(2*np.pi * w * tt) + np.random.normal(size = len(tt), scale = 1)

ff_welch, spectrum = welch(data, nperseg = 250, fs = 1/dt, scaling = 'spectrum')
ff_welch, density = welch(data, nperseg = 250, fs = 1/dt, scaling = 'density')
#df_welch = (f1_welch[1]-f1_welch[0])
#test =  spectrum / (density * df_welch)
print ('density divided by spectrum is:', density/spectrum)

# Spectrum is for power per bin, density is for power per Hz.
# Scaling factor for spectrum is: 
#     sum(ff_welch)**2
# Scaling factor for density is:  
#     sum(ff_welch)**2 * fs/NFFT * sum(ff_welch) / sum(ff_welch**2) 
# If you divide density by spectrum, the sum(ff_welch)**2 cancels and you get: 
#    fs/NFFT * sum(ff_welch) / sum(ff_welch**2) 
fs = 1.0/dt
NFFT = len(density)
scaling_factor_density = sum(ff_welch)**2 * fs/NFFT * sum(ff_welch) / sum(ff_welch**2)
scaling_factor_spectrum = sum(ff_welch)**2
scaling_factor_comparison = 1 / (fs/NFFT * sum(ff_welch) / sum(ff_welch**2))
print ('scaling_factor_density:', scaling_factor_density)
print ('scaling_factor_spectrum:', scaling_factor_spectrum)
print ('spectrum scaling divided by density scaling:', scaling_factor_comparison)

freq_step_size = (fs/2.0) / NFFT
x_axis_Hz = np.arange(0,fs/2.0,freq_step_size)
plt.plot(x_axis_Hz, 10*np.log10(density), label='Density')
plt.plot(x_axis_Hz, 10*np.log10(spectrum), label='Spectrum')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Output from scipy.signal.welch in dB')
plt.legend()
plt.show()

Density vs Spectrum outputs converted to dB

scopil
  • 1