import numpy as np
import matplotlib.pyplot as plt
import wave
import math

pltimg=True

# 出力する声の高さ [Hz]
# 自分の声の特性にあわせて適宜調整する(100〜400程度)
VOICEFREQ = ?

# 自己相関関数 R[n] を求める関数
# 自己相関関数の演習で作成した ACR 関数の中身をコピーする
def ACR( n, f, N ):
?

# 自己回帰モデル
# 演習 1-1 からコピー
def AR(a,e):
?

# レビンソン・ダービン(Levinson-Durbin)・アルゴリズム
# 演習 1-4 からコピー
def Levinson(x,L):
?


# Waveファイル読み込み関数
def read_wavedata(name):
    print('read: '+name)
    with wave.open(name,'r') as wavein:
        print(wavein.getparams())
        rate = wavein.getframerate()
        nframes = wavein.getnframes()
        data = wavein.readframes(wavein.getnframes())
    return np.frombuffer(data, dtype=np.int16),rate,nframes

# Waveファイル書き込み関数
def write_wavedata(name,rate,data):
    print('write: '+name)
    with wave.open( name, 'wb') as waveout:
        rawdata=np.array(data,dtype='int16').tobytes()
        waveout.setparams((1, 2, rate, len(rawdata), 'NONE', 'NONE'))
        waveout.writeframes(rawdata)

# 音量ノーマライズ関数
def normalize_wave(rms,data):
    r = 0;
    for i in range(data.size):
        r += data[i]**2
    r = math.sqrt(r/data.size)
    gain = rms/r
    print(f'RMS = {r:.02f}, gain = {gain:.02f}')
    for i in range(data.size):
        data[i] *= gain

# スペクトログラム描画関数
def plot_spectrogram(data,file):
    cmap = plt.get_cmap('magma')
    plt.specgram(data/32768,Fs=fs,NFFT=512,scale='dB',mode='magnitude',vmax=0,vmin=-120,cmap=cmap)
    plt.colorbar()
    plt.xlabel('sec')
    plt.ylabel('Hz')
    if pltimg == True:
        plt.show()
    else:
        plt.savefig(file)

L = 120 # 次数
VOICELNG = 1.5 # 声ひとつ分の長さ(秒)
RMS = 3500 # 出力音声の音量(RMS)
STARTPOS = 1.0 # wave ファイルの読み込み開始位置 (秒)
DURATION = 0.5 # wave ファイルの読み込みサイズ (秒)
AMP = 500.0 # 三角波の大きさ
SIGMA = AMP/4.0  # ノイズの標準偏差

files = ['a','i','u','e','o']
fs = 0
hat_as = []

# 各音声の LPC 係数を推定
for file in files:
    voice,fs,N = read_wavedata(file+'.wav')
    duration = N/fs
    print(f'duration={duration:.02f} sec')
    if (duration < STARTPOS + DURATION ):
        raise Exception('録音時間が短いです')
    voice = voice[int(STARTPOS*fs):int((STARTPOS+DURATION)*fs)].astype(np.float64)
    hat_a,E = Levinson(voice,L)
    hat_as.append(hat_a)
    
# 入力パルス音 e[i] (周期 T [点] の三角波) の作成
e = np.zeros(int(VOICELNG*fs))
T = int( float(fs)/VOICEFREQ )  # 周期
for i in range(e.size):
    e[i] = -AMP/2 + AMP/T*(i%T) + np.random.normal(0,SIGMA)

# e[i] を自己回帰モデルに入力してxを出力
x = []
for hat_a in hat_as:
    tmpx = AR(hat_a,e)
    normalize_wave(RMS,tmpx)
    x = np.r_[x,tmpx]

normalize_wave(RMS,e)

write_wavedata('AR_1_6_e.wav',fs,e)
write_wavedata('AR_1_6_x.wav',fs,x)

print('')
print('入力')
plot_spectrogram(e,'AR_1_6_e.png')

print('')
print('出力')
plot_spectrogram(x,'AR_1_6_x.png')