#範例12-3：D.T.F.T：使用spicy轉換離散資料，非週期數據的傅立葉轉換法FFT(Fast Fourier Transform)
#reference:https://www.twblogs.net/a/5b8d116d2b717718833a60d8
#numpy與spicy都可以出來傅立葉轉換，而且提供很多支援的函數
#建議：因為spicy執行的效率較快，所以建議用spicy來處理快速傅立葉轉換
#(1).spicy的FFT語法：scipy.fftpack.fft(sig)
#或是：scipy.fftpack.fft(sig, n = 10)
#第一個參數sig：輸入的訊號序列或陣列list
#第二個參數n：輸入的訊號點數目 = len(sig) = sig.size
#注意：n最好是2*n倍，若不夠，系統會補0

#(2).把轉換後fft，由小到大排列，語法：scipy.fftpack.fftshift(sig_fft)
#注意：若用plot()畫frequency spectrum，則會自動大小排序，不需奧在fftshift()

#(3).計算頻譜圖的頻率陣列，語法：fftpack.fftfreq(n, d = time_step)
#其中 n = 訊號點數 = len(sig) = sig.size
#time step = d = Δt
#採樣點頻率 = fs = 1/d
#採樣點的頻率分辨率 = k*fs/N = k*1/(d*N)
#範例：sample_freq = fftpack.fftfreq(len(sig), d = time_step)

#(4).把頻率陣列重新排序(讓0頻率在中央)，語法：sample_freq_shift = fftpack.fftshift(sample_freq)
#不但頻率可以把0放在中央，fft後的複數函數也會依據頻率的排序而調整順序（會讓0頻率在中央）

#(5).計算頻譜圖的振幅陣列：因為頻率的順序已經改變，所以amplitude也要改版
#amplitude 公式：amplitude_shift = np.sqrt(np.real(sig_fft)**2 + np.imag(sig_fft)**2)

#(6).手動計算頻譜圖的振幅陣列(若是沒有shift頻率，那麼就可以自行手動計算amplitude)：
#方法二：Amplitude = np.abs(sig_fft)
#方法二：
# amplitude = []
# for item in sig_fft:
#     amplitude.append(N(Abs(N(item))))
#sympy數值估值函數：N()，或是.evalf()
#sympy計算複數絕對值 = sqrt(a**2+b**2) = Abs() = 振幅  

#(7).Inverse FT 語法： sig_inverse_fft = ifft(sig_fft)

from sympy import *
import numpy as np
from scipy import fftpack
import matplotlib.pyplot as plt
#x = symbols('x')

#(1)原始量測數據
time_step = 0.1
time_vec = np.arange(0, 0.8, time_step)
print('(1).time=', time_vec)

sig = [0, 0.707, 1, 0.707, 0, -0.707, -1, -0.707]
#sig = [-0.2, 0.757, 0.95, 0.687, 0.02, -0.727, -0.95, -0.727]
print('(2).signal data =', sig)
plt.subplot(3,1,1)
plt.plot(time_vec, sig, label='raw data')
plt.legend()

#(2). FFTFFT：把時域，轉成頻域
#使用sympy的FFT快速傅立葉轉換函數fft，轉換後是個複數
sig_fft = fftpack.fft(sig)
print('sig_fft=', sig_fft)
#依據頻率的排序而調整順序（會讓0頻率在中央）
sig_shift =fftpack.fftshift(sig_fft)
print('sig_shift=', sig_shift)

#算出fft的水平頻率軸
#因為sympy這裡提供的函數不夠，所以使用scipy的函數
sample_freq = fftpack.fftfreq(len(sig), d = time_step)
print('sample_freq=',sample_freq)
#把頻率陣列重新排序(讓0頻率在中央)，語法：fftpack.fftshift(sample_freq)
sample_freq_shift =fftpack.fftshift(sample_freq)
print('sample_freq_shift=',sample_freq_shift)

#算出fft的垂直頻率軸振幅
#amplitude = np.abs(sig_fft)
#發現：np.abs()，有些複數沒有計算，所以必須要手動計算
#print('abs=', amplitude[0],sample_freq[0],len(amplitude),len(sample_freq))

#另外，因為你把frequency的順序shift了，所以amplitude也必須依照sample_freq_shift的順序，手動計算
#amplitude 公式：amplitude_shift = np.sqrt(np.real(sig_fft)**2 + np.imag(sig_fft)**2)
amplitude_shift = np.sqrt(np.real(sig_shift)**2 + np.imag(sig_shift)**2)
#amplitude_shift = []
#for item in sig_shift:
#    amplitude_shift.append(N(Abs(N(item))))
#sympy數值估值函數：N()，或是.evalf()
#sympy計算複數絕對值 = sqrt(a**2+b**2) = Abs() = 振幅  
print('amplitude=',amplitude_shift)
plt.subplot(3,1,2)
plt.plot(sample_freq_shift, amplitude_shift, label='FFTFFT(frequency spectrum)')
plt.legend()

#(3).inver FFT：：把頻域，反轉成時域
#注意要invers，不可以針對sig_shift，因為它的順序以及被改變了
sig_inverse_fft = fftpack.ifft(sig_fft)
sig_inverse_fft_shift =fftpack.fftshift(sig_inverse_fft)
print('inverse fft=', sig_inverse_fft_shift)
#這些複數都是未被簡化或估值的數據
plt.subplot(3,1,3)
plt.plot(time_vec, sig_inverse_fft,'r-', label='inverse fft data')
plt.legend()