#範例12-5：從2個sin合成波訊號中濾除雜訊
#訊號的頻率=1 kHz，訊號間的總量測時間是=1.5 seconds.
#訊號的頻率(Sampling frequency)：Fs = 1000 Hz
#訊號的間隔時間= T = 1/Fs = 0.001
#訊號數目(Length of signal)= L = 1500
#時間矩陣(time_vector)= t = (0:L-1)*T
#模擬訊號：由兩個sin(x)波組合而成（50 Hz sinusoid of amplitude 0.7）+（120 Hz sinusoid of amplitude 1）
#合成波 S = 0.7*sin(2πt*50) + 1*sin(2πt*120)
#加上干擾雜訊： 2*random.ranint(時間矩陣點數的數目) =  2*random.ranint(time_vector.size)
#轉換成頻域，會發現不止50Hz, 120Hz，還有很多頻率的雜訊，要將之濾除

import numpy as np
from scipy import fftpack
import matplotlib.pyplot as plt
#(1).設定時間軸矩陣向量，頻率，訊號數目
Fs = 1000
T =  1/Fs
time_vec = np.arange(0, 1.5, T)
L = time_vec.size
print('L=',L, ', T=', T, ', Fs=', Fs)
print('time_vec=', time_vec)

#(2).建立2個sin的組合波（50 Hz sinusoid of amplitude 0.7）+（120 Hz sinusoid of amplitude 1）
#sig2 = 0.7*np.sin(2*np.pi*time_vec*50) + 1*np.sin(2*np.pi*time_vec*120)
sig2 = 0.7*np.sin(2*np.pi*time_vec*50) + 1*np.sin(2*np.pi*time_vec*120)
#np.random.randn(n)：建立一個一維n元素亂數矩陣, 內容是自然分佈，平均=0， 變化量=1(distribution of mean 0 and variance 1)
#np.random.randn(m, n)：建立一個一維mxn亂數矩陣
sig_nose = 2*np.random.randn(L)
sig =  sig2 + sig_nose
#sig =  sig2
print('sig=', sig)

#(3).畫出時域訊號圖(show the first 50 signal)
plt.subplot(3,1,1)
plt.plot(1000*time_vec[0:50], sig[0:50], label='raw signal')
plt.legend()
#(4).用FFT=轉成頻域(用scipy的fft函數，效率較好)
sig_fft = fftpack.fft(sig)

#(5).計算頻域的水平軸頻率矩陣
#scipy語法：fftpack.fftfreq(sig_fft.size, d=訊號間隔時間)
sample_freq = fftpack.fftfreq(sig_fft.size, d=T)

#(6).計算頻域的振幅矩陣
#amp2 = np.array()
amplitude = np.abs(sig_fft)/L
#只看頻率為正值部分，振幅正負部位合在一起，2倍
amplitude2 = 2*amplitude

#(7).畫圖：頻域頻譜圖
plt.subplot(3,1,2)
#頻率報復只顯示一半:0～L/2，這是正值部分
plt.plot(sample_freq[0:int(L/2)], amplitude2[0:int(L/2)], color='red',label='fft(frequency)')
plt.legend()

#(8).印出振幅最大1st的頻率
index_peak_1stmax = amplitude2[:].argmax()
#print('index_peak_1stmax=', index_peak_1stmax)
freq_peak_1smax = sample_freq[index_peak_1stmax]
print('freq_peak_1smax=', freq_peak_1smax)
amplitude_peak_1smax = amplitude2[index_peak_1stmax]
print('amplitude_peak_1smax=', amplitude_peak_1smax)

#(9).濾波：把最大振幅波頻率以外的全部的波，其振幅都設定為0（沒有雜訊）
#最大漲幅波頻率, A < 0.5的波，其sig_fft，都設定為0
high_sig_ft = sig_fft.copy()
high_sig_ft[np.abs(amplitude2) < 0.4] = 0
#print(high_freq_fft)

#(10).反向Inverse推算原函數（由已經濾波後的fft=high_freq_fft來反推）
sig_inverse_fft = fftpack.ifft(high_sig_ft)

#(11).在下方畫出sub plot的方法
plt.subplot(3,1,3)
plt.plot(time_vec[0:50], sig_inverse_fft[0:50],'r-', label='inverse fft data')
plt.legend()
plt.show()