 # prefactors computing psd of a signal with numpy.fft vs. scipy.signal.welch

• Last Update :
• Techknowledgy :

The correct sampling frequency is to be provided as argument`fs` to retreive the correct frequencies and an accurate power spectral density. To recover a power spectrum similar to that computed using `np.multiply(u_fft, np.conj(u_fft))`, the length of the fft frame and the applied window must respectively be provided as the length of the frame and `boxcar` (equivalent to no window at all). The fact that `scipy.signal.welch` applies the correct scaling can be checked by testing a sine wave:

```import numpy as np
import scipy.signal

import matplotlib.pyplot as plt

def abs2(x):
return x.real ** 2 + x.imag ** 2

if __name__ == '__main__':
framelength = 1.0
N = 1000
x = np.linspace(0, framelength, N, endpoint = False)
y = np.sin(44 * 2 * np.pi * x)
#y = y - np.mean(y)
ffty = np.fft.fft(y)
#power spectrum, after real2complex transfrom(factor)
scale = 2.0 / (len(y) * len(y))
power = scale * abs2(ffty)
freq = np.fft.fftfreq(len(y), framelength / len(y))

# power spectrum, via scipy welch.
'boxcar'
means no window, nperseg = len(y) so that fft computed on the whole signal.
freq2, power2 = scipy.signal.welch(y, fs = len(y) / framelength, window = 'boxcar', nperseg = len(y), scaling = 'spectrum', axis = -1, average = 'mean')

for i in range(len(freq2)):
print i, freq2[i], power2[i], freq[i], power[i]
print np.sum(power2)

plt.figure()
plt.plot(freq[0: len(y) / 2 + 1], power[0: len(y) / 2 + 1], label = 'np.fft.fft()')
plt.plot(freq2, power2, label = 'scipy.signal.welch()')
plt.legend()
plt.xlim(0, np.max(freq[0: len(y) / 2 + 1]))

plt.show()```

Suggestion : 2

Welch’s method  computes an estimate of the power spectral density by dividing the data into overlapping segments, computing a modified periodogram for each segment and averaging the periodograms.,P. Welch, “The use of the fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms”, IEEE Trans. Audio Electroacoust. vol. 15, pp. 70-73, 1967.,Estimate power spectral density using Welch’s method.,If we now introduce a discontinuity in the signal, by increasing the amplitude of a small portion of the signal by 50, we can see the corruption of the mean average power spectral density, but using a median average better estimates the normal behaviour.

```>>> from scipy
import signal
>>>
import matplotlib.pyplot as plt >>>
rng = np.random.default_rng()```
```>>> fs = 10e3 >>>
N = 1e5 >>>
amp = 2 * np.sqrt(2) >>>
freq = 1234.0 >>>
noise_power = 0.001 * fs / 2 >>>
time = np.arange(N) / fs >>>
x = amp * np.sin(2 * np.pi * freq * time) >>>
x += rng.normal(scale = np.sqrt(noise_power), size = time.shape)```
```>>> f, Pxx_den = signal.welch(x, fs, nperseg = 1024) >>>
plt.semilogy(f, Pxx_den) >>>
plt.ylim([0.5e-3, 1]) >>>
plt.xlabel('frequency [Hz]') >>>
plt.ylabel('PSD [V**2/Hz]') >>>
plt.show()```
```>>> np.mean(Pxx_den[256: ])
0.0009924865443739191```
```>>> f, Pxx_spec = signal.welch(x, fs, 'flattop', 1024, scaling = 'spectrum') >>>
plt.figure() >>>
plt.semilogy(f, np.sqrt(Pxx_spec)) >>>
plt.xlabel('frequency [Hz]') >>>
plt.ylabel('Linear spectrum [V RMS]') >>>
plt.show()```
```>>> np.sqrt(Pxx_spec.max())
2.0077340678640727```

Suggestion : 3

If True, divide by log2(psd.size) to normalize the spectral entropy between 0 and 1. Otherwise, return the spectral entropy in bit.,'fft' : Fourier Transform (scipy.signal.periodogram()),Inouye, T. et al. (1991). Quantification of EEG irregularity by use of the entropy of the power spectrum. Electroencephalography and clinical neurophysiology, 79(3), 204-210.,'welch' : Welch periodogram (scipy.signal.welch())

```>>>
import numpy as np
>>>
import entropy as ent >>>
sf, f, dur = 100, 1, 4 >>>
N = sf * dur # Total number of discrete samples >>>
t = np.arange(N) / sf # Time vector >>>
x = np.sin(2 * np.pi * f * t) >>>
np.round(ent.spectral_entropy(x, sf, method = 'fft'), 2)
0.0```
```>>> np.random.seed(42) >>>
x = np.random.rand(3000) >>>
ent.spectral_entropy(x, sf = 100, method = 'welch')
6.980045662371389```
```>>> ent.spectral_entropy(x, sf = 100, method = 'welch', normalize = True)
0.9955526198316071```
```>>> np.random.seed(42) >>>
x = np.random.normal(size = (4, 3000)) >>>
np.round(ent.spectral_entropy(x, sf = 100, normalize = True), 4)
array([0.9464, 0.9428, 0.9431, 0.9417])```
```>>>
import stochastic.processes.noise as sn >>>
rng = np.random.default_rng(seed = 42) >>>
x = sn.FractionalGaussianNoise(hurst = 0.5, rng = rng).sample(10000) >>>
print(f "{ent.spectral_entropy(x, sf=100, normalize=True):.4f}")
0.9505```
```>>> rng = np.random.default_rng(seed = 42) >>>
x = sn.FractionalGaussianNoise(hurst = 0.9, rng = rng).sample(10000) >>>
print(f "{ent.spectral_entropy(x, sf=100, normalize=True):.4f}")
0.8477```