40

301 フレームのムービー クリップから収集された 301 値の配列があります。これは、1 フレームから 1 つの値を意味します。ムービー クリップは 30 fps で実行されているため、実際の長さは 10 秒です

ここで、この「信号」のパワー スペクトルを取得したいと思います (右の軸を使用)。私は試した：

 X = fft(S_[:,2]);
 pl.plot(abs(X))
 pl.show()

私も試しました：

 X = fft(S_[:,2]);
 pl.plot(abs(X)**2)
 pl.show()

これが本当のスペクトルだとは思いませんが。

シグナル： ここに画像の説明を入力

スペクトル:ここに画像の説明を入力

パワースペクトル :

ここに画像の説明を入力

誰でもこれについて助けてもらえますか? Hz でプロットしたいと思います。

4

5 に答える 5

60

Numpy には、np.fft.fftfreqFFT コンポーネントに関連付けられた周波数を計算するための便利な関数があります。

from __future__ import division
import numpy as np
import matplotlib.pyplot as plt

data = np.random.rand(301) - 0.5
ps = np.abs(np.fft.fft(data))**2

time_step = 1 / 30
freqs = np.fft.fftfreq(data.size, time_step)
idx = np.argsort(freqs)

plt.plot(freqs[idx], ps[idx])

ここに画像の説明を入力

あなたのケースで見られる最大の周波数は 30 Hz ではありませんが、

In [7]: max(freqs)
Out[7]: 14.950166112956811

パワー スペクトルにサンプリング周波数が表示されることはありません。サンプル数が偶数の場合、ナイキスト周波数に達したことになりますが、この場合は 15 Hz です (ただし、numpy では -15 と計算されます)。

于 2013-03-13T14:39:01.677 に答える
21

rate がサンプリング レート (Hz) の場合、np.linspace(0, rate/2, n)は fft の各ポイントの周波数配列です。rfftデータの fft を計算するために使用できるのは実数値です。

import numpy as np
import pylab as pl
rate = 30.0
t = np.arange(0, 10, 1/rate)
x = np.sin(2*np.pi*4*t) + np.sin(2*np.pi*7*t) + np.random.randn(len(t))*0.2
p = 20*np.log10(np.abs(np.fft.rfft(x)))
f = np.linspace(0, rate/2, len(p))
plot(f, p)

ここに画像の説明を入力

信号 x には 4Hz と 7Hz の正弦波が含まれているため、4Hz と 7Hz に 2 つのピークがあります。

于 2013-03-13T12:37:53.100 に答える
10

scipy.signal.welchを使用して、ウェルチ法を使用してパワー スペクトル密度を推定することもできます。np.fft.fft と scipy.signal.welch の比較は次のとおりです。

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

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 += np.random.normal(scale=np.sqrt(noise_power), size=time.shape)

# np.fft.fft
freqs = np.fft.fftfreq(time.size, 1/fs)
idx = np.argsort(freqs)
ps = np.abs(np.fft.fft(x))**2
plt.figure()
plt.plot(freqs[idx], ps[idx])
plt.title('Power spectrum (np.fft.fft)')

# signal.welch
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.title('Power spectrum (scipy.signal.welch)')
plt.show()

[fft[2] ウェルチ

于 2018-07-31T13:26:31.143 に答える
6

numpy fftページからhttp://docs.scipy.org/doc/numpy/reference/routines.fft.html：

入力aが時間領域信号でA=fft（a）の場合、np.abs（A）はその振幅スペクトルであり、np.abs（A）**2はそのパワースペクトルです。位相スペクトルはnp.angle（A）によって取得されます。

于 2013-03-13T10:06:00.320 に答える