#!/usr/bin/env python
# coding: utf-8

# In[1]:


import numpy as np
import pandas as pd


# In[2]:


# 図化用のモジュール
import holoviews as hv 
hv.extension('bokeh', logo=False)


# ### 元データ
# 
# とある地形データです。Lが平面距離(m）、Zが標高(m)です。Lは0.5m間隔に正規化されています。

# In[3]:


df = pd.read_csv('test_bedElevation.csv', index_col=0)
df


# #### 縦断図

# In[4]:


L, Z = df.L.values, df.Z.values
g = hv.Curve((df.L,df.Z)).options(xlabel='L[m]', ylabel='Z[m]',width=500)
g


# 複雑なデータですが、地形の表面に出現する細かい波長のデータ、つまり、高周波成分の特徴を分析します。

# ### FFTによるフーリエ変換
# 
# numpyのFFTを使います。データの個数とかも気にせず上手に変換してくれます。細かい中身は見てません。

# In[5]:


dx = 0.5
N = len(Z) 
freq = np.fft.fftfreq(N, d=dx)
F = np.fft.fft(Z)


# #### FFT結果

# In[6]:


g = hv.Curve((freq[1:int(len(F)/2)], np.abs(F/N*2)[1:int(len(F)/2)])).options(xlabel='波長の逆数[1/m]', ylabel='振幅[m]',width=500)
g.redim.range(x=(0,0.2)) 


# ### ハイパスフィルター
# 
# 高周波成分のみを使うため、ハイパスフィルターを使います。上記のFFTの結果の低周波成分を0として、フーリエ逆変換を行います。低周波成分の分類は波長が50m以上（逆数で0.02）としました。

# #### フーリエ逆変換

# In[7]:


F2 = np.copy(F)
ind = (freq < 0.02) & (freq > -0.02) 
F2[ind] = 0.0
Z2 = np.fft.ifft(F2)


# In[8]:


g = hv.Curve((freq[1:int(len(F2)/2)], np.abs(F2/N*2)[1:int(len(F2)/2)])).options(xlabel='波長の逆数[1/m]', ylabel='振幅[m]',width=500)
g.redim.range(x=(0,0.2)) 


# In[9]:


g = hv.Curve((1/freq[1:int(len(F2)/2)], np.abs(F2/N*2)[1:int(len(F2)/2)])).options(xlabel='波長[m]', ylabel='振幅[m]',width=500)
g.redim.range(x=(0,50)) 


# #### ハイパスフィルタ後の縦断図

# In[10]:


g = hv.Curve((L,Z2.real)).options(xlabel='L[m]', ylabel='Z[m]',width=500)
g


# いい感じ高周波成分のみ抽出できました。
# 
# こんな感じでみると、0～900m間は、0.2-0.4m程度の波高の成分、900m以降は、0.1-0.2m程度の成分が主となっています。
# 
# 場所によって卓越する成分が異なっていそうな感じです。

# ### おまけ
# 
# 上記を参照に2区間に分けて、同じ処理をしてみます。

# In[11]:


def highpass(Z, freq):
    F = np.fft.fft(Z)
    F2 = np.copy(F)
    ind = (freq < 0.02) & (freq > -0.02) 
    F2[ind] = 0.0
    Z2 = np.fft.ifft(F2)
    return Z2, F2


# In[12]:


Z1 = Z[:1800]
Z2 = Z[1800:]

dx = 0.5
freq1 = np.fft.fftfreq(len(Z1), d=dx)
Z1h, F1h = highpass(Z1, freq1)
freq2 = np.fft.fftfreq(len(Z2), d=dx)
Z2h, F2h = highpass(Z2, freq2)


# #### FFT結果
# 
# 卓越する波長の差は無さそうです。

# In[13]:


g1 = hv.Curve((1/freq1[1:int(len(Z1)/2)], np.abs(F1h/len(Z1)*2)[1:int(len(Z1)/2)]), label = 's1')
g2 = hv.Curve((1/freq2[1:int(len(Z2)/2)], np.abs(F2h/len(Z2)*2)[1:int(len(Z2)/2)]), label = 's2')
g1*g2.redim.range(x=(0,50)).options(xlabel='波長[m]', ylabel='振幅[m]',width=500)


# #### ハイパスフィルタ後の縦断図

# In[14]:


g1 = hv.Curve((np.array(range(len(Z1h)))*dx,Z1h.real), label = 's1').options(xlabel='L[m]', ylabel='Z[m]',width=500)
g2 = hv.Curve((np.array(range(len(Z2h)))*dx,Z2h.real), label = 's2').options(xlabel='L[m]', ylabel='Z[m]',width=500)
g1 * g2