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

Question

前提・実現したいこと

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

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

該当のソースコード

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from scipy import signal

N = 64 # サンプル数
dt = 0.01 # サンプリング周期
freq = 10 # 周波数
amp = 1 # 振幅
t = np.arange(0, N*dt, dt) # 時間軸
f = amp * np.sin(2*np.pi * freq*t) + 0.5 * freq*t # 実験用の正弦波
f0 = amp * np.sin(2*np.pi * freq*t) # 目標とする正弦波

F = np.fft.fft(f) # 高速フーリエ変換(FFT)
F_abs = np.abs(F)
F_abs_amp = F_abs / N * 2
F_abs_amp[0] = F_abs_amp[0] / 2

angle =  F_abs[0]/N/N*20 # 傾きを取得しているつもり

fq = np.linspace(0, 1.0/dt, N)
fig = plt.figure(figsize=(12, 4))
ax2 = fig.add_subplot(121) ## 信号のグラフ
plt.xlabel('time(sec)')
plt.ylabel('amp')
plt.plot(t, f)
plt.plot(t, f0)
ax2 = fig.add_subplot(122) ## FFTのグラフ（周波数軸）
plt.xlabel('freqency(Hz)')
plt.ylabel('amp')
plt.plot(fq[:int(N/2)+1], F_abs_amp[:int(N/2)+1])
plt.show()

試したこと

F_ifft = np.fft.ifft(F) # IFFT
F_ifft_real = F_ifft.real * 2 # 振幅を元スケールに戻している

F_minus = F_ifft_real - angle*t # ここで傾きを除去しているつもり

plt.plot(t, F_ifft_real)
plt.plot(t, F_minus, c="r")
plt.show()

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

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

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

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

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

python 3（Colaboratory）
参考にしたサイト

Answer 1

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

ざっくりこんな感じ

import numpy as np
import matplotlib.pyplot as plt
#matplotlib inline
from scipy import signal

N = 2048 # サンプル数
dt = 0.01 # サンプリング周期
freq = 10 # 周波数
amp = 1 # 振幅
t = np.arange(0, N*dt, dt) # 時間軸
f = amp * np.sin(2*np.pi * freq*t) + 0.5 * freq*t # 実験用の正弦波
f0 = amp * np.sin(2*np.pi * freq*t) # 目標とする正弦波

F = np.fft.fft(f) # 高速フーリエ変換(FFT)
plt.subplot( 311)
plt.plot( f)

plt.subplot( 312)
plt.plot( np.abs( F))
F[ 0: 40] = 0.0			# 周波数低めをカット
plt.plot( np.abs( F))

rfft = np.fft.ifft( F)	# 逆FFT

plt.subplot( 313)
plt.plot( np.abs( rfft))

plt.show()

Comments

[NCC1701]

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

Answer 2

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