scipy进行谱分析

关注

原创最新推荐文章于 2026-10-08 09:31:24 发布 1.1k 阅读 AI 写同款· GEO 优化 ›

Pytorch 实现 PSD 功率谱计算：periodogram 和 welch 方法前言由于跑代码需要，本人用 Pytorch 简单实现了 PSD 估计中的最常用的 periodogram 和 welch 函数。这里把我的实现分享给大家。在实现方案和参数的设计上，本人参考了 scipy 的 API 设计。 periodogram 的实现函数体如下。 def _periodogram(X: torch.Tensor, fs, detrend, scaling): if X.dim() > 2: X = torch.squeeze(X) elif X 阅读详情

periodogram

scipy.signal.periodogram(x, fs=1.0, window='boxcar', nfft=None, detrend='constant', return_onesided=True, scaling='density', axis=- 1)[source]
使用周期图估计功率谱密度。

参数
x:array_like
测量值的时间序列

fs:float,可选
x时间序列的采样频率。默认为1.0。

Window:str或tuple或array_like，可选
想要使用的窗口。如果window是字符串或元组，则将它传递给get_window以生成窗口值，即使在默认情况下也是dft。有关窗口和所需参数的列表，请参阅get_window。如果window是array_like，它将直接用作窗口，其长度必须为nperseg。默认为' boxcar '。

nfft:int,可选
FFT的长度。如果为None，则使用x的长度。

detrend:str或function或False，可选
指定如何为每个段趋势化。如果detrend是一个字符串，它将作为类型参数传递给detrend函数。如果它是一个函数，它接受一个段并返回一个有趋势的段。如果detrend为False，则不进行去趋势处理。默认为' constant '。

return_onesided:bool,可选
如果为True，返回真实数据的单边频谱。如果False返回一个双面频谱。默认为True，但是对于复杂的数据，总是返回一个双面谱。

scaling{ ‘density’, ‘spectrum’ }, optional
在计算功率谱密度(‘密度’)和计算功率谱(‘频谱’)之间进行选择，其中Pxx的单位是V**2/Hz，其中Pxx的单位是V**2，如果x的单位是V, fs的单位是Hz。默认为' density '

axis:int,可选
计算周期图的轴;默认值是在最后一个轴上(即轴=-1)。

返回
f :ndarray
采样频率数组。

Pxx :ndarray
x的功率谱密度或功率谱。

示例：

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

rng = np.random.default_rng()
# Generate a test signal, a 2 Vrms sine wave at 1234 Hz, corrupted by 0.001 V**2/Hz of white noise sampled at 10 kHz.

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)
# Compute and plot the power spectral density.

f, Pxx_den = signal.periodogram(x, fs)
plt.semilogy(f, Pxx_den)
plt.ylim([1e-7, 1e2])
plt.xlabel('frequency [Hz]')
plt.ylabel('PSD [V**2/Hz]')
plt.show()

welch

scipy.signal.welch(x, fs=1.0, window='hann', nperseg=None, noverlap=None, nfft=None, detrend='constant', return_onesided=True, scaling='density', axis=- 1, average='mean')
用韦尔奇方法估计功率谱密度。
Welch的方法[1]通过将数据划分为重叠的段，计算每个段的修正周期图并对周期图进行平均来计算功率谱密度的估计值。

参数
x:array_like
测量值的时间序列

fs:float,可选
x时间序列的采样频率。默认为1.0。

Window:str或tuple或array_like，可选
想要使用的窗口。如果window是字符串或元组，则将它传递给get_window以生成窗口值，即使在默认情况下也是dft。有关窗口和所需参数的列表，请参阅get_window。如果window是array_like，它将直接用作窗口，其长度必须为nperseg。默认为Hann窗口。

nperseg:int,可选
每段的长度。默认为None，但如果window是str或tuple，则设置为256，如果window是array_like，则设置为窗口的长度。

noverlap:int,可选
分段之间重叠的点数。如果为None, noverlap = nperseg // 2。默认为None。

nfft:int,可选
如果需要补零FFT，则使用的FFT的长度。如果为None，则FFT长度为nperseg。默认为None。

detrend:str或function或False，可选
指定如何为每个段趋势化。如果detrend是一个字符串，它将作为类型参数传递给detrend函数。如果它是一个函数，它接受一个段并返回一个有趋势的段。如果detrend为False，则不进行去趋势处理。默认为' constant '。


return_onesided:bool,可选
如果为True，返回真实数据的单边频谱。如果False返回一个双面频谱。默认为True，但是对于复杂的数据，总是返回一个双面谱。

scaling{ ‘density’, ‘spectrum’ }, optional
在计算功率谱密度(‘密度’)和计算功率谱(‘频谱’)之间进行选择，其中Pxx的单位是V**2/Hz，其中Pxx的单位是V**2，如果x的单位是V, fs的单位是Hz。默认为' density '

axis:int,可选
计算周期图的轴;默认值是在最后一个轴上(即轴=-1)。

average{ ‘mean’, ‘median’ }, optional
周期图平均时使用的方法。默认为' mean '。

返回
f:ndarray
采样频率数组。

Pxx:ndarray
x的功率谱密度或功率谱。

note:
适当的重叠量取决于窗口的选择和您的需求。对于默认的Hann窗口，50%的重叠是一个合理的权衡，既能准确估计信号功率，又不会过多计算任何数据。较窄的窗口可能需要较大的重叠。
如果超叠为0，此方法等价于巴特利特方法[2]。

示例：

from scipy import signal
import matplotlib.pyplot as plt
rng = np.random.default_rng()
# Generate a test signal, a 2 Vrms sine wave at 1234 Hz, corrupted by 0.001 V**2/Hz # of white noise sampled at 10 kHz.

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)
Compute and plot the power spectral density.

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()

csd

scipy.signal.csd(x, y, fs=1.0, window='hann', nperseg=None, noverlap=None, nfft=None, detrend='constant', return_onesided=True, scaling='density', axis=- 1, average='mean')
使用Welch的方法估计交叉功率谱密度Pxy。


参数
x:array_like
测量值的时间序列

y:array_like
测量值的时间序列

fs:float,可选
x和y时间序列的采样频率。默认为1.0。

Window:str或tuple或array_like，可选
想要使用的窗口。如果window是字符串或元组，则将它传递给get_window以生成窗口值，即使在默认情况下也是dft。有关窗口和所需参数的列表，请参阅get_window。如果window是array_like，它将直接用作窗口，其长度必须为nperseg。默认为Hann窗口。

nperseg:int,可选
每段的长度。默认为None，但如果window是str或tuple，则设置为256，如果window是array_like，则设置为窗口的长度。

Noverlap: int，可选
分段之间重叠的点数。如果为None, noverlap = nperseg // 2。默认为None。

nfft:int,可选
如果需要补零FFT，则使用的FFT的长度。如果为None，则FFT长度为nperseg。默认为None。

detrend:str或function或False，可选
指定如何为每个段趋势化。如果detrend是一个字符串，它将作为类型参数传递给detrend函数。如果它是一个函数，它接受一个段并返回一个有趋势的段。如果detrend为False，则不进行去趋势处理。默认为' constant '。

return_onesided:bool,可选
如果为True，返回真实数据的单边频谱。如果False返回一个双面频谱。默认为True，但是对于复杂的数据，总是返回一个双面谱。

scaling{ ‘density’, ‘spectrum’ }, optional
选择计算交叉谱密度(' density ')，其中Pxy的单位是V**2/Hz，以及计算交叉谱(' spectrum ')，其中Pxy的单位是V**2，如果x和y的单位是V, fs的单位是Hz。默认为' density '

axis:int,可选
计算两个输入的CSD的轴;默认值是在最后一个轴上(即轴=-1)。

average{ ‘mean’, ‘median’ }, optional
周期图平均时使用的方法。如果谱是复的，则分别计算实部和虚部的平均值。默认为' mean '。

返回
f:ndarray
采样频率数组。

Pxy:ndarray
x,y的交叉谱密度或交叉功率谱。

note:
By convention, Pxy is computed with the conjugate FFT of X multiplied by the FFT of Y.

If the input series differ in length, the shorter series will be zero-padded to match.

按照惯例，Pxy是用X的共轭FFT乘以Y的共轭FFT来计算的。
如果输入序列长度不同，则较短的序列将补零以匹配。
适当的重叠量取决于窗口的选择和您的需求。对于默认的Hann窗口，50%的重叠是一个合理的权衡，既能准确估计信号功率，又不会过多计算任何数据。较窄的窗口可能需要较大的重叠。

示例：

from scipy import signal
import matplotlib.pyplot as plt
rng = np.random.default_rng()
#Generate two test signals with some common features.

fs = 10e3
N = 1e5
amp = 20
freq = 1234.0
noise_power = 0.001 * fs / 2
time = np.arange(N) / fs
b, a = signal.butter(2, 0.25, 'low')
x = rng.normal(scale=np.sqrt(noise_power), size=time.shape)
y = signal.lfilter(b, a, x)
x += amp*np.sin(2*np.pi*freq*time)
y += rng.normal(scale=0.1*np.sqrt(noise_power), size=time.shape)
# Compute and plot the magnitude of the cross spectral density.

f, Pxy = signal.csd(x, y, fs, nperseg=1024)
plt.semilogy(f, np.abs(Pxy))
plt.xlabel('frequency [Hz]')
plt.ylabel('CSD [V**2/Hz]')
plt.show()