numpy和scipy模块计算出的数据的PSD结果差异

问题描述 投票:0回答:0
  1. 对于一维数据的计算,我对比了PSD numpy和scipy的结果,发现计算结果 np.fft.fft() 和 scipy.signal.periodogram() 是相同的,但是 scipy.signal.welch() 的结果并不完全不同。
  2. 对于一维时间序列数据,不像 scipy.signal() 可以直接获取PSD值的模块,np.ft.ft()需要 计算输入序列的快速傅立叶变换 (FFT) 和 返回复值谱,并计算谱 将值平方并除以归一化参数以获得 PSD 值。但是,在这个规范化过程中,我不确定 用于规范化的参数。归一化 我看到的参数包括样本总数n,平方 样本总数n ** 2,和总数 samples乘以sampling Frequency n*fs,没找到一个太 权威定义
  3. 另外,二维矩阵数据的PSD计算 也涉及相同的相同归一化参数 尺寸数据。

这里是相关的实验代码: 科学

import numpy as np
from scipy.signal import periodogram, welch
import matplotlib.pyplot as plt

# Generate random data
fs = 4000
t = np.arange(0, 3, 1/fs)
y = np.sin(2*np.pi*50*t) + np.cos(2*np.pi*140*t) + np.cos(2*np.pi*300*t) * t ** 2

# Calculate PSD using signal.welch
# f, psd = welch(x, fs=fs, nperseg=256, scaling='density')
f, psd = periodogram(y, fs, scaling='density')

# Plot the PSD versus frequencies
plt.plot(f, psd)
plt.xlabel('Frequency (Hz)')
plt.ylabel('PSD')
plt.show()

numpy

# Calculate PSD using np.fft
n = len(y)
y_FFT = np.fft.fft(y)
# psd = np.abs(y_FFT)**2 / (n*fs)
freq = np.fft.fftfreq(n, d=1/fs)
# freq = np.fft.fftshift(freq)
# print(freq)
# Find the indices corresponding to positive frequencies
positive_freq_indices = np.where(freq >= 0)
# Extract the PSD for positive frequencies
psd = np.abs(y_FFT[positive_freq_indices])**2 / (n * fs)
# Plot PSD
plt.plot(freq[positive_freq_indices], psd)
plt.xlabel('Frequency')
plt.ylabel('PSD')
plt.show()

根据不同的PSD结果进行频谱分析,会得出完全不同的结论。 结果绘图可视化 y值随时间变化如图 y value changes

scipy.周期图 scipy.periodogram

scipy.welch scipy.welch

numpy.fft numpy.fft

不确定PSD计算和归一化参数选择,二维空间坐标轴的计算是否正确

from scipy.fft import fft2, fftfreq

fs = 8000
fft_coefficients = fft2(Two_dim_data)  # Compute the 2D FFT

# Compute the power spectrum by squaring the magnitude of the Fourier coefficients
power_spectrum = np.abs(fft_coefficients) ** 2

# Scale the power spectrum by the appropriate normalization factor
# Normalize by the number of elements in the matrix
two_dim_psd = power_spectrum / (Two_dim_data.shape[0] * fs * Two_dim_data.shape[1] * fs)

# Calculate the spatial frequencies
freq_x = fftfreq(Two_dim_data.shape[0], 1/fs)
freq_y = fftfreq(Two_dim_data.shape[1], 1/fs)

reference_1 reference_2

我想知道为什么得到的PSD结果差异如此之大,以及在从功率计算功率密度时使用什么归一化参数。

python numpy scipy fft psd
© www.soinside.com 2019 - 2024. All rights reserved.