4

numpy と matplotlib を使用して、シミュレーションからのデータ出力を分析しています。根源を見つけることができない (明らかな) 矛盾が 1 つあります。それは次のとおりです。

特定のエネルギー a^2~1 を持つ信号があります。rfft を使用して FFT を取得し、フーリエ空間でエネルギーを計算すると、かなり大きくなります。私のデータなどの詳細を無効にするために、単純な正弦波の例を次に示します。

from pylab import *
xx=np.linspace(0.,2*pi,128)
a=np.zeros(128)
for i in range(0,128):
    a[i]=sin(xx[i])
aft=rfft(a)
print mean(abs(aft)**2),mean(a**2)

原則として、両方の数値は（少なくとも数値的な意味では）同じでなければなりませんが、これがこのコードから得られるものです。

62.523081632 0.49609375

numpy.fft のドキュメントを調べてみましたが、何も見つかりませんでした。ここで検索すると次のことがわかりましたが、そこの説明を理解できませんでした。

既存の (合成された) 信号とフィルター処理された信号の間の大きな FFT 振幅差

私は何を見逃していますか/誤解していますか? この点に関するヘルプ/ポインタは大歓迎です。

ありがとう！

4

2 に答える 2

9

rfftヘンリーは非正規化の部分で正しいですが、あなたはではなくを使用しているため、もう少しありますfft。以下は彼の答えと一致しています。

>>> x = np.linspace(0, 2 * np.pi, 128)
>>> y = 1 - np.sin(x)
>>> fft = np.fft.fft(y)
>>> np.mean((fft * fft.conj()).real)
191.49999999999991
>>> np.mean(y**2)
1.4960937500000004
>>> fft = fft / np.sqrt(len(fft))
>>> np.mean((fft * fft.conj()).real)
1.4960937499999991

しかし、 で同じことを試してみると、うまくいきrfftません:

>>> rfft = np.fft.rfft(y)
>>> np.mean((rfft * rfft.conj()).real)
314.58462009358772
>>> rfft /= np.sqrt(len(rfft))
>>> np.mean((rfft * rfft.conj()).real)
4.8397633860551954
65
>>> np.mean((rfft * rfft.conj()).real) / len(rfft)
4.8397633860551954

ただし、以下は適切に機能します。

>>> (rfft[0] * rfft[0].conj() +
...  2 * np.sum(rfft[1:] * rfft[1:].conj())).real / len(y)
1.4960937873636722

あなたが得ているものを使用するrfftと、データのDFTは適切ではありませんが、負は対称であるため、データの正の半分のみです。平均を計算するには、DC コンポーネント以外のすべての値を 2 回考慮する必要があります。これは、コードの最後の行で行われます。

于 2013-03-01T00:17:19.383 に答える
3

ほとんどの FFT ライブラリでは、さまざまな DFT フレーバーが直交していません。numpy.fftライブラリは、逆変換中にのみ必要な正規化を適用します。

DFT のウィキペディアの説明を検討してください。逆 DFT には、DFT にはない 1/N 項があります (N は変換の長さです)。DFT の直交バージョンを作成するには、正規化されていない DFT の結果を 1/sqrt(N) でスケーリングする必要があります。この場合、変換は直交します (つまり、直交 DFT を F と定義すると、逆 DFT は F の共役またはエルミート転置になります)。

あなたの場合、単純にスケーリングすることで正しい答えを得ることができますaft( 1.0/sqrt(len(a))Nは変換の長さから見つかることに注意してください。実際のFFTは値の約半分を捨てるだけなので、a重要なのはその長さです)。

正規化を最後まで残す理由は、ほとんどの状況では問題にならないため、正規化を 2 回行う計算コストを節約できるためだと思います。実際、非常に高速なFFTW ライブラリは、どちらの方向にも正規化を行わず、完全にユーザーに任せています。

編集：明確にするために、上記の説明は完全に正しくありません。あなたの場合、DC成分は2回追加されますが1.0/sqrt(len(a))、ユニタリ変換を生成するための正しいスケーリングであるため、その単純なスケーリングでは正しい答えは得られません。

于 2013-02-28T23:42:34.317 に答える