#範例12-2：D.T.F.T：離散資料，非週期數據，的傅立葉轉換法FFT(Fast Fourier Transform)
#sympy要處理：離散性，非週期數據，乃是用快速傅立葉轉換法FFT(Fast Fourier Transform)
#語法：fft(seq, dps=None)
#注意：sympy要給定定的離散數據量，要是2*n倍，若不夠，系統會補0
#缺點：sympy的fft()函數，不好用。若給定太多的x數據，計算會變得很複雜，最後程式會當機
#使用sympy的FFT函數來做快速傅立葉轉換，有很多個嚴重缺點：
#（1）.讀入的離散資料，不能太多，讀入20個還可以計算，讀入200個數據就會無法計算
#（2）.fft轉換後的數據資料都是未簡化的數，若要將之轉成估值（numerical evaluation)會很麻煩
#（3）fft後續的陣列，sympy不提供對整個array陣列list的直接萃取數值（例如：只能對單一值做估值：.evalf()或是N(..)
#建議：使用：numpy或scipy，它們都有直接對陣列list操作的各種函數，而且轉換後的fft陣列，也是自動進行估值計算了

#Inverse FT： 非週期性的連續波之傅立葉轉換語法： InverseFourierTransform(ft, k, x).doit()
from sympy import *
import numpy as np
import matplotlib.pyplot as plt
#from sympy.plotting import plot, PlotGrid
x = symbols('x')
k = symbols('k')
f = Function('f')(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 = fft(sig)
#print(sig_fft)
#這些複數都是未被簡化或估值的數據

#算出fft的水平頻率軸
#因為sympy這裡提供的函數不夠，所以使用scipy的函數
from scipy import fftpack
sample_freq = fftpack.fftfreq(len(sig), d = time_step)
print('sample_freq=',sample_freq)

#算出fft的垂直頻率軸振幅
#amplitude = np.abs(sig_fft)
#發現：np.abs()，有些複數沒有計算，所以必須要手動計算
#print('abs=', amplitude[0],sample_freq[0],len(amplitude),len(sample_freq))
#amplitude = np.sqrt(np.real(sig_fft)**2 + np.imag(sig_fft)**2)
amplitude = []
for item in sig_fft:
    amplitude.append(N(Abs(N(item))))
#sympy數值估值函數：N()，或是.evalf()
#sympy計算複數絕對值 = sqrt(a**2+b**2) = Abs() = 振幅  
print('amplitude=',amplitude)
plt.subplot(3,1,2)
plt.plot(sample_freq, amplitude, label='FFTFFT(frequency spectrum)')
plt.legend()

#(3).inver FFT：：把頻域，反轉成時域
sig_inverse_fft = ifft(sig_fft)
print('inverse fft=', sig_inverse_fft)
#這些複數都是未被簡化或估值的數據
plt.subplot(3,1,3)
plt.plot(time_vec, sig_inverse_fft,'r-', label='inverse fft data')
plt.legend()