import numpy as np
import math
from scipy.fftpack import fft,ifft
import matplotlib.pyplot as plt
from matplotlib.pylab import mpl
import os,sys
from scipy.signal import butter, lfilter, freqz
import matplotlib.pyplot as plt


def butter_lowpass(cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='low', analog=False)
    return b, a


def butter_lowpass_filter(data, cutoff, fs, order=5):
    b, a = butter_lowpass(cutoff, fs, order=order)
    y = lfilter(b, a, data)
    return y  # Filter requirements.

def plot_fft(file_name):
    order = 6
    fs = 114*1000 # sample rate, Hz
    cutoff = 20*1000 # desired cutoff frequency of the filter, Hz # Get the filter coefficients so we can check its frequency response.

    dir = '/run/media/root/35b2fb91-52d3-4e4e-a8ec-4d3a9e587ece/fanzaiyun/plot/'
    #f = open(r'/run/media/root/35b2fb91-52d3-4e4e-a8ec-4d3a9e587ece/fanzaiyun/PGC/25KM-2.txt')             # 返回一个文件对象
    f = open(file_name)             # 返回一个文件对象

    postfix = file_name.rfind(".txt")
    prefix = file_name.rfind("/")
    path = dir + file_name[prefix:postfix] 
    if not os.path.isdir(path):
        os.mkdir(path)
    print("current path : ", path)
    line = f.readline()             # 调用文件的 readline()方法
    while line:
        #print line,                 # 在 Python 2中，后面跟 ',' 将忽略换行符
        #print(line, end = '')       # 在 Python 3中使用
        datay = line.split(" ")
        print(len(datay), datay[0])
        line = f.readline()
    f.close()

    data = []
    for m in datay:
        m = m.replace(" ","")
        m = m.replace("\n","")
        if len(m):
            k = int(m)
            data.append(k)

    #for j in data:
    #    print(j)

    print("len=",len(data))
    #mpl.rcParams['font.sans-serif'] = ['SimHei']   #显示中文
    mpl.rcParams['axes.unicode_minus']=False       #显示负号
    
    
    #采样点选择1400个，因为设置的信号频率分量最高为600赫兹，根据采样定理知采样频率要大于信号频率2倍，所以这里设置采样频率为1400赫兹（即一秒内有1400个采样点，一样意思的）
    x=np.linspace(0,30000,len(data))      
    
    #设置需要采样的信号，频率分量有200，400和600
    plt.figure()
    plt.ylim(2000, 4500)
    plt.plot(x,data, ls="-")   
    plt.title('raw wave')
    plt.savefig(path+"/raw.png")
    plt.show()
    
    plt.figure()
    plt.plot(x[00:5200],data[00:5200])   
    plt.title('raw wave(0~5200 header data)')
    plt.savefig(path+"/raw-0-5200.png")
    plt.show()

    data_pad = []
    for m in data:
        data_pad.append(m)
    i = 30000
    while 1:
        data_pad.append(0)
        i =  i + 1
        if i >= 32768:
            break
    i = 0
    for m in data_pad:
        i = i + 1
        #print(i, ":", m)

    data_lowpass = butter_lowpass_filter(data_pad, cutoff, fs, order)
    
    plt.figure()
    plt.plot(x[00:30000],data_lowpass[00:30000])   
    plt.title('lowpass wave(0~30000 header data)')
    plt.savefig(path+"/lowpass-0-30000.png")
    plt.show()

    bw = 4096
    start = 0
    print("len of data_lowpass : " , len(data_lowpass))
    print("fft data")
    fft_y=fft(data_lowpass[start:start+bw-1])
    i = 0   
    print("raw data")
    for m in fft_y:
        i = i + 1
        #print( i, ":", m)

    #complex square
    sq = [ ]
    j = 0
    for m in fft_y:
        r = m.real
        i = m.imag
        k = math.log(r*r+i*i)*10
        #if j < 100:
        #    sq.append(0)
        #else:
        sq.append(k)
        j = j + 1

    i = 0
    print("sq data :")
    for m in sq:
        i = i + 1
        #print(i, ":", m)

    N=32768
    x = np.arange(N)           # 频率个数
    
    print("len of x" , len(x))
    #abs_y=np.abs(fft_y)                # 取复数的绝对值，即复数的模(双边频谱)
    #angle_y=np.angle(fft_y)              #取复数的角度
    
    plt.figure()
    plt.plot(x[0:500],sq[00:500])   
    plt.title('fft result fft()')
    plt.savefig(path+"/fft.png")
    plt.show()
    
    #plt.figure()
    #plt.plot(x,angle_y)   
    #plt.title('angle fft result')
    #plt.show()

    #N=32768
    #x = np.arange(N)             # 频率个数
    #half_x = x[range(int(N/2))]  #取一半区间
    
    #abs_y=np.abs(fft_y)                # 取复数的绝对值，即复数的模(双边频谱)
    #angle_y=np.angle(fft_y)            #取复数的角度
    #normalization_y=abs_y/N            #归一化处理（双边频谱）                              
    #normalization_half_y = normalization_y[range(int(N/2))]      #由于对称性，只取一半区间（单边频谱）
    #y = data 
    #plt.subplot(231)
    #plt.plot(x,y)   
    #plt.title('raw wave ')
    
    #plt.subplot(232)
    #plt.plot(x,fft_y,'black')
    #plt.title('shuangbian zhenfu',fontsize=9,color='black') 
    
    #plt.subplot(233)
    #plt.plot(x,abs_y,'r')
    #plt.title('zhenfu(wei guiyihua)',fontsize=9,color='red') 
    
    #plt.subplot(234)
    #plt.plot(x,angle_y,'violet')
    #plt.title('xiangweipu (wei guiyihua)',fontsize=9,color='violet')
    
    #plt.subplot(235)
    #plt.plot(x,normalization_y,'g')
    #plt.title('zhenfu (guiyihua)',fontsize=9,color='green')
    
    #plt.subplot(236)
    #plt.plot(half_x,normalization_half_y,'blue')
    #plt.title('danbian zhenfu(guiyihua',fontsize=9,color='blue')
    
    #plt.show()

#plot_fft('/run/media/root/35b2fb91-52d3-4e4e-a8ec-4d3a9e587ece/fanzaiyun/PGC/25KM-4.txt')

plot_fft('/run/media/root/35b2fb91-52d3-4e4e-a8ec-4d3a9e587ece/fanzaiyun/PGC/50KM-4.txt')
#plot_fft('/run/media/root/35b2fb91-52d3-4e4e-a8ec-4d3a9e587ece/fanzaiyun/PGC/50KM-5.txt')
#for start in range(1,5):
#    plot_fft('/run/media/root/35b2fb91-52d3-4e4e-a8ec-4d3a9e587ece/fanzaiyun/PGC/25KM-'+str(start)+'.txt')
##    plot_fft('/run/media/root/35b2fb91-52d3-4e4e-a8ec-4d3a9e587ece/fanzaiyun/PGC/50KM-'+str(start)+'.txt')