FFT後に正規化、交流成分2倍をしても振幅が改善されませんでした

Question

前提・実現したいこと

csvファイルで，2SIN(2 * PI() * 10^7 * A7)+5SIN(2 * PI() * 3.5*10^7 * A7)のような正弦波を入力し（A列には時間を定義しています），pythonで読み込みをしました．表示・FFT処理を行ったところ，ピーク時の周波数は10MHz，35MHzと意図する周波数になったのですが，振幅の方が思ったようになりませんでした．上の関数は2と5の振幅の正弦波の組み合わせなので，以下のコードのように正規化（ 1/(データ数)を掛ける），交流成分2 倍を行ったら，周波数特性のグラフで10MHzのところに振幅の値2，35MHzのところに振幅の値5のピークが出るはずなのではないでしょうか?
（ちなみに以下のプログラムは超音波パルスエコーの切り出し処理を前提に書いたもののため，時間軸で3か所の切り出しを行い3つそれぞれの周波数特性を求めAverageWaveとして表示させているのでグラフがいくつも出てしまっています．ご了承ください．）
![
＊拡大バージョン

実行プログラム

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# # データの読み込み
##csvファイル用正弦波関数：Amp * SIN(2.0 * PI() * f * t)
df = pd.read_csv('C:/ホットプレート電気炉焼成16コ比較/16csv/C1_00000test.csv', header=4) #ヘッダーは0から数える
t = np.array(df['Time'])
volt = np.array(df['Ampl'])
use_index = np.where(t > 0) #0以上のtを取得
t = t[use_index]
volt = volt[use_index]
src_fs = 1 / (t[1] - t[0]) #サンプリング周波数
src_wav_data_length = t.shape[0] #全波形データ長（〇.shapeで〇の行と列の数が取り出せる（今回は行の数））
src_wav_time_length = src_wav_data_length / src_fs #全波形当たりの時間
aux_line_length = max(abs(volt)) / 2 
all_arr = np.zeros((5, src_wav_data_length)) #平均処理のための行列
N = src_wav_data_length


##切り出し範囲指定
################C1_00000.csv##################
start_No_1 = 1.3E-6 #us
end_No_1 = 2.35E-6 #us
start_No_2 = 2.56E-6
end_No_2 = 3.5E-6
start_No_3 = 3.6E-6
end_No_3 = 4.6E-6

######グラフ作成########
fig1 = plt.figure()
wav_data_length_us = src_wav_data_length
target_arr = np.array([start_No_1, end_No_1, start_No_2, end_No_2, start_No_3, end_No_3])

###################時間軸データ1###################
ax1 = fig1.add_subplot(2, 1, 1)
plt.plot(t * 1000000, volt, linewidth=1)
plt.vlines([start_No_1*1E6, start_No_1*1E6, end_No_1*1E6], -aux_line_length, aux_line_length, color = "blue")
plt.vlines([start_No_2*1E6, end_No_2*1E6], -aux_line_length, aux_line_length, color = "green")
plt.vlines([start_No_3*1E6, end_No_3*1E6], -aux_line_length, aux_line_length, color = "red")
plt.axis("tight") #すべてのデータが見えるように最大・最小を変更する
plt.xlabel("time [us]")
plt.ylabel("amplitude [V]")
plt.grid()

###################時間軸データ2###################
ax2 = fig1.add_subplot(2, 1, 2)
plt.plot(t * 1000000, volt, linewidth=1)
plt.axis("tight") #すべてのデータが見えるように最大・最小を変更する
plt.xlabel("time [us]")
plt.ylabel("amplitude [V]")
plt.grid()
plt.tight_layout()

plt.show()

print(src_fs)

fig, ax = plt.subplots()  # figure と axes を定義
###################周波数軸データ1~3###################
for a in range(0, 5, 2):
    target_gate = np.zeros(wav_data_length_us, dtype='float32')
    target_gate[int(target_arr[a]*src_fs):int(target_arr[a+1]*src_fs)] = 1.0
    target_wav = np.zeros([wav_data_length_us], dtype='float32')
    target_wav = volt * target_gate
    target_wav = np.roll(target_wav, -1*int(target_arr[a]*src_fs))
    
    # 信号の作成
    # src_wav_data_length = t.shape[0] #全波形データ長（〇.shapeで〇の行と列の数が取り出せる（今回は行の数）
    dt = 1 / src_fs
    freq = src_fs # 周波数(10Hz) =>正弦波の周期0.1sec
    t = np.arange(0, N*dt, dt) # 時間軸
    f = target_wav # 信号（周波数10、振幅1の正弦波）
    F = np.fft.fft(f) 
    F_abs = np.abs(F)# FFTの複素数結果を絶対に変換
    M = int(target_arr[a+1]*src_fs) - int(target_arr[a]*src_fs) + 1
    # 振幅をもとの信号に揃える
    F_abs_amp = F_abs / M * 2 # 交流成分はデータ数で割って2倍
    F_abs_amp[0] = F_abs_amp[0] / 2 # 直流成分は2倍不要

    fq = np.linspace(0, 1.0/dt, N) # 周波数軸　linspace(開始,終了,分割数)


    if a == 0:
        ax.set_xlabel("freqency[Hz]")
        ax.set_xlim(0, 10E6)
        ax.set_ylabel("amplitude[arb.unit]")
        ax.plot(fq[:int(N/2)+1], F_abs_amp[:int(N/2)+1], label='1st-Wave_FrequencyResponse', color="blue") # ナイキスト定数まで表示
        ax.legend()
        ax.grid()
        ax.axis("tight") #すべてのデータが見えるように最大・最小を変更する
        ax.set_xlim(0, 4.0E7)


    elif a == 2:
        ax.plot(fq[:int(N/2)+1], F_abs_amp[:int(N/2)+1], label='2nd-Wave_FrequencyResponse', color="green") # ナイキスト定数まで表示
        ax.legend()
        
    else :
        ax.plot(fq[:int(N/2)+1], F_abs_amp[:int(N/2)+1], label='3rd-Wave_FrequencyResponse', color="red") # ナイキスト定数まで表示
        ax.legend()
        
    for i in range(N):
        all_arr[a][i] = all_arr[a][i] + F_abs_amp[i]


###################周波数軸データ平均###################
all_ave = (all_arr[0][0:N] + all_arr[2][0:N] + all_arr[4][0:N]) / 3 
ax.plot(fq[:int(N/2)+1], all_ave[:int(N/2)+1], label='AverageWave', color = "orange") # ナイキスト定数まで表示
ax.legend()

plt.show()

該当のソースコード

 # 振幅をもとの信号に揃える
    F_abs_amp = F_abs / M * 2 # 交流成分はデータ数で割って2倍
    F_abs_amp[0] = F_abs_amp[0] / 2 # 直流成分（今回は扱わないけど）は2倍不要

試したこと

以下のリンクを始めとする分かりやすいページで，FFT処理した後の正規化や交流成分2 倍を調べました．ですが，上のプログラムで，周波数特性のグラフで10MHzのところに振幅の値2（こちらの方は解決いたしました，ありがとうございました），35MHzのところに振幅の値5が出ませんでした．webにあげられているものを組み合わさせていただきながらいろいろ書き替えたため非常に分かりにくいプログラムとなっていますが，上記のプログラムの問題点が分かられる方ご回答よろしくお願いします．
https://org-technology.com/posts/fft-03.html

以下に単体の正弦波の周波数を変えてプログラムに通した結果を示します．（振幅は2で試しています）
．10MHzの場合の時間軸グラフと周波数軸グラフ

．15MHzの場合の時間軸グラフと周波数軸グラフ

．20MHzの場合の時間軸グラフと周波数軸グラフ

．30MHzの場合の時間軸グラフと周波数軸グラフ

．31MHzの場合の時間軸グラフと周波数軸グラフ

．32MHzの場合の時間軸グラフと周波数軸グラフ

．33MHzの場合の時間軸グラフと周波数軸グラフ

．34MHzの場合の時間軸グラフと周波数軸グラフ

．35MHzの場合の時間軸グラフと周波数軸グラフ

csvファイルでは以下のように正弦波を作成しています．（csvファイルは，保存してもいったんファイルを閉じてしまうと数値しか保存されていておらず定義した数式が保存されなかったので（保存方法が分かりませんでした），以下にはファイルではなくスクショ画面を載せておきます）印をつけた部分の値を，上のグラフで定義した周波数に変更してそれぞれプログラムに通してみました．

・実際に処理したいデータ

Comments

[jbpb0]

> F = np.fft.fft(f)
F_abs = np.abs(F)# FFTの複素数結果を絶対に変換
F_abs_amp = F_abs / N * 2 # 交流成分はデータ数で割って2倍

「f」のデータ数は「N」ですか？

> src_wav_data_length = t.shape[0] #全波形データ長（〇.shapeで〇の行と列の数が取り出せる（今回は行の数））

> N = src_wav_data_length

「N」は元データ全体の長さですよね

コードよく読んでないのですが、「f」がもし元データ全体ではなく、そこから切り出した一部なら、「N」ではなく「f」のデータ数で正規化しないといけないですけど、そこ大丈夫でしょうか？

[ma2_718]

すみません，その通りですm(_ _)m
ありがとうございました！
Nをfのデータ数に変更して修正しました<(_ _)>

変更したところ10MHzのピーク値は指定通りの2となりました！本当にありがとうございました．ですが，35MHzのピーク値（振幅）が作成した正弦波の振幅の5ではなく4になってしまいました．作成した波形の振幅を書き換えていろいろ試したのですが，35MHzのピーク値がどうしても指定の値より1小さくなってしまいました．何度も申し訳ありません，お暇があられるとき教えて頂けると幸いです．

[jbpb0]

二つの周波数の合成ではなく、周波数を一つだけでデータを作ったら、どうなりますか？
・元データが10MHz振幅2だけの場合
・元データが35MHz振幅2だけの場合
のそれぞれで、f特はどうなりますか？
どちらもf特のピークは2になりますか？

もし、元データの振幅は同じなのに、周波数によってf特のピークが異なるのなら、10MHzから周波数を少しずつ変えて、10, 15, 20, 25, 30, 35MHzみたいに徐々に変えたら、どうなりますか？

[jbpb0]

あと、質問の横軸timeのグラフは、横軸が0〜1あたりと、それよりも右で様子が違っているように見えるのですが、実際そうなのでしょうか？
最初から最後まで

> 2*SIN(2 * PI() * 10^7 * A7)+5*SIN(2 * PI() * 3.5*10^7 * A7)のような正弦波を入力

なら、そうはならないはずですが

[jbpb0]

もう一つ、サンプリング周波数はいくつで、「f」のデータ数はいくつでしょうか？

> A列には時間を定義しています

の時間(時刻)の間隔は、全データで一定ですよね？

> 時間軸で3か所の切り出しを行い3つそれぞれの周波数特性を求めAverageWaveとして表示

の3ヵ所の切り出しで、「f」の長さは全て同じですよね？

[ma2_718]

ご指摘ありがとうございます．上部の「試したこと」の続きに周波数を少しずつ変えてみた時の変化を追加しています．おっしゃる通り，横軸timeのグラフで横軸1以前と1以降とで様子が違っており，おそらくそれが原因ではないかと思われました．

[jbpb0]

> 上部の「試したこと」の続きに周波数を少しずつ変えてみた時の変化を追加しています．

が見当たりません

> 横軸timeのグラフで横軸1以前と1以降とで様子が違っており

元データが

> 上の関数は2と5の振幅の正弦波の組み合わせ

となってないなら、予想通りにはなりませんよね
元データの素性を正確に把握するところからですね

あと、私の一つ前のコメントで質問してる下記も、差し支えなければ教えてください
・サンプリング周波数はいくつか？
・「f」のデータ数はいくつか？
・サンプリング周波数は、データ全体で同一か？
・3ヵ所の切り出しで、「f」のデータ数は全て同じか？

[jbpb0]

> 上部の「試したこと」の続きに周波数を少しずつ変えてみた時の変化を追加しています．

が見えました
横軸1以降が、意味不明な形してますね

エクセルで正しいCSVが作れないなら、元データもPythonで作ったらいかがですか？

[ma2_718]

・サンプリング周波数はいくつか？
・・・プログラムの「##データの読み込み」の欄にあるsrc_fs = 1 / (t[1] - t[0])で定義したところです．print(src_fs)と入れてみたら2500000000.96875と返ってきました．約2500MHz．．って大きすぎる気がしました(..)

・「f」のデータ数はいくつか？
・サンプリング周波数は、データ全体で同一か？
・3ヵ所の切り出しで、「f」のデータ数は全て同じか？
・・・「f」のデータ数に関しては，target_wav すなわち volt * target_gateのデータ数から来ています．target_gate は target_gate[int(target_arr[a]*src_fs):int(target_arr[a+1]*src_fs)] = 1.0 としているので，プログラムの「##切り出し範囲指定」の部分で指定した切り出しの値で決まります．たいへん見えにくいですが，時間軸グラフの上の段にはすべて順に青緑赤の切り出しゲートが表示されています．この切り出し範囲については，実際に実験データを処理する場合にパルスの位置（「試したこと」に追加した「・実際に処理したいデータ」）が変わってくるので得られた実験データに応じてその都度変更するようにプログラムを書いています．（ですので，3ヵ所の切り出しで「f」のデータ数は全て同じ，とは必ずしもならないです．）例えば，for文でaが0の時には青色部分で切り出しているので，
start_No_1 = 1.3E-6 #us
end_No_1 = 2.35E-6 #us
で指定した時間範囲内のデータだけが切り出されます．（start_No_1からend_No_1 までは値1，それ以外は0となるtarget_gateを元信号に乗することでtarget_wav（切り出し波形）を生成しています．）そして，その範囲のデータ'数'を考えるとサンプリング周期であるsrc_fsを，切り出し時間の両端の時間（start_No_1とend_No_1）に乗して，ベクトルtarget_gateのインデックスとして用いている数値の差を求めた（個数なので+1）M
M = int(target_arr[a+1]*src_fs) - int(target_arr[a]*src_fs) + 1
が青色部分の切り出しデータ数となります．for文内でaが2，4と増えていき，緑色切り出し，赤色きりだしとなっても同様の操作です．
また，サンプリング周波数に関しては，オシロスコープでとった時間間隔（(t[1] - t[0])）の逆数としていてデータ処理するうえで全体で同一にしています．

アドバイスいただいた通りcsvファイルを変更してみたら，最初予想していた通りにはならないはずだと気づくことができました．csvファイルの方もいじりながら元データ素性を正確に把握できるように頑張ります！本当にありがとうございました！

[jbpb0]

> 3ヵ所の切り出しで「f」のデータ数は全て同じ，とは必ずしもならないです．

それだと、FFTの結果のf特の周波数が3回で一致しないかもしれません
周波数が違えば、単純に平均できません

FFTの結果は周波数が等間隔に並びますが、その周波数の間隔はFFTしたデータの全長の逆数です
たとえば全体で10秒のデータの場合は1/10=0.1Hz間隔なので、
0, 0.1, 0.2, 0.3…Hz
となり、全体で11秒の場合は1/11=0.091Hz間隔なので、
0, 0.091, 0.182, 0.273…Hz
となり、異なる周波数の結果になります

コードでは三つのf特のプロットで共通の「fq」を用いるようですが、三つそれぞれのFFTデータ「f」の秒の長さで異なる周波数になるはずで、
fq = np.linspace(0, 1.0/dt, N) # 周波数軸　linspace(開始,終了,分割数)
は、(正規化のデータ数と同様に)「N」ではなく「f」のデータ数で計算しないといけないと思います
もちろん、3回のそれぞれで計算し直す必要があります

そうしたら、「f」のデータ数が違う場合は周波数が一致しないので、そのままではf特が平均できません
周波数を補間するとかして周波数を一致させたら平均できますが、そんな怪しいことするよりも「f」の長さを3回で揃える方がいいと思います (絶対に平均が必要なら)

[ma2_718]

確かに，周波数特性を出す際切り出しデータの処理なので「f」のデータ数で計算する必要がありました．ご指摘のfq の部分とそれに付随する部分のプログラムを書き直してみたいと思います！（平均処理はあきらめるかもしれませんが．．．）
ありがとうございましたm(__)m

[jbpb0]

平均が必須ではなく、3回の結果を同一のグラフで比較できたらよいのであれば、3回それぞれで周波数の計算を別個に正しく行ってプロットすれば、もちろんOKです

[jbpb0]

訂正します

コードを詳しく読んだら、「f」は必要ないところは0にしてるだけで、サイズは3回とも全データのサイズと同じですね
その場合は、f特の周波数は3回とも同じになり、全データ長さで算出すればいいので、現状の質問のコードの
fq = np.linspace(0, 1.0/dt, N) # 周波数軸　linspace(開始,終了,分割数)
で大丈夫です
失礼しました

なお、f特の振幅の正規化は、実際にデータが入ってるところの長さでやるので、現状の質問のコードの
M = int(target_arr[a+1]*src_fs) - int(target_arr[a]*src_fs) + 1
F_abs_amp = F_abs / M * 2 # 交流成分はデータ数で割って2倍
で合ってます

[ma2_718]

すみません，私の説明が良くなかったです．「f」の”サイズ”は3回とも全データサイズと同じですね．
合っているようで良かったです．ありがとうございました．