Scipy寻峰函数和距离参数

velaa5lx  于 2023-06-23  发布在  其他
关注(0)|答案(1)|浏览(101)

我正在与scipy的find peaks功能作斗争。我有一个大约20 k点的数据集。X轴是从0秒到20秒的以秒计的时间,y轴表示从-2mb到5 mb的压力测量值。我的信号是一个正弦,几乎是一个完美的,除了一些振动(不是噪音,只是在正弦最大值之前和/或之后的压力的快速波动)。这意味着我对两个极大值之间的距离有了很好的估计。为了提取峰值,我使用scipy.signal.find_peaks(y, height=0, distance=500, prominence=(0.1, 1)),高度等于0,因为我想要正峰值,突出是有用的去除振动峰值和距离,据我所知,它是有用的去除两个正弦峰值之间的局部最大值。不幸的是,我不明白为什么,但似乎距离不是两个最大值之间的点数。实际上,对于第一个命令,结果是array([ 3928, 5047, 6243, 11451, 11955, 13737, 19300]。如果我将距离设置为等于1000,第一个极大值应该仍然出现,因为它与第二个极大值之间的距离超过1000点,对吗?但在这种情况下,我得到以下结果array([11955, 19300]
我错过了什么?
谢谢你的帮助:)

EDIT. FFT代码此代码返回等于1.58的幅度,这显然是不正确的,因为最大值约为4.2,最小值约为-2.2。

import numpy as np
import matplotlib.pyplot as plt

# Generate the time values
duration = 1.0  # Duration of the signal in seconds
sampling_rate = 1000  # Number of samples per second
t = x
signal = y

# Apply the FFT to the signal
fft_result = np.fft.fft(signal)

#number of samples per second
sampling_rate = 1000

# Calculate the amplitudes of the frequencies
frequencies = np.fft.fftfreq(len(t), 1.0 / sampling_rate)
amplitudes = np.abs(fft_result) / len(t)

# Find the index of the main frequency component
main_frequency_index = np.argmax(amplitudes)

# Get the amplitude of the main frequency component
main_frequency_amplitude = amplitudes[main_frequency_index]

# Print the result
print("Amplitude of the main frequency component:", main_frequency_amplitude)

# Plot the signal and its frequency spectrum
fig, (ax1, ax2) = plt.subplots(2, 1)
ax1.plot(t, signal)
ax1.set_xlabel('Time')
ax1.set_ylabel('Amplitude')
ax1.set_title('Signal')
ax2.stem(frequencies, amplitudes)
ax2.set_xlabel('Frequency')
ax2.set_ylabel('Amplitude')
ax2.set_title('Frequency Spectrum')
plt.tight_layout()
plt.show()
js81xvg6

js81xvg61#

一个非常基本的FFT分析如下所示

import matplotlib.pyplot as plt
import numpy as np

irregular_time, irregular_amplitude = np.loadtxt(
    fname='Data - Feuille 1.csv',
    delimiter=',',
).T

time = np.linspace(
    start=irregular_time.min(),
    stop=irregular_time.max(),
    num=irregular_time.size,
)
timestep = time[1] - time[0]
amplitude = np.interp(
    x=time,
    xp=irregular_time,
    fp=irregular_amplitude,
)

spectrum = np.fft.rfft(amplitude, norm='forward')
spectral_amplitude = np.abs(spectrum)
freqs = np.fft.rfftfreq(n=time.size, d=timestep)

peak_index = spectral_amplitude.argmax()
peak_freq = freqs[peak_index]
peak_amplitude = 2*spectral_amplitude[peak_index]
peak_phase = np.angle(spectrum[peak_index])

print('Amplitude:', peak_amplitude)
print('Frequency (approx, Hz):', peak_freq)
print('Phase (deg):', np.rad2deg(peak_phase))

reconstructed = peak_amplitude * np.cos(
    2*np.pi*peak_freq*time + peak_phase
) + np.abs(spectrum[0])

ax_top: plt.Axes
ax_bottom: plt.Axes
fig, (ax_top, ax_bottom) = plt.subplots(nrows=2)

ax_top.plot(time, amplitude, label='experiment')
ax_top.plot(time, reconstructed, label='reconstructed')
ax_top.set_xlabel('Interpolated time (s)')
ax_top.set_ylabel('Amplitude')
ax_top.legend(loc='lower right')

ax_bottom.semilogx(freqs, spectral_amplitude)
ax_bottom.set_xlabel('Frequency (Hz)')
ax_bottom.set_ylabel('Amplitude')

plt.show()
Amplitude: 3.183337345334171
Frequency (approx, Hz): 0.1499400170118757
Phase (deg): 86.4919931694836

该频率仅为近似值,因为在样本中循环次数较少的情况下,该频率相当低。这可以通过例如改进。谱峰估计,但不清楚这是需要的;也可以用更长的样本来改进。

相关问题