やりたいこと
①csvファイルの数値を読み込んで,二次元フーリエ変換(2DFFT)を行い,計算結果をcsvファイルで出力したいです.
②①で得られた結果から動径方向分布 p(r) と角度方向分布 q(θ)のスペクトル分布を算出したいです.
困っている事
②の2DFFTから動径分布を計算を算出する部分.
CSVファイルの中身
CSVファイルの中身はこのようになっています.
前回の質問のcsvファイルとは異なります.https://teratail.com/questions/ds77z65aj8c652
余分な部分を削除し,データの縦×横を64×64に揃えました.
試したこと
①2DFFTのプログラム(python)の改造
pythonで書かれた画像の2DFFTのプログラム https://watlab-blog.com/2020/03/22/2d-fft/
画像を読み込む部分をcsvに置き換え, 角度方向分布,動軸方向分布のグラフを算出するプログラムを作りました.
Python
import pandas as pd import numpy as np import matplotlib.pyplot as plt def main(): # csvを読み込み、2Dフーリエ変換をする df = pd.read_csv('pivdata.csv') # csvを読み込み f = np.fft.fft2(df) # 2Dフーリエ変換 f_shift = np.fft.fftshift(f) # 直流成分を中心に移動させるためN/2シフトさせる mag = 20 * np.log(np.abs(f_shift)) # パワースペクトルの計算 #pd.DataFrame(f_shift).to_csv('out.csv') # フーリエ変換の結果をcsvに保存 # パワースペクトルの表示 plt.figure() plt.plot(mag) # 角度方向分布の表示 (Angular Distribution Function) adf = average_angle(mag, 10) # average_angleに渡す plt.figure() plt.plot(range(-180,180,10),adf) plt.show() # 動径方向分布の表示 (Radial Distribution Function) rdf = average_radius(mag, 1) # average_radiusに渡す print(rdf) plt.figure() plt.plot(range(0, 32, 1),rdf) plt.show # 複素平面を構成する値の配列を生成(原点は中心) 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] min_radius, max_radius = 0, width // 2 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 average_angle(spectrum, angle_gap): means = [aggregate_in_angle(np.mean, spectrum, angle, angle + angle_gap) for angle in range(-180, 180, angle_gap)] return means # スペクトルを与えられた半径の範囲内だけ集計 def aggregate_in_radius(agg_fun, spectrum, min_radius, max_radius): width = spectrum.shape[0] min_radius, max_radius = 0, width // 2 cp = complex_plane(width) cp_mag = np.abs(cp) in_radius = np.logical_and(min_radius < cp_mag, cp_mag < max_radius) return agg_fun(spectrum[np.logical_and(in_radius)]) # 指定半径ごとの平均を求める def average_radius(spectrum, radius_gap): means = [aggregate_in_angle(np.mean, spectrum, radius, radius + radius_gap) for radius in range(0, 32, radius_gap)] return means main()
下記の角度方向分布を求めるプログラムを元に動軸方向分布のプログラムを作りました.
python スペクトル分布https://teratail.com/questions/157585
グラフは以下のようになりました.角度方向分布は欲しい結果を得る事ができたのですが,動軸方向分布の結果は異なるものでした.
- パワースペクトル分布
- 角度方向分布(指定角度ごとの平均)
- 動軸方向分布(指定半径ごとの平均)
②2DFFTのプログラム(Excel)との比較
得られた結果が正しいかを確認する為,「エクセルで作る地形解析図」にあるExcelファイルで確認しました.
こちらのファイルは数値を入力して,計算すると入力された数値を二次元フーリエ変換し,角度方向の平均,動軸方向の平均が算出されます.
https://husvsakai.cup.com/tikei/data/2D_Fourier_19.zip
角度方向分布は同様の結果を得る事ができ,問題はなかったのですが,動軸方向分布の結果は異なるものでした.
例)動軸方向分布のグラフ
上記のように左肩上がりのグラフを表示したいのですが,pythonの方の動軸方向分布は大きく振動したグラフになっています.
参考にしたもの一覧
- python スペクトル分布
https://teratail.com/questions/157585
- エクセルで作る地形解析図
https://husvsakai.cup.com/tikei/tikei_6.html
- 2DFFTを行った後に求めたいグラフが乗っているレポート
https://www.edu.kobe-u.ac.jp/eng-arch-sled/dat/research/adachi-kigami-1/1.pdf
- csvファイルの数値を読み込み,二次元フーリエ変換を行いたい
まだ回答がついていません
会員登録して回答してみよう