welch#
- scipy.signal.welch(x, fs=1.0, window='hann', nperseg=None, noverlap=None, nfft=None, detrend='constant', return_onesided=True, масштабирование='density', ось=-1, среднее='mean')[источник]#
Оценка спектральной плотности мощности с использованием метода Уэлча.
Метод Уэлча [1] вычисляет оценку спектральной плотности мощности путем разделения данных на перекрывающиеся сегменты, вычисления модифицированного периодограмма для каждого сегмента и усреднения периодограмм.
- Параметры:
- xarray_like
Временной ряд значений измерений
- fsfloat, опционально
Частота дискретизации x временной ряд. По умолчанию 1.0.
- windowstr или tuple или array_like, опционально
Желаемое окно для использования. Если window является строкой или кортежем, он передаётся в
get_windowдля генерации значений окна, которые по умолчанию являются DFT-чётными. См.get_windowдля списка окон и требуемых параметров. Если window если это array_like, он будет использоваться напрямую как окно, и его длина должна быть nperseg. По умолчанию используется окно Ханна.- npersegint, необязательный
Длина каждого сегмента. По умолчанию None, но если window является строкой или кортежем, устанавливается в 256, а если window является array_like, устанавливается в длину окна.
- noverlapint, необязательный
Количество точек перекрытия между сегментами. Если None,
noverlap = nperseg // 2. По умолчанию None.- nfftint, необязательный
Длина БПФ, используемая, если требуется дополненное нулями БПФ. Если None, длина БПФ равна nperseg. По умолчанию None.
- detrendstr или функция или False, опционально
Определяет, как удалять тренд в каждом сегменте. Если
detrendявляется строкой, она передается как тип аргумент дляdetrendфункция. Если это функция, она принимает сегмент и возвращает сегмент с удалённым трендом. Еслиdetrendявляется False, детрендинг не выполняется. По умолчанию 'constant'.- return_onesidedbool, необязательно
Если True, возвращает односторонний спектр для реальных данных. Если False возвращает двусторонний спектр. По умолчанию True, но для комплексных данных всегда возвращается двусторонний спектр.
- масштабирование{ 'density', 'spectrum' }, опционально
Выбор между вычислением спектральной плотности мощности ('density') где Pxx имеет единицы измерения V**2/Гц и вычисление квадрата амплитудного спектра ('spectrum'), где Pxx имеет единицы измерения V**2, если x измеряется в В, а fs измеряется в Гц. По умолчанию ‘density’
- осьint, необязательный
Ось, вдоль которой вычисляется периодограмма; по умолчанию по последней оси (т.е.
axis=-1).- среднее{ ‘mean’, ‘median’ }, опционально
Метод усреднения периодограмм. По умолчанию 'mean'.
Добавлено в версии 1.2.0.
- Возвращает:
- fndarray
Массив частот выборки.
- Pxxndarray
Спектральная плотность мощности или спектр мощности x.
Смотрите также
csdПерекрёстная спектральная плотность мощности методом Уэлча
periodogramПростой, опционально модифицированный периодограмм
lombscargleПериодограмма Ломба-Скаргла для неравномерно дискретизированных данных
Примечания
Подходящее количество перекрытия будет зависеть от выбора окна и ваших требований. Для окна Ханна по умолчанию перекрытие 50% является разумным компромиссом между точной оценкой мощности сигнала и избеганием повторного учета данных. Более узкие окна могут требовать большего перекрытия. Если noverlap равен 0, этот метод эквивалентен методу Бартлетта [2].
Отношение квадрата величины (
scaling='spectrum') деленная на спектральную плотность мощности (scaling='density') является постоянным множителемsum(abs(window)**2)*fs / abs(sum(window))**2. Если return_onesided являетсяTrue, значения отрицательных частот добавляются к значениям соответствующих положительных.Обратитесь к Спектральный анализ раздел Руководство пользователя SciPy для обсуждения масштабирования спектральной плотности мощности и (квадрата) амплитудного спектра.
Добавлено в версии 0.12.0.
Ссылки
[1]P. Welch, "Использование быстрого преобразования Фурье для оценки спектров мощности: метод, основанный на усреднении по времени коротких модифицированных периодограмм", IEEE Trans. Audio Electroacoust. том 15, стр. 70-73, 1967.
[2]M.S. Bartlett, “Periodogram Analysis and Continuous Spectra”, Biometrika, vol. 37, pp. 1-16, 1950.
Примеры
>>> import numpy as np >>> from scipy import signal >>> import matplotlib.pyplot as plt >>> rng = np.random.default_rng()
Сгенерировать тестовый сигнал, синусоиду 2 Vrms на частоте 1234 Гц, искажённую белым шумом 0.001 V**2/Гц, дискретизированную на 10 кГц.
>>> 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
Если теперь ввести разрыв в сигнал, увеличив амплитуду небольшой части сигнала на 50, можно увидеть искажение средней спектральной плотности мощности, но использование медианного среднего лучше оценивает нормальное поведение.
>>> x[int(N//2):int(N//2)+10] *= 50. >>> f, Pxx_den = signal.welch(x, fs, nperseg=1024) >>> f_med, Pxx_den_med = signal.welch(x, fs, nperseg=1024, average='median') >>> plt.semilogy(f, Pxx_den, label='mean') >>> plt.semilogy(f_med, Pxx_den_med, label='median') >>> plt.ylim([0.5e-3, 1]) >>> plt.xlabel('frequency [Hz]') >>> plt.ylabel('PSD [V**2/Hz]') >>> plt.legend() >>> plt.show()
