前提・実現したいこと
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ファイルは,保存してもいったんファイルを閉じてしまうと数値しか保存されていておらず定義した数式が保存されなかったので(保存方法が分かりませんでした),以下にはファイルではなくスクショ画面を載せておきます)印をつけた部分の値を,上のグラフで定義した周波数に変更してそれぞれプログラムに通してみました.
・実際に処理したいデータ
まだ回答がついていません
会員登録して回答してみよう