質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

新規登録して質問してみよう
ただいま回答率
87.20%
CSV

CSV(Comma-Separated Values)はコンマで区切られた明白なテキスト値のリストです。もしくは、そのフォーマットでひとつ以上のリストを含むファイルを指します。

関数

関数(ファンクション・メソッド・サブルーチンとも呼ばれる)は、はプログラムのコードの一部であり、ある特定のタスクを処理するように設計されたものです。

Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

解決済

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

ma2_718
ma2_718

総合スコア2

CSV

CSV(Comma-Separated Values)はコンマで区切られた明白なテキスト値のリストです。もしくは、そのフォーマットでひとつ以上のリストを含むファイルを指します。

関数

関数(ファンクション・メソッド・サブルーチンとも呼ばれる)は、はプログラムのコードの一部であり、ある特定のタスクを処理するように設計されたものです。

Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

2回答

0評価

0クリップ

1244閲覧

投稿2021/07/20 06:20

編集2021/07/20 16:25

前提・実現したいこと

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ファイルは,保存してもいったんファイルを閉じてしまうと数値しか保存されていておらず定義した数式が保存されなかったので(保存方法が分かりませんでした),以下にはファイルではなくスクショ画面を載せておきます)印をつけた部分の値を,上のグラフで定義した周波数に変更してそれぞれプログラムに通してみました.
イメージ説明

・実際に処理したいデータ
イメージ説明
イメージ説明

良い質問の評価を上げる

以下のような質問は評価を上げましょう

  • 質問内容が明確
  • 自分も答えを知りたい
  • 質問者以外のユーザにも役立つ

評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

気になる質問をクリップする

クリップした質問は、後からいつでもマイページで確認できます。

またクリップした質問に回答があった際、通知やメールを受け取ることができます。

teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。

  • プログラミングに関係のない質問
  • やってほしいことだけを記載した丸投げの質問
  • 問題・課題が含まれていない質問
  • 意図的に内容が抹消された質問
  • 過去に投稿した質問と同じ内容の質問
  • 広告と受け取られるような投稿

評価を下げると、トップページの「アクティブ」「注目」タブのフィードに表示されにくくなります。

jbpb0

2021/07/20 06:57 編集

> 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

2021/07/20 08:49

すみません,その通りですm(_ _)m ありがとうございました! Nをfのデータ数に変更して修正しました<(_ _)> 変更したところ10MHzのピーク値は指定通りの2となりました!本当にありがとうございました.ですが,35MHzのピーク値(振幅)が作成した正弦波の振幅の5ではなく4になってしまいました.作成した波形の振幅を書き換えていろいろ試したのですが,35MHzのピーク値がどうしても指定の値より1小さくなってしまいました.何度も申し訳ありません,お暇があられるとき教えて頂けると幸いです.
jbpb0

2021/07/20 11:36 編集

二つの周波数の合成ではなく、周波数を一つだけでデータを作ったら、どうなりますか? ・元データが10MHz振幅2だけの場合 ・元データが35MHz振幅2だけの場合 のそれぞれで、f特はどうなりますか? どちらもf特のピークは2になりますか? もし、元データの振幅は同じなのに、周波数によってf特のピークが異なるのなら、10MHzから周波数を少しずつ変えて、10, 15, 20, 25, 30, 35MHzみたいに徐々に変えたら、どうなりますか?
jbpb0

2021/07/20 10:06 編集

あと、質問の横軸timeのグラフは、横軸が0〜1あたりと、それよりも右で様子が違っているように見えるのですが、実際そうなのでしょうか? 最初から最後まで > 2*SIN(2 * PI() * 10^7 * A7)+5*SIN(2 * PI() * 3.5*10^7 * A7)のような正弦波を入力 なら、そうはならないはずですが
jbpb0

2021/07/20 10:13 編集

もう一つ、サンプリング周波数はいくつで、「f」のデータ数はいくつでしょうか? > A列には時間を定義しています の時間(時刻)の間隔は、全データで一定ですよね? > 時間軸で3か所の切り出しを行い3つそれぞれの周波数特性を求めAverageWaveとして表示 の3ヵ所の切り出しで、「f」の長さは全て同じですよね?
ma2_718

2021/07/20 14:38

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

2021/07/20 14:53

> 上部の「試したこと」の続きに周波数を少しずつ変えてみた時の変化を追加しています. が見当たりません > 横軸timeのグラフで横軸1以前と1以降とで様子が違っており 元データが > 上の関数は2と5の振幅の正弦波の組み合わせ となってないなら、予想通りにはなりませんよね 元データの素性を正確に把握するところからですね あと、私の一つ前のコメントで質問してる下記も、差し支えなければ教えてください ・サンプリング周波数はいくつか? ・「f」のデータ数はいくつか? ・サンプリング周波数は、データ全体で同一か? ・3ヵ所の切り出しで、「f」のデータ数は全て同じか?
jbpb0

2021/07/20 15:28

> 上部の「試したこと」の続きに周波数を少しずつ変えてみた時の変化を追加しています. が見えました 横軸1以降が、意味不明な形してますね エクセルで正しいCSVが作れないなら、元データもPythonで作ったらいかがですか?
ma2_718

2021/07/20 16:21

・サンプリング周波数はいくつか? ・・・プログラムの「##データの読み込み」の欄にある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

2021/07/20 22:50

> 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

2021/07/21 01:51

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

2021/07/21 03:47

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

2021/07/21 04:56 編集

訂正します コードを詳しく読んだら、「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

2021/07/21 06:01

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

まだ回答がついていません

会員登録して回答してみよう

15分調べてもわからないことは
teratailで質問しよう!

ただいまの回答率
87.20%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

関連した質問

同じタグがついた質問を見る

CSV

CSV(Comma-Separated Values)はコンマで区切られた明白なテキスト値のリストです。もしくは、そのフォーマットでひとつ以上のリストを含むファイルを指します。

関数

関数(ファンクション・メソッド・サブルーチンとも呼ばれる)は、はプログラムのコードの一部であり、ある特定のタスクを処理するように設計されたものです。

Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。