3

Python で 2D FFT を実行するときの正規化に関する簡単な質問があります。私の理解では、正規化係数は、配列を 1 で埋めることから決定できます。

たとえば、1d では、[1,1,1,1] の FFT で [4+0j,0+0j,0+0j,0+0j] が得られるため、正規化係数は 1/N=1/4 にする必要があります。

2D では、[[1,1],[1,1]] の FFT は [[4+0j,0+0j],[0+0j,0+0j]] を与えるので、正規化は 1/MN になるはずです=1/(2*2)=1/4。

ここで、3000 x 3000 の行列があり、各要素が平均 0 のガウス分布値を持つとします。これを FFT して正規化すると (正規化係数 = 1/(3000*3000))、次数 10^- の平均べき乗が得られます。 7。

ここで、1000 x 1000 要素のサブ領域 (正規化係数 = 1/(1000*1000)) を使用してこれを繰り返します。これから得られる平均電力は、10^-6 のオーダーです。なぜ10倍の差があるのか​​ 疑問に思っています。平均電力は同じであるべきではありませんか？追加の正規化係数が不足していますか?

係数の差が事実上 9 であると言う場合、これは要素の数 (3000 x 3000 には 1000 x 1000 の 9 倍の要素があります) から来ていると推測できますが、この余分な要素の直感的な理由は何でしょうか? また、「真の」基礎となる平均検出力を得るために、絶対正規化係数をどのように決定すればよいでしょうか?

どんな洞察も大歓迎です。前もって感謝します！

サンプルコード:

import numpy as np
a   = np.random.randn(3000,3000)
af  = np.fft.fft2(a)/3000.0/3000.0
aP  = np.mean(np.abs(af)**2)

b   = a[1000:2000,1000:2000]
bf  = np.fft.fft2(b)/1000.0/1000.0
bP  = np.mean(np.abs(bf)**2)

print aP,bP

>1.11094908545e-07 1.00226264535e-06
4

1 に答える 1

3

まず、この問題は 1D FFT と 2D FFT の違いとは関係がなく、総電力と平均電力が配列内の要素数にどのように対応するかに関係していることに注意することが重要です。

a9 の因数は の 9倍の要素に由来すると言うとき、あなたは正確に正しいですb。np.fft.fft2(a)/3000./3000.紛らわしいのは、 と で除算して既に正規化していることに気付いたことでしょう。np.fft.fft2(b)/1000./1000.実際、これらの正規化は、空間ドメインと周波数ドメインで合計 (平均ではない) パワーを等しくするために必要です。平均を得るには、配列のサイズで再度割る必要があります。

あなたの質問は、2 つのドメイン (空間/時間と周波数) の総電力が等しいというパーセヴァルの定理に関するものです。そのステートメント、DFT はthisです。右側の 1/N にもかかわらず、これは平均電力ではなく、総電力であることに注意してください。1/N の理由は、DFT の正規化規則です。

Python で言えば、これは、時間/空間信号sigの場合、Parseval の等価性は次のように記述できることを意味します。

np.sum(np.abs(sig)**2) == np.sum(np.abs(np.fft.fft(sig))**2)/sig.size

以下は完全な例です。おもちゃのケース (1 が 1 で埋められた 1 次元および 2 次元の配列) から始まり、独自のケースで終わります。.size配列内の要素の総数を返す numpy.ndarrayのプロパティを使用したことに注意してください。あなたの/1000./1000.etc と同等です。これが役立つことを願っています!

import numpy as np

print 'simple examples:'

# 1-d, 4 elements:
ones_1d = np.array([1.,1.,1.,1.])
ones_1d_f = np.fft.fft(ones_1d)

# compute total power in space and frequency domains:
space_power_1d = np.sum(np.abs(ones_1d)**2)
freq_power_1d = np.sum(np.abs(ones_1d_f)**2)/ones_1d.size
print 'space_power_1d = %f'%space_power_1d
print 'freq_power_1d = %f'%freq_power_1d

# 2-d, 4 elements:
ones_2d = np.array([[1.,1.],[1.,1.]])
ones_2d_f = np.fft.fft2(ones_2d)

# compute and print total power in space and frequency domains:
space_power_2d = np.sum(np.abs(ones_2d)**2)
freq_power_2d = np.sum(np.abs(ones_2d_f)**2)/ones_2d.size
print 'space_power_2d = %f'%space_power_2d
print 'freq_power_2d = %f'%freq_power_2d

# 2-d, 9 elements:
ones_2d_big = np.array([[1.,1.,1.],[1.,1.,1.],[1.,1.,1.]])
ones_2d_big_f = np.fft.fft2(ones_2d_big)

# compute and print total power in space and frequency domains:
space_power_2d_big = np.sum(np.abs(ones_2d_big)**2)
freq_power_2d_big = np.sum(np.abs(ones_2d_big_f)**2)/ones_2d_big.size
print 'space_power_2d_big = %f'%space_power_2d_big
print 'freq_power_2d_big = %f'%freq_power_2d_big
print


# asker's example array a and fft af:
print 'askers examples:'
a = np.random.randn(3000,3000)
af = np.fft.fft2(a)

# compute the space and frequency total powers:
space_power_a = np.sum(np.abs(a)**2)
freq_power_a = np.sum(np.abs(af)**2)/af.size

# mean power is the total power divided by the array size:
freq_power_a_mean = freq_power_a/af.size

print 'space_power_a = %f'%space_power_a
print 'freq_power_a = %f'%freq_power_a
print 'freq_power_a_mean = %f'%freq_power_a_mean
print
# the central 1000x1000 section of the 3000x3000 original array:
b = a[1000:2000,1000:2000]
bf = np.fft.fft2(b)

# we expect the total power in the space and frequency domains 
# to be about 1/9 of the total power in the space frequency domains 
# for matrix a:
space_power_b = np.sum(np.abs(b)**2)
freq_power_b = np.sum(np.abs(bf)**2)/bf.size
# we expect the mean power to be the same as the mean power from
# matrix a:
freq_power_b_mean = freq_power_b/bf.size

print 'space_power_b = %f'%space_power_b
print 'freq_power_b = %f'%freq_power_b
print 'freq_power_b_mean = %f'%freq_power_b_mean
于 2014-09-08T18:49:21.930 に答える