２DFFTから動径方向分布を算出したい

Question

やりたいこと

①csvファイルの数値を読み込んで，二次元フーリエ変換(2DFFT)を行い，計算結果をcsvファイルで出力したいです．

②①で得られた結果から動径方向分布 p(r) と角度方向分布 q(θ)のスペクトル分布を算出したいです．

困っている事

②の２DFFTから動径分布を計算を算出する部分．

CSVファイルの中身

CSVファイルの中身はこのようになっています．

前回の質問のcsvファイルとは異なります．https://teratail.com/questions/ds77z65aj8c652
余分な部分を削除し，データの縦×横を64×64に揃えました．

試したこと

①２DFFTのプログラム(python)の改造

pythonで書かれた画像の２DFFTのプログラム　https://watlab-blog.com/2020/03/22/2d-fft/

画像を読み込む部分をcsvに置き換え, 角度方向分布，動軸方向分布のグラフを算出するプログラムを作りました．

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

グラフは以下のようになりました．角度方向分布は欲しい結果を得る事ができたのですが，動軸方向分布の結果は異なるものでした．

• パワースペクトル分布

• 角度方向分布（指定角度ごとの平均）

• 動軸方向分布（指定半径ごとの平均）

②２DFFTのプログラム(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

• ２DFFTを行った後に求めたいグラフが乗っているレポート

　https://www.edu.kobe-u.ac.jp/eng-arch-sled/dat/research/adachi-kigami-1/1.pdf

• csvファイルの数値を読み込み，二次元フーリエ変換を行いたい

　https://teratail.com/questions/ds77z65aj8c652

Answer 1

 # スペクトルを与えられた半径の範囲内だけ集計
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)
    in_radius = np.logical_and(min_radius <= cp_mag, cp_mag < max_radius)
    #return agg_fun(spectrum[np.logical_and(in_radius)])
    return agg_fun(spectrum[in_radius])

 # 指定半径ごとの平均を求める
def average_radius(spectrum, radius_gap):
    #means = [aggregate_in_angle(np.mean, spectrum, radius, radius + radius_gap)
    means = [aggregate_in_radius(np.mean, spectrum, radius, radius + radius_gap)
             for radius in range(0, 32, radius_gap)]
    return means 

で、それっぽい値が計算されると思います
細部は、自分で直してください

　
【追記】

余分な部分を削除し，データの縦×横を64×64に揃えました．

csvファイルの数値を読み込み，二次元フーリエ変換を行いたい
の回答にも書きましたけど、CSVファイルの1行目をヘッダーとして扱われないようにしないといけません

 df = pd.read_csv('pivdata.csv')

↓ 変更

 df = pd.read_csv('pivdata.csv', header=None)

Answer 2

average_radius内でaggregate_in_radiusが使われていませんがこれでいいのでしょうか?