# -*- coding: utf-8 -*-
"""
Created on Tue Dec 20 11:36:31 2022

@author: camdr
"""

import numpy as np
#import pylab      
from scipy.io.wavfile import write
import random
from scipy import signal
import matplotlib.pyplot as plt
from scipy.io import wavfile
from scipy.interpolate import interp1d
import librosa
global alp
global alp2
global feedback1 
global estimulo1
global destimulodt1
global b
global b2
    

# --------------------
# Function definitions
# --------------------
def ecuaciones(v, dv):
    x,y,i1,i2,i3, z, w = v
    dv[0]=y
    dv[1]=gamma*gamma*alp+b*gamma*gamma*x-gamma*gamma*x*x*x-gamma*x*x*y+gamma*gamma*x*x-gamma*x*y+0.*(w-y)
    dv[2]= i2
    dv[3] = - i1/(Lg*Ch) - Rh*i2*(1/Lb + 1/Lg) + i3*( 1/(Lg*Ch) - (Rb*Rh)/(Lb*Lg) ) + destimulodt/Lg + estimulo*Rh / (Lg*Lb)
    dv[4] = - Lg*i2/Lb - Rb*i3 / Lb + estimulo/Lb
    dv[5]=w
    dv[6]=gamma*gamma*alp2+b2*gamma*gamma*z-gamma*gamma*z*z*z-gamma*z*z*w+gamma*gamma*z*z-gamma*z*w+0.*(y-w)
    return dv


def expo(ti,tf,wi,wf,factor,fch,frequencias,amplitudes,ch):
    i=int(ti/dt)
    j=int(tf/dt)
    for k in range((j-i)):
        t=ti+k*dt
        frequencias[i+k]=wf+(wi-wf)*np.exp(-3*(t-ti)/((tf-ti)))
        alpha[i+k]=-0.15
        amplitudes[i+k]=1.*factor 
        ch[i+k]=fch*ch[i+k]
    return frequencias,amplitudes
def expo2(ti,tf,wi,wf,factor,fch,frequencias2,amplitudes2,ch):
    i=int(ti/dt)
    j=int(tf/dt)
    for k in range((j-i)):
        t=ti+k*dt
        frequencias2[i+k]=wf+(wi-wf)*np.exp(-3*(t-ti)/((tf-ti)))
        alpha[i+k]=-0.15
        amplitudes2[i+k]=1.*factor 
        ch[i+k]=fch*ch[i+k]
    return frequencias2,amplitudes2



def rectas(ti,tf,wi,wf,factor,fch,frequencias,amplitudes,ch):
    i=int(ti/dt)
    j=int(tf/dt)
    for k in range((j-i)):
        t=ti+k*dt
        frequencias[i+k]=wi+(wf-wi)*(t-ti)/(tf-ti)
        alpha[i+k]=-0.15
        amplitudes[i+k]= 1.*factor 
        ch[i+k]=fch*ch[i+k]
    return frequencias,amplitudes
def rectas2(ti,tf,wi,wf,factor,fch,frequencias2,amplitudes2,ch):
    i=int(ti/dt)
    j=int(tf/dt)
    for k in range((j-i)):
        t=ti+k*dt
        frequencias2[i+k]=wi+(wf-wi)*(t-ti)/(tf-ti)
        alpha[i+k]=-0.15
        amplitudes2[i+k]= 1.*factor 
        ch[i+k]=fch*ch[i+k]
    return frequencias2,amplitudes2    

def senito(ti,tf,media,amplitud,alphai,alphaf,factor,fch,frequencias,amplitudes,ch):
    i=int(ti/dt)
    j=int(tf/dt)
    for k in range((j-i)):
        t=ti+k*dt
        frequencias[i+k]=media+amplitud*np.sin(alphai+(alphaf-alphai)*(t-ti)/(tf-ti))
        alpha[i+k]=-0.15
        amplitudes[i+k]=1.*factor 
        ch[i+k]=fch*ch[i+k]
    return frequencias,amplitudes
def senito2(ti,tf,media,amplitud,alphai,alphaf,factor,fch,frequencias2,amplitudes2,ch):
    i=int(ti/dt)
    j=int(tf/dt)
    for k in range((j-i)):
        t=ti+k*dt
        frequencias2[i+k]=media+amplitud*np.sin(alphai+(alphaf-alphai)*(t-ti)/(tf-ti))
        alpha[i+k]=-0.15
        amplitudes2[i+k]=1.*factor 
        ch[i+k]=fch*ch[i+k]
    return frequencias2,amplitudes2 

for lazo in range(1):
    gamma=40000
    Ch = (.70 * 1e-8) /350 # Capacitancia Helmholtz
    Lg = 60 # Inductancia glotis Helmholtz
    Rh = 33000 # Resistencia disipación Helmholtz
    Lb = 1e4
    Rb = (2.5)*1e7
    beta, dt, t0, L= -0.0292, 1/882000.0, 0, 0.013
    t=0
    fsamp=1/dt
    samplerate, data_real = wavfile.read('E:/CAMI/Python 2023-2025/Musico/3555 (23-24)/3_mono.wav')
    tiempo_total=20*len(data_real)*(1/882000.)
    frequencias=np.zeros(int(tiempo_total/(dt)))
    tiempos=np.zeros(int(tiempo_total/(dt)))
    alpha=np.zeros(int(tiempo_total/(dt)))
    amplitudes=np.zeros(int(tiempo_total/(dt)))
    amplitudes2=np.zeros(int(tiempo_total/(dt)))
    frequencias2=np.zeros(int(tiempo_total/(dt)))
    beta=np.zeros(int(tiempo_total/(dt)))  
    beta2=np.zeros(int(tiempo_total/(dt)))
    ch=np.zeros(int(tiempo_total/(dt)))
    ch=ch+5800000000.0
    fch=1.0
    cont=0
    tau=0.05 # el tiempo que se ajusta
    env=np.zeros(len(data_real))
    eve=[]
    a=max(data_real)
    for i in range(2,len(data_real)):
        env[i]=env[i-1]+(-tau*env[i-1]+np.abs(data_real[i]/a))
        eve.append(env[i])


    tf=(1.0/samplerate)*len(data_real)
    index_data_max=len(data_real)

    tiempitos=np.linspace(0, tf, num=index_data_max, endpoint=True)
    tiempos=t=np.linspace(0, tf, num=20*index_data_max, endpoint=True)
    envelope=interp1d(tiempitos, env,kind='linear')

 
    bes,was=np.loadtxt('C:/Users/camdr/OneDrive/Desktop/b_w/b_w_40000.txt',unpack=True)
    bes2,was2=np.loadtxt('C:/Users/camdr/OneDrive/Desktop/b_w/b_w_40000.txt',unpack=True)
    
    for i in range(int(tiempo_total/(dt))):
        tiempos[i]=i*dt
        beta[i]=0.15
        beta2[i]=0.15
        amplitudes[i]=0.
        alpha[i]=0.15

    v=np.zeros(7)
    v[0], v[1], v[2], v[3], v[4], v[5], v[6] =0.01,0.001,0.001, 0.0001, 0.0001, 0.00012, 0.0012

# 3555:
   
    rectas(0.04436,0.05354,6292,6930,1,7,frequencias,amplitudes,ch) #3555:3
    senito(0.05354,0.05919,6930,300,1*np.pi,2*np.pi-0.1*np.pi,1,7,frequencias,amplitudes,ch) #3555:3
    senito(0.05919,0.06507,6445,350,1*np.pi/2+0.1*np.pi,1*np.pi-0.0*np.pi,1,7,frequencias,amplitudes,ch) #3555:3
    rectas(0.06507,0.07289,6445,5820,1,7,frequencias,amplitudes,ch) #3555:3
    rectas(0.07289,0.07995,5820,4881,1,7,frequencias,amplitudes,ch) #3555:3
    rectas2(0.04512,0.04893,7071,7321,1,7,frequencias2,amplitudes2,ch) #3555:3
    rectas2(0.04893,0.05919,7321,7384,1,7,frequencias2,amplitudes2,ch) #3555:3
    expo2(0.05919,0.06446,7384,7132,1,7,frequencias2,amplitudes2,ch) #3555:3
    rectas2(0.06446,0.07189,7132,6633,1,7,frequencias2,amplitudes2,ch) #3555:3
    rectas2(0.07189,0.07859,6633,6070,1,7,frequencias2,amplitudes2,ch) #3555:3
    
    rectas(0.1249,0.13108,7071,7008,1,7,frequencias,amplitudes,ch) #3555:3
    rectas(0.13108,0.1375,7008,6320,1,7,frequencias,amplitudes,ch) #3555:3
    rectas2(0.1249,0.1302,7071,6445,1,7,frequencias2,amplitudes2,ch) #3555:3
    rectas2(0.1302,0.1358,6445,5820,1,7,frequencias2,amplitudes2,ch) #3555:3
    rectas2(0.1358,0.1410,5820,5006,1,7,frequencias2,amplitudes2,ch) #3555:3 
    
    rectas(0.1886,0.1972,6821,6445,1,7,frequencias,amplitudes,ch) #3555:3
    rectas(0.1972,0.2005,6445,6070,1,7,frequencias,amplitudes,ch) #3555:3
    rectas(0.2005,0.2033,6070,5820,1,7,frequencias,amplitudes,ch) #3555:3
    
    senito(0.2754,0.3,5192,1076,0*np.pi/2+0*np.pi,1*np.pi+0*np.pi,1,7,frequencias,amplitudes,ch) #3555:3
    senito2(0.2794,0.2831,6545,300,1*np.pi+0*np.pi,2*np.pi+0*np.pi,1,7,frequencias2,amplitudes2,ch) #3555:3
    senito2(0.2831,0.2937,6570,200,-5*np.pi/2+0*np.pi,7*np.pi/2+0*np.pi,1,7,frequencias2,amplitudes2,ch) #3555:3
    rectas2(0.2937,0.2958,6320,6527,1,7,frequencias2,amplitudes2,ch) #3555:3
  
    
    z = np.polyfit(was, bes, 5) #genera una funcion de las frecuencias en funcion de los betas que ya tengo para poder despues meter otros betas y que me calcule la frecuencia
    p = np.poly1d(z)
    z2 = np.polyfit(was2, bes2, 5) #genera una funcion de las frecuencias en funcion de los betas que ya tengo para poder despues meter otros betas y que me calcule la frecuencia
    p2 = np.poly1d(z2)#se queda solo con un vector de los betas
    asuhg = []
    for i in range(int(tiempo_total/(dt))):
        if(alpha[i]<0):
            beta[i]=p(frequencias[i])
            beta2[i]=p2(frequencias2[i])
            asuhg.append(beta[i])


#%%


    def rk4(dv,v,n,t,dt):
        v1=[]
        k1=[]
        k2=[]
        k3=[]
        k4=[]
        for x in range(0, n):
            v1.append(x)
            k1.append(x)
            k2.append(x)
            k3.append(x)
            k4.append(x)
            
        dt2=dt/2.0
        dt6=dt/6.0
        for x in range(0, n):
            v1[x]=v[x]
        dv(v1, k1)
        for x in range(0, n):
            v1[x]=v[x]+dt2*k1[x]
        dv(v1, k2)     
        for x in range(0, n):
            v1[x]=v[x]+dt2*k2[x]
        dv(v1, k3)
        for x in range(0, n):
            v1[x]=v[x]+dt*k3[x]
        dv(v1, k4)
        for x in range(0, n):
            v1[x]=v[x]+dt*k4[x]        
        for x in range(0, n):
            v[x]=v[x]+dt6*(2.0*(k2[x]+k3[x])+k1[x]+k4[x])
        return v


# una integracion

    n=7 #Cantidad de variables   
    x1=[]
    y1=[]
    z1=[]
    w1=[]
    tiempo1=[]
    sonido=[]
    sonido_total=[]
    amplitud1=[]
    forzado1=[]
    dforzadodt1=[]
    elbeta1=[]
    
    cont1=0
    N=int((L/(350*dt))//1)
    fil1=np.zeros(N)
    back1=np.zeros(N)
    feedback1=0

#while t<tf:
    for i in range(int(tiempo_total/(dt))):
        alp=alpha[i]*(1+0*random.normalvariate(0,0.01))
        b=beta[i]*(1+0*random.normalvariate(0,0.01))
        alp2=alpha[i]*(1+0*random.normalvariate(0,0.01))
        b2=beta2[i]*(1+0*random.normalvariate(0,0.01))        
        uoch=ch[i]
        t=i*dt
        estimulo=fil1[N-1]
        destimulodt=(fil1[N-1]-fil1[N-2])/dt
        rk4(ecuaciones,v,n,t,dt)
        fil1[0]=(v[1]+v[6])*(1+6*random.normalvariate(0,0.05))+back1[N-1]
        back1[0]=-0.85*fil1[N-1]
        fil1[1:]=fil1[:-1]
        back1[1:]=back1[:-1]
        feedback1=back1[N-1]
        x1.append(cont1)  #ACÁ ARMO LOS ARREGLOS DE X Y Z CON LOS RESULTADOS QUE VA LARGANDO "V"
        y1.append(cont1)
        z1.append(cont1)
        w1.append(cont1)
        tiempo1.append(cont1)
        sonido.append(cont1)
        sonido_total.append(cont1)
        amplitud1.append(cont1)
        forzado1.append(cont1)
        dforzadodt1.append(cont1)
        elbeta1.append(cont1)
        x1[cont1]=v[0]
        y1[cont1]=v[1]
        z1[cont1]=v[5]
        w1[cont1]=v[6]
        tiempo1[cont1]=t
        sonido[cont1]=(fil1[0])*envelope(tiempos[i])*(1+2*random.normalvariate(0,0.05))
        sonido_total[cont1]=0
        amplitud1[cont1]=amplitudes[i]
        forzado1[cont1]=estimulo
        dforzadodt1[cont1]=destimulodt
        elbeta1[cont1]=beta[i]

        cont1=cont1+1
        
#%%

    sonido=np.asarray(sonido) 
    f, t, Sxx = signal.spectrogram(sonido, 882000, window=('gaussian',20*128),nperseg=10*1024,noverlap=18*512,scaling='spectrum')
    Sxx = np.clip(Sxx, a_min=np.amax(Sxx)/10**3, a_max=np.amax(Sxx))
    scaled = np.int16(sonido/np.max(np.abs(sonido)) * 32767)
    write('test_juvenil_segunda_alta{}.wav'.format(lazo), 882000, scaled)