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回とも全データサイズと同じですね．
合っているようで良かったです．ありがとうございました．

Answer 1

データをPythonで作るようにしました
データの全長さは適当ですが

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]
##
dt0 = 1 / 2500000000  # データ間隔
at0 = 7.314 * 10**-6  # 全データ長
f10 = 10 * 10**6  # 周波数1
f20 = 35 * 10**6  # 周波数2
a10 = 2  # 周波数1の振幅
a20 = 5  # 周波数2の振幅
##
t = np.arange(0, at0, dt0)  # arange(開始, 終了の次, 間隔)
volt = a10 * np.sin(2 * np.pi * f10 * t) + a20 * np.sin(2 * np.pi * f20 * t)
##
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特振幅ピークは、ほぼ狙い通りです

CSVファイルのデータが正しくできてれば、こうなるはず

次に、「f」を、必要な部分のみ切り出して他は捨てる(0埋めではなく)ようにコードを変えます
その場合は、3回の「f」の長さが違うと周波数が合わないので、平均はしません

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]
##
dt0 = 1 / 2500000000  # データ間隔
at0 = 7.314 * 10**-6  # 全データ長
f10 = 10 * 10**6  # 周波数1
f20 = 35 * 10**6  # 周波数2
a10 = 2  # 周波数1の振幅
a20 = 5  # 周波数2の振幅
##
t = np.arange(0, at0, dt0)  # arange(開始, 終了の次, 間隔)
volt = a10 * np.sin(2 * np.pi * f10 * t) + a20 * np.sin(2 * np.pi * f20 * t)
##
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) # 時間軸
    #M = int(target_arr[a+1]*src_fs) - int(target_arr[a]*src_fs) + 1
    M = int(target_arr[a+1]*src_fs) - int(target_arr[a]*src_fs)
    #f = target_wav # 信号（周波数10、振幅1の正弦波）
    f = target_wav[:M]
    F = np.fft.fft(f) 
    F_abs = np.abs(F)# FFTの複素数結果を絶対に変換
    # 振幅をもとの信号に揃える
    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(開始,終了,分割数)
    fq = np.linspace(0, 1.0/dt, M) # 周波数軸　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.plot(fq[:int(M/2)+1], F_abs_amp[:int(M/2)+1], label='1st-Wave_FrequencyResponse', color="blue", marker="o") # ナイキスト定数まで表示
        ax.legend()
        ax.grid()
        ax.axis("tight") #すべてのデータが見えるように最大・最小を変更する
        ax.set_xlim(0, 4.0E7)
        ax.set_ylim(-0.2, 5.2)


    elif a == 2:
        #ax.plot(fq[:int(N/2)+1], F_abs_amp[:int(N/2)+1], label='2nd-Wave_FrequencyResponse', color="green") # ナイキスト定数まで表示
        ax.plot(fq[:int(M/2)+1], F_abs_amp[:int(M/2)+1], label='2nd-Wave_FrequencyResponse', color="green", marker="o") # ナイキスト定数まで表示
        ax.legend()

    else :
        #ax.plot(fq[:int(N/2)+1], F_abs_amp[:int(N/2)+1], label='3rd-Wave_FrequencyResponse', color="red") # ナイキスト定数まで表示
        ax.plot(fq[:int(M/2)+1], F_abs_amp[:int(M/2)+1], label='3rd-Wave_FrequencyResponse', color="red", marker="o") # ナイキスト定数まで表示
        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()

0埋めがなくなり、データがあるところと0埋めしてるところの段差が無くなるので、f特のたくさんあった細かい振動は無くなります
3回の切り出しの内で、1, 2回目はデータの両端が不連続なので、f特振幅ピークが狙いと合わなくなります
(35MHzだけじゃなく、10MHzも)

今回の質問のデータの場合は、周波数が10, 35MHzと整数MHzしかないので、「f」の切り出し長さを

start_No_1 = 1.3E-6 #us
#end_No_1 = 2.35E-6 #us
end_No_1 = 2.3E-6 #us
#start_No_2 = 2.56E-6
start_No_2 = 2.5E-6
end_No_2 = 3.5E-6
start_No_3 = 3.6E-6
end_No_3 = 4.6E-6

のようにして1μsecの整数倍に合わせればデータの両端が連続になるので、f特の10, 35MHz以外の振幅は消え、振幅ピークはピッタリ狙い通りになります

あくまでも、今回の質問のデータのようなきれいなデータの場合はそうできる、というだけですが、ご参考までに

Comments

[ma2_718]

ありがとうございます！周波数や振幅を変えてみても適当な値が出てきました．
pythonで直接定義する方が早かったですね(><)
CSVファイルデータの方も今後時間があれば検討していきたいと思います！
本当にお世話になりました！ありがとうございましたm(__)m

[ma2_718]

ありがとうございます！助かります(´▽｀)

Answer 2

質問者やコメントされている方は納得がいかないかもしれませんが，
結果は私が想定した通りの結果となっています。特に間違っているところはありません。

実際のところ，これはプログラムの問題ではなくて数学の問題です。
想定通りにならない理由の一つは，波形の一部分のみを矩形窓で切り出してFFTを行っているところです。

DFT(離散フーリエ変換)は，演算に用いる有限時間の波形が無限に繰り返すとし
てフーリエ級数の係数を求めます。

したがって矩形窓の0と1のところで不連続点が発生し，単一のはずの周波数スペクトルに広がりが
みられるようになります。理想的には，単一周波数の場合周波数スペクトルは1本になるはずですが，
グラフを見るとsinc関数[(sin x)/x]の絶対値のような特性になっているのはこれが理由です。

ただ，元波形がsin波形なので，矩形窓で1にしている時間幅が周期と同じ時はsin波形の両端の
値はほぼ0なので，他の周波数に洩れる成分は小さいです。これが10MHzの成分の
中心周波数スペクトルがほぼ想定通りの値になっている理由です。

35MHz成分の場合は，おそらく矩形窓の1にしている両端で振幅が大きくなっていて，
不連続点の影響が大きく，35MHz以外のスペクトルに分散しているため中心周波数スペクトルが
小さく表示されているのだと思います。

じゃあどうすればいいねんというところですが，
35MHzの周期の整数倍の時間の波形を切り出して(どれぐらい切り出せばよいかは自分で計算して下さい。)，切り出したデータそのままFFTをすると洩れスペクトルはなくなります。
ソースリストでいうと，Fを代入しているところで元データのvoltを直接範囲を指定して
そのままFFTすればいいと思います。次のような感じですかね。

F = np.fft.fft(volt[int(target_arr[a]*src_fs):int(target_arr[a+1]*src_fs)])

グラフにする時の周波数分解能(刻み周波数)は切り出したデータ数に従って可変となりますが，
「ΔF=サンプリング周波数/FFTに突っ込むデータ数」で求まるので，それを利用すればよいでしょう。

なお，実際の分析ではサンプリング周波数や演算窓の関係から漏れスペクトルの影響は避けられない
ことも多いので，大きさの評価は中心周波数の前後数本のスペクトルのピタゴラス和をとって行ったりもします。

Comments

[ma2_718]

なるほど！切出しタイミングの影響で振幅が小さくなっていたのですね．．．
ご指摘いただいたように，voltの範囲を直接指定してみて試してみようと思います．
ありがとうございましたm(__)m