matplotlib.mlab.psd(...) और scipy.signal.welch(...) सिग्नल के पावर स्पेक्ट्रल घनत्व (PSD) का अनुमान लगाने के लिए वेल्च की औसत अवधि अवधि विधि लागू करते हैं। मान लें कि बराबर पैरामीटर प्रत्येक फ़ंक्शन को पास किए जाते हैं, लौटाया गया PSD बराबर होना चाहिए।matplotlib.mlab.psd और scipy.signal.welch से पावर स्पेक्ट्रल घनत्व का अनुमान क्यों भिन्न होता है जब प्रति विंडो बिंदुओं की संख्या भी होती है?
हालांकि, एक सरल sinusoidal परीक्षण समारोह का उपयोग कर, वहाँ दो तरीकों के बीच व्यवस्थित मतभेद हैं जब खिड़की प्रति इस्तेमाल किया अंकों की संख्या (scipy.signal.welch के लिए nperseg कीवर्ड; mlab.psd के लिए NFFT कीवर्ड) भी है, जिन्हें आप नीचे खिड़की प्रति 64 अंक के मामले के लिए
शीर्ष साजिश, जबकि नीचे की साजिश उनके रिश्तेदार त्रुटि (उनके निरपेक्ष averag से विभाजित दो PSDs का यानी पूर्ण अंतर को प्रदर्शित करता है, PSD दोनों विधियों के माध्यम से गणना की चलता ई)। एक बड़ी सापेक्ष त्रुटि दो तरीकों के बीच बड़ी असहमति का संकेत है।
दो काम करता है, जब खिड़की प्रति इस्तेमाल किया अंकों की संख्या अजीब है ज्यादा बेहतर समझौता है के रूप में एक ही संकेत के लिए नीचे दिखाया गया है, लेकिन खिड़की
जोड़ना प्रति 65 अंक के साथ कार्रवाई सिग्नल के लिए अन्य विशेषताएं (जैसे शोर) कुछ हद तक इस प्रभाव को कम करती है, लेकिन यह अभी भी मौजूद है, दो विधियों के बीच ~ 10% की सापेक्ष त्रुटियों के साथ जब प्रति विंडो की संख्या भी उपयोग की जाती है। (मुझे ऊपर की प्रति विंडो 65 अंक के साथ गणना की गई PSDs के लिए उच्चतम आवृत्ति पर सापेक्ष त्रुटि का एहसास है। हालांकि, यह सिग्नल की निक्विस्ट आवृत्ति पर है, और मैं इस तरह की उच्च आवृत्तियों पर सुविधाओं के बारे में चिंतित नहीं हूं। सिग्नल की बैंडविड्थ के बहुमत में बड़ी और व्यवस्थित सापेक्ष त्रुटि के बारे में अधिक चिंतित हूं, जब प्रति विंडो बिंदुओं की संख्या भी होती है)।
क्या इस विसंगति का कोई कारण है? सटीकता के मामले में एक विधि दूसरे के लिए बेहतर है? मैं scipy संस्करण 0.16.0 और matplotlib संस्करण 1.4.3 का उपयोग कर रहा हूँ।
उपरोक्त भूखंडों को उत्पन्न करने के लिए उपयोग किया गया कोड नीचे शामिल है। शुद्ध साइनसॉइड सिग्नल के लिए, A_noise = 0 सेट करें; एक शोर संकेत के लिए, A_noise को एक सीमित मूल्य पर सेट करें।
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import welch
from matplotlib import mlab
# Sampling information
Fs = 200.
t0 = 0
tf = 10
t = np.arange(t0, tf, 1./Fs)
# Pure sinusoidal signal
f0 = 27.
y = np.cos(2 * np.pi * f0 * t)
# Add in some white noise
A_noise = 1e-3
y += (A_noise * np.random.randn(len(y)))
nperseg = 64 # even
# nperseg = 65 # odd
if nperseg > len(y):
raise ValueError('nperseg > len(y); preventing unintentional zero padding')
# Compute PSD with `scipy.signal.welch`
f_welch, S_welch = welch(
y, fs=Fs, nperseg=nperseg, noverlap=(nperseg // 2),
detrend=None, scaling='density', window='hanning')
# Compute PSD with `matplotlib.mlab.psd`, using parameters that are
# *equivalent* to those used in `scipy.signal.welch` above
S_mlab, f_mlab = mlab.psd(
y, Fs=Fs, NFFT=nperseg, noverlap=(nperseg // 2),
detrend=None, scale_by_freq=True, window=mlab.window_hanning)
fig, axes = plt.subplots(2, 1, sharex=True)
# Plot PSD computed via both methods
axes[0].loglog(f_welch, S_welch, '-s')
axes[0].loglog(f_mlab, S_mlab, '-^')
axes[0].set_xlabel('f')
axes[0].set_ylabel('PSD')
axes[0].legend(['scipy.signal.welch', 'mlab.psd'], loc='upper left')
# Plot relative error
delta = np.abs(S_welch - S_mlab)
avg = 0.5 * np.abs(S_welch + S_mlab)
relative_error = delta/avg
axes[1].loglog(f_welch, relative_error, 'k')
axes[1].set_xlabel('f')
axes[1].set_ylabel('relative error')
plt.suptitle('nperseg = %i' % nperseg)
plt.show()
बहुत बढ़िया! धन्यवाद! मुझे "सममित" बनाम "आवधिक" विंडोिंग कार्यों के बारे में पता नहीं था, इसलिए पृष्ठभूमि की जानकारी का लिंक विशेष रूप से उपयोगी था :) –