質問: 'gray') plt.show() # 表示 plt.plot(phase, ampr) plt.grid() plt.show() plt.plot(amp) plt.grid() ``` パワースペクトルの角度方向はどのように抽出すればよいでしょうか? アドバイスよろしくお願いします.
回答: ご質問の図からの想像ですが、多分以下のレポートを参照しておられるのではないでしょうか? この中に > > q(θ)は水平軸から角度θの微小な線形領域内のパワースペクトルの和を表す。 > > > とありますので、DFTの結果の中からq(θ)の範囲を示す扇型の内側の点のみを選択的に集計すればよいと思います。 以下、煩雑さを避けるために解析対象の画像が正方形であり、垂直方向の符号は数学でおなじみの座標系(複素平面)に合わせる前提とします。 画像の辺の画素数をWとすると、dftshift後のスペクトル配列spectrumの要素spectrumW // 2 + re, W // 2 - im]は、複素数 c = re + 1j * im を用いて以下のような波の成分(スペクトル)となります。 ・方向: cの偏角 ・波長: W/|c| (単位は画素) numpyには複素数の偏角を求めるangle、複素数の絶対値(複素平面原点からの距離)を求めるabsなどの関数があります。これらの機能を利用すると次のように集計できます。(本件では同じ波長範囲を集計すればよく、コードでは波長ではなく周波数で範囲をチェックしています) ``` import numpy as np # 複素平面を構成する値の配列を生成(原点は中心) def complex_plane(width): half_width = width // 2 re = np.array(range(width)) - half_width im = - re re, im = np.meshgrid(re, im) return re + im * 1j # スペクトルを与えられた角度の範囲内だけ集計 def aggregate_in_angle(agg_fun, spectrum, min_angle, max_angle): width = spectrum.shape[0] assert spectrum.shape[1] == width min_radius, max_radius = 0, width // 2 assert 0 <= min_angle < max_angle <= 180 # 単位はdegree cp = complex_plane(width) cp_mag = np.abs(cp) in_radius = np.logical_and(min_radius < cp_mag, cp_mag < max_radius) cp_angle = np.angle(cp, deg=True) in_angle = np.logical_and(min_angle <= cp_angle, cp_angle < max_angle) return agg_fun(spectrum[np.logical_and(in_radius, in_angle)]) def test(spectrum, angle_gap): # 指定角度ごとの平均を求める means = [aggregate_in_angle(np.mean, spectrum, angle, angle + angle_gap) for angle in range(0, 180, angle_gap)] return means ``` ちなみにレポートでは「パワースペクトルの和」と書いてありますが・・・ * 集計方法 DFTの結果を集計するなら和ではなく平均の方がよい気がします。連続関数の積分ならまだしも離散的スペクトルを集計するのですから扇型の範囲に含まれるデータ数が角度によって増減することを考慮した方がよいのではないでしょうか。 * パワースペクトル 質問者さんのコードにあるデシベル(?)換算は「ダイナミックレンジが大きなデータの傾向を目視で確認しやすくするための処理」であって「集計などのデータ処理をする前に施すべき変換ではない」気がします。 自然画像では結果を評価しにくいので単純なスペクトルを持つサンプル画像を生成して上記を試してみました。 ![イメージ説明 ちなみに、いろいろな波の成分で試してみたところ、fft結果のスペクトルの絶対値に対して平均や最大値をとるとそこそこの結果が出るように見えました。ただ、目的によっては集計方法を変えたほうがよいかも知れません。そのあたりの知識が不十分なので自分にははっきりしたことが言えませんが。