The correct sampling frequency is to be provided as argumentfs 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()Welch’s method [1] 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
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