2012-12-23 36 views
11

私はPythonのfftの機能を理解しようとしています。私が遭遇した奇妙なことの1つは、Parseval's theoremは約50私はこの機能については見つけることができるすべての情報があるとして、今、それは0PythonでのParsevalの定理

import numpy as np 
import matplotlib.pyplot as plt 
import scipy.fftpack as fftpack 

pi = np.pi 

tdata = np.arange(5999.)/300 
dt = tdata[1]-tdata[0] 

datay = np.sin(pi*tdata)+2*np.sin(pi*2*tdata) 
N = len(datay) 

fouriery = abs(fftpack.rfft(datay))/N 

freqs = fftpack.rfftfreq(len(datay), d=(tdata[1]-tdata[0])) 

df = freqs[1] - freqs[0] 

parceval = sum(datay**2)*dt - sum(fouriery**2)*df 
print parceval 

plt.plot(freqs, fouriery, 'b-') 
plt.xlim(0,3) 
plt.show() 

する必要がありますが、私はそれが正規化因子だとかなり確信して、私はそれを見つけることができるようには見えません。 scipy.fftpack.rfft documentation

答えて

15

あなたの正規化係数は、離散シーケンスへの連続信号のフーリエ変換のためのParsevalの定理を適用しようとしているから来ています。 the wikipedia article on the Discrete Fourier transformのサイドパネルには、フーリエ変換、フーリエ級数、離散フーリエ変換、およびDirac combsによるサンプリングの関係についての議論がいくつかあります。

は、 Parseval's theorem, when applied to DFTs、長い話を短くするために、統合を必要としませんが、合計: 2*piあなたは dtdfあなたの総和でmultipliyingによって作成されています。

scipy.fftpack.rfftを使用しているため、データのDFTは適切ではありませんが、負の値は対称であるため、データの正の半分だけが正しいことに注意してください。したがって、データの半分だけを追加しているため、DC用語の0には、2があり、4*piには、@unutbuが見つかりました。いずれの場合においても

datayはあなたの順序を保持している場合、次のように、あなたはパーセバルの定理を確認することができますが:

fouriery = fftpack.rfft(datay) 
N = len(datay) 
parseval_1 = np.sum(datay**2) 
parseval_2 = (fouriery[0]**2 + 2 * np.sum(fouriery[1:]**2))/N 
print parseval_1 - parseval_2 

二総和は、このような奇妙な形を取る必要はありませんscipy.fftpack.fftまたはnumpy.fft.fftを使用する:

fouriery_1 = fftpack.fft(datay) 
fouriery_2 = np.fft.fft(datay) 
N = len(datay) 
parseval_1 = np.sum(datay**2) 
parseval_2_1 = np.sum(np.abs(fouriery_1)**2)/N 
parseval_2_2 = np.sum(np.abs(fouriery_2)**2)/N 
print parseval_1 - parseval_2_1 
print parseval_1 - parseval_2_2 
+0

訂正ありがとうございます! – unutbu

+0

注記:これは、浮動小数点数が実数と同じではないという一般的な問題の一部です。 – Marcin

+0

@Marcin Yep、チェックで '=='を '-'に変更しました... – Jaime

関連する問題