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

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

新規登録して質問してみよう
ただいま回答率
85.45%
Python 3.x

Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。

Q&A

解決済

2回答

1391閲覧

pythonで逆フーリエ変換

Pinkman1112

総合スコア4

Python 3.x

Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。

0グッド

0クリップ

投稿2022/10/11 02:25

編集2022/10/11 02:31

前提

信号処理はほとんど初心者です。
pythonでノイズ周波数スペクトルから逆fftでノイズ信号を得たいと考えております。

実現したいこと

・ジョンソンノイズスペクトルからノイズ信号を取得
・シミュレーションで得た周波数信号のデータが実数列で単純に逆フーリエ変換しても正しい信号が得られない(自分では周波数信号を取得していなく、スペクトルデータしかない)
・ネットで調べていると逆FFTはFFTした虚数データから逆FFTを行なっている例はよく見られるのですが、実数(パワースペクトル?)だけのデータから逆フーリエ変換は可能でしょうか?

発生している問題・エラーメッセージ

上が周波数スペクトル、下が逆フーリエ変換した信号
イメージ説明

イメージ説明

該当のソースコード

import numpy as np import os import glob from natsort import natsorted from scipy.fftpack import fft,ifft import matplotlib.pyplot as plt list = natsorted(glob.glob("pulse_sim_*.txt")) noise_spe_sim = np.loadtxt("noise_spe.txt") #ノイズスペクトルのデータ #------parameter取得------ with open("input.txt","r") as f: para = np.loadtxt(f,comments="#",delimiter=",") samples = para[16] rate = para[17] eta =105 time_max =500.0e-3 N_time_max = 10000 fre_max = 3.0e+6 N_fre_max = 2**18+1 ptfn_Flink = 0.5 #------------------------ time = np.linspace(0.0,time_max,N_time_max) fq = np.linspace(0.0,fre_max,N_fre_max) plt.plot(fq[:int(N_fre_max/2)+1],noise_spe_sim[:int(N_fre_max/2)+1],linestyle = '-',linewidth = 0.7) plt.xscale('log') plt.yscale('log') plt.xlabel('Frequency [Hz]',fontsize = 16) plt.ylabel('Noise [fA/rtHz]',fontsize = 16) plt.show() plt.cla() r = ifft(noise_spe_sim,len(time)) plt.plot(time,r) plt.xlabel('time[s]',fontsize = 16) plt.show()

補足情報(FW/ツールのバージョンなど)

noise_spe.txtの中身

3.915279528816821193e-11
3.997824114377672742e-11
4.080145862759802740e-11
4.162145186680246524e-11
4.243723294845935490e-11
4.324782879527359946e-11
4.405228809852697713e-11
4.484968810654406152e-11
4.563914107104591467e-11

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

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

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

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

PondVillege

2022/10/11 02:54

> 逆FFTはFFTした虚数データから逆FFTを行なっている のところ,正しくは「複素数データから」ではないでしょうか?
jbpb0

2022/10/11 03:01

実数の周波数スペクトルが、本来の複素数のスペクトルの絶対値(log変換とかしてない)であり、周波数が等間隔(周波数のlogの等間隔ではなく、普通の周波数で等間隔)で並んでるのなら、 https://watlab-blog.com/2019/04/21/python-fft/#350363627565306125111252112540125221253112464123922160827874259683660012395123881235612390 の「補足:ミラーリングと周波数軸について」の「ミラーリングとナイキスト周波数」の、横軸が周波数のグラフの「ミラーリングされている」と書かれてる、ナイキスト周波数よりも高い周波数のスペクトルデータを追加すれば、逆フーリエ変換の計算はできます 上記のグラフを見ると何となく分かるかもしれませんが、元のスペクトルデータから0Hzのデータのみ除いたものを、順番を逆にして、元のスペクトルデータ(の後ろ)に追加します ただし、スペクトルデータが複素数ではなく実数なので、それから逆フーリエ変換で求まった信号は、全周波数の信号の位相が全て合ってるという状態になります > 逆fftでノイズ信号を得たい の「ノイズ信号」は、全周波数の位相が合っててもいいのでしょうか?
jbpb0

2022/10/11 03:29

どんな「ノイズ信号」を作りたいのかがよく分かりませんが、周波数特性の振幅が質問のグラフにあるような所望のもので、位相がランダムなものが作りたいのなら、実数の振幅データから、乱数でランダムな各周波数の位相をでっち上げて複素数にしてしまって、それを逆フーリエ変換するとか ただしその場合は、私の一つ前のコメントに書いたようなナイキスト周波数よりも高い周波数のデータを追加する際に、位相の符号を逆にする必要があります たとえば、ある周波数で元データの実数から複素数化したデータが「1.10841641+3.28462627j」だった場合は、追加するデータは「1.10841641-3.28462627j」になります(虚部の±が逆)
Pinkman1112

2022/10/11 04:07

> どんな「ノイズ信号」を作りたいのかがよく分かりませんが、周波数特性の振幅が質問のグラフにあるような所望のもので、位相がランダムなものが作りたいのなら、実数の振幅データから、乱数でランダムな各周波数の位相をでっち上げて複素数にしてしまって、それを逆フーリエ変換するとか やりたいことはおっしゃる通りで位相がランダムなノイズを作りたいと考えております。正直あまり理解が追いついていないのですが、乱数で位相をでっち上げるというのは0~1の適当な乱数でいいということでしょうか?
jbpb0

2022/10/11 04:33 編集

> 乱数で位相をでっち上げるというのは0~1の適当な乱数でいいということでしょうか? http://www.ee.t-kougei.ac.jp/tuushin/lecture/math1/htdocs/complex/complexPlane.html#complexPlane の「3. 複素平面 (complex plane)」で、「r」が振幅(元の実数データ)、「θ」が位相です 「θ」の範囲は、単位がradianの場合は0〜2πです (元実数データそのままの)振幅「r」と、0〜2πの範囲の乱数から生成した「θ」で、上記Webページの「x」と「y」を算出して、複素数「z」にします
Pinkman1112

2022/10/11 05:53

できました!ありがとうございます!
guest

回答2

0

自己解決

上がミラーリングした周波数スペクトル、下が逆フーリエ変換したノイズ

非常に勉強になりました。

イメージ説明

イメージ説明

投稿2022/10/11 05:54

Pinkman1112

総合スコア4

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

0

不可能です。
パワースペクトルは振幅の情報はありますが位相の情報が欠落しています。

例えば元の波形がsin(wt)cos(wt)では位相が異なりますが、
パワースペクトル上ではどちらも各周波数wにピークを持ち区別できません。

投稿2022/10/11 02:46

ozwk

総合スコア13532

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

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

ただいまの回答率
85.45%

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

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

質問する

関連した質問