import wave
import numpy as np
import math
import matplotlib.pyplot as plt
from scipy.fftpack import dct

def read(data_path):
 '''讀取語(yǔ)音信號(hào)
 '''
 wavepath = data_path
 f = wave.open(wavepath,'rb')
 params = f.getparams()
 nchannels,sampwidth,framerate,nframes = params[:4] #聲道數(shù)、量化位數(shù)、采樣頻率、采樣點(diǎn)數(shù)
 str_data = f.readframes(nframes) #讀取音頻，字符串格式
 f.close()
 wavedata = np.fromstring(str_data,dtype = np.short) #將字符串轉(zhuǎn)化為浮點(diǎn)型數(shù)據(jù)
 wavedata = wavedata * 1.0 / (max(abs(wavedata))) #wave幅值歸一化
 return wavedata,nframes,framerate

def enframe(data,win,inc):
 '''對(duì)語(yǔ)音數(shù)據(jù)進(jìn)行分幀處理
 input:data(一維array):語(yǔ)音信號(hào)
   wlen(int):滑動(dòng)窗長(zhǎng)
   inc(int):窗口每次移動(dòng)的長(zhǎng)度
 output:f(二維array)每次滑動(dòng)窗內(nèi)的數(shù)據(jù)組成的二維array
 '''
 nx = len(data) #語(yǔ)音信號(hào)的長(zhǎng)度
 try:
  nwin = len(win)
 except Exception as err:
  nwin = 1 
 if nwin == 1:
  wlen = win
 else:
  wlen = nwin
 nf = int(np.fix((nx - wlen) / inc) + 1) #窗口移動(dòng)的次數(shù)
 f = np.zeros((nf,wlen)) #初始化二維數(shù)組
 indf = [inc * j for j in range(nf)]
 indf = (np.mat(indf)).T
 inds = np.mat(range(wlen))
 indf_tile = np.tile(indf,wlen)
 inds_tile = np.tile(inds,(nf,1))
 mix_tile = indf_tile + inds_tile
 f = np.zeros((nf,wlen))
 for i in range(nf):
  for j in range(wlen):
   f[i,j] = data[mix_tile[i,j]]
 return f

def point_check(wavedata,win,inc):
 '''語(yǔ)音信號(hào)端點(diǎn)檢測(cè)
 input:wavedata(一維array)：原始語(yǔ)音信號(hào)
 output:StartPoint(int):起始端點(diǎn)
   EndPoint(int):終止端點(diǎn)
 '''
 #1.計(jì)算短時(shí)過零率
 FrameTemp1 = enframe(wavedata[0:-1],win,inc)
 FrameTemp2 = enframe(wavedata[1:],win,inc)
 signs = np.sign(np.multiply(FrameTemp1,FrameTemp2)) # 計(jì)算每一位與其相鄰的數(shù)據(jù)是否異號(hào)，異號(hào)則過零
 signs = list(map(lambda x:[[i,0] [i>0] for i in x],signs))
 signs = list(map(lambda x:[[i,1] [i<0] for i in x], signs))
 diffs = np.sign(abs(FrameTemp1 - FrameTemp2)-0.01)
 diffs = list(map(lambda x:[[i,0] [i<0] for i in x], diffs))
 zcr = list((np.multiply(signs, diffs)).sum(axis = 1))
 #2.計(jì)算短時(shí)能量
 amp = list((abs(enframe(wavedata,win,inc))).sum(axis = 1))
# # 設(shè)置門限
# print('設(shè)置門限')
 ZcrLow = max([round(np.mean(zcr)*0.1),3])#過零率低門限
 ZcrHigh = max([round(max(zcr)*0.1),5])#過零率高門限
 AmpLow = min([min(amp)*10,np.mean(amp)*0.2,max(amp)*0.1])#能量低門限
 AmpHigh = max([min(amp)*10,np.mean(amp)*0.2,max(amp)*0.1])#能量高門限
 # 端點(diǎn)檢測(cè)
 MaxSilence = 8 #最長(zhǎng)語(yǔ)音間隙時(shí)間
 MinAudio = 16 #最短語(yǔ)音時(shí)間
 Status = 0 #狀態(tài)0:靜音段,1:過渡段,2:語(yǔ)音段,3:結(jié)束段
 HoldTime = 0 #語(yǔ)音持續(xù)時(shí)間
 SilenceTime = 0 #語(yǔ)音間隙時(shí)間
 print('開始端點(diǎn)檢測(cè)')
 StartPoint = 0
 for n in range(len(zcr)):
  if Status ==0 or Status == 1:
   if amp[n] > AmpHigh or zcr[n] > ZcrHigh:
    StartPoint = n - HoldTime
    Status = 2
    HoldTime = HoldTime + 1
    SilenceTime = 0
   elif amp[n] > AmpLow or zcr[n] > ZcrLow:
    Status = 1
    HoldTime = HoldTime + 1
   else:
    Status = 0
    HoldTime = 0
  elif Status == 2:
   if amp[n] > AmpLow or zcr[n] > ZcrLow:
    HoldTime = HoldTime + 1
   else:
    SilenceTime = SilenceTime + 1
    if SilenceTime < MaxSilence:
     HoldTime = HoldTime + 1
    elif (HoldTime - SilenceTime) < MinAudio:
     Status = 0
     HoldTime = 0
     SilenceTime = 0
    else:
     Status = 3
  elif Status == 3:
   break
  if Status == 3:
   break
 HoldTime = HoldTime - SilenceTime
 EndPoint = StartPoint + HoldTime
 return FrameTemp1[StartPoint:EndPoint]


def mfcc(FrameK,framerate,win):
 '''提取mfcc參數(shù) 
 input:FrameK(二維array):二維分幀語(yǔ)音信號(hào)
   framerate:語(yǔ)音采樣頻率
   win:分幀窗長(zhǎng)（FFT點(diǎn)數(shù)）
 output:
 '''
 #mel濾波器
 mel_bank,w2 = mel_filter(24,win,framerate,0,0.5)
 FrameK = FrameK.T
 #計(jì)算功率譜
 S = abs(np.fft.fft(FrameK,axis = 0)) ** 2
 #將功率譜通過濾波器
 P = np.dot(mel_bank,S[0:w2,:])
 #取對(duì)數(shù)
 logP = np.log(P)
 #計(jì)算DCT系數(shù)
# rDCT = 12
# cDCT = 24
# dctcoef = []
# for i in range(1,rDCT+1):
#  tmp = [np.cos((2*j+1)*i*math.pi*1.0/(2.0*cDCT)) for j in range(cDCT)]
#  dctcoef.append(tmp)
# #取對(duì)數(shù)后做余弦變換 
# D = np.dot(dctcoef,logP)
 num_ceps = 12
 D = dct(logP,type = 2,axis = 0,norm = 'ortho')[1:(num_ceps+1),:]
 return S,mel_bank,P,logP,D
 


def mel_filter(M,N,fs,l,h):
 '''mel濾波器
 input:M(int)：濾波器個(gè)數(shù)
   N(int)：FFT點(diǎn)數(shù)
   fs(int)：采樣頻率
   l(float)：低頻系數(shù)
   h(float)：高頻系數(shù)
 output:melbank(二維array):mel濾波器
 '''
 fl = fs * l #濾波器范圍的最低頻率
 fh = fs * h #濾波器范圍的最高頻率
 bl = 1125 * np.log(1 + fl / 700) #將頻率轉(zhuǎn)換為mel頻率
 bh = 1125 * np.log(1 + fh /700) 
 B = bh - bl #頻帶寬度
 y = np.linspace(0,B,M+2) #將mel刻度等間距
 print('mel間隔',y)
 Fb = 700 * (np.exp(y / 1125) - 1) #將mel變?yōu)镠Z
 print(Fb)
 w2 = int(N / 2 + 1)
 df = fs / N
 freq = [] #采樣頻率值
 for n in range(0,w2):
  freqs = int(n * df)
  freq.append(freqs)
 melbank = np.zeros((M,w2))
 print(freq)
 
 for k in range(1,M+1):
  f1 = Fb[k - 1]
  f2 = Fb[k + 1]
  f0 = Fb[k]
  n1 = np.floor(f1/df)
  n2 = np.floor(f2/df)
  n0 = np.floor(f0/df)
  for i in range(1,w2):
   if i >= n1 and i <= n0:
    melbank[k-1,i] = (i-n1)/(n0-n1)
   if i >= n0 and i <= n2:
    melbank[k-1,i] = (n2-i)/(n2-n0)
  plt.plot(freq,melbank[k-1,:])
 plt.show()
 return melbank,w2

if __name__ == '__main__':
 data_path = 'audio_data.wav'
 win = 256
 inc = 80
 wavedata,nframes,framerate = read(data_path)
 FrameK = point_check(wavedata,win,inc)
 S,mel_bank,P,logP,D = mfcc(FrameK,framerate,win)