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

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

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

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

Q&A

解決済

2回答

5716閲覧

フーリエ変換でグラフの傾きを除去する方法

NCC1701

総合スコア1680

Python

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

0グッド

1クリップ

投稿2019/08/04 09:44

前提・実現したいこと

傾いているグラフを水平にする方法を、フーリエ変換を利用してやろうとしてます(直流成分を除く?)が、フーリエ変換の理解が浅いようで、うまくいきません。

図1

例えば上図1で青い信号が得られたとします。これをオレンジのグラフのように修正したいのです。

該当のソースコード

python

1import numpy as np 2import matplotlib.pyplot as plt 3%matplotlib inline 4from scipy import signal 5 6N = 64 # サンプル数 7dt = 0.01 # サンプリング周期 8freq = 10 # 周波数 9amp = 1 # 振幅 10t = np.arange(0, N*dt, dt) # 時間軸 11f = amp * np.sin(2*np.pi * freq*t) + 0.5 * freq*t # 実験用の正弦波 12f0 = amp * np.sin(2*np.pi * freq*t) # 目標とする正弦波 13 14F = np.fft.fft(f) # 高速フーリエ変換(FFT) 15F_abs = np.abs(F) 16F_abs_amp = F_abs / N * 2 17F_abs_amp[0] = F_abs_amp[0] / 2 18 19angle = F_abs[0]/N/N*20 # 傾きを取得しているつもり 20 21fq = np.linspace(0, 1.0/dt, N) 22fig = plt.figure(figsize=(12, 4)) 23ax2 = fig.add_subplot(121) ## 信号のグラフ 24plt.xlabel('time(sec)') 25plt.ylabel('amp') 26plt.plot(t, f) 27plt.plot(t, f0) 28ax2 = fig.add_subplot(122) ## FFTのグラフ(周波数軸) 29plt.xlabel('freqency(Hz)') 30plt.ylabel('amp') 31plt.plot(fq[:int(N/2)+1], F_abs_amp[:int(N/2)+1]) 32plt.show()

試したこと

python

1F_ifft = np.fft.ifft(F) # IFFT 2F_ifft_real = F_ifft.real * 2 # 振幅を元スケールに戻している 3 4F_minus = F_ifft_real - angle*t # ここで傾きを除去しているつもり 5 6plt.plot(t, F_ifft_real) 7plt.plot(t, F_minus, c="r") 8plt.show()

上記のようにフーリエ変換で直流成分をangleとして取り出し、元の信号から差し引いてやれば傾きが直ると思ったのですが、うまく行きませんでした。下図のようにほんの少ししか修正されません。
図2

python

1F_ifft = np.fft.ifft(F[1:]) 2F_ifft_real = F_ifft.real * 2 3plt.plot(t[1:], F_ifft_real)

また、上述のようにFFT配列の0番目を除去して試したものの下図のようにこれもうまくいきません
図3

傾きが直っている気配すらなく、どの辺がまずいのかもわからなくなってます。どうぞよろしく。

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

python 3(Colaboratory)
参考にしたサイト

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

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

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

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

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

guest

回答2

0

ベストアンサー

fftでDC近傍(上のグラフの傾きはとても周波数の遅い成分と捉えます)を落とす必要があります。
この周波数成分をを除きたいのならコードのF[0] から少しの範囲を0にして、逆fftすれば行けると思います。 どこまで0にする必要があるかは、周波数で決めてください。

ざっくりこんな感じ

python

1import numpy as np 2import matplotlib.pyplot as plt 3#matplotlib inline 4from scipy import signal 5 6N = 2048 # サンプル数 7dt = 0.01 # サンプリング周期 8freq = 10 # 周波数 9amp = 1 # 振幅 10t = np.arange(0, N*dt, dt) # 時間軸 11f = amp * np.sin(2*np.pi * freq*t) + 0.5 * freq*t # 実験用の正弦波 12f0 = amp * np.sin(2*np.pi * freq*t) # 目標とする正弦波 13 14F = np.fft.fft(f) # 高速フーリエ変換(FFT) 15plt.subplot( 311) 16plt.plot( f) 17 18plt.subplot( 312) 19plt.plot( np.abs( F)) 20F[ 0: 40] = 0.0 # 周波数低めをカット 21plt.plot( np.abs( F)) 22 23rfft = np.fft.ifft( F) # 逆FFT 24 25plt.subplot( 313) 26plt.plot( np.abs( rfft)) 27 28plt.show()

イメージ説明

投稿2019/08/05 06:22

編集2019/08/05 07:07
glp

総合スコア102

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

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

NCC1701

2019/08/06 11:27

なるほど、超低周波が合成されていると解釈するのは気づきませんでした。勉強になります。
guest

0

実験用と目標とする波形のそれぞれのFFTの結果を比べてみてはどうでしょうか?

イメージ説明

投稿2019/08/04 12:19

TaroToyotomi

総合スコア1430

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

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

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

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

ただいまの回答率
85.50%

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

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

質問する

関連した質問