Autor: new.jmanza@gmail.com

Grabando un archivo wav

[\atexpage]

https://www.dropbox.com/scl/fi/6z98nsbu3wuzfffaqjqjf/g.wav?rlkey=ilxj1ygpzawlks88x9tvpnw8k&dl=0

Para leer el arhcivo anterior, usamos

import numpy as np
import matplotlib.pyplot as plt
from scipy.io import wavfile

# Leer el archivo WAV
sample_rate, data = wavfile.read('g.wav')

# Comprobar si el archivo es estereofónico o monofónico
if len(data.shape) == 2:
    # Si es estereofónico, tomamos solo uno de los canales
    data = data[:, 0]

# Crear un vector de tiempo
time = np.linspace(0, len(data) / sample_rate, num=len(data))

# Graficar la señal
plt.figure(figsize=(12, 6))
plt.plot(time, data)
plt.title('Señal acústica')
plt.xlabel('Tiempo (s)')
plt.ylabel('Amplitud')
plt.grid()
plt.show()

Para obervar que tenemos un seno, hacemos un zoom a la grafica

La instruccion sample_rate define la frecuencia a la cual se toman los datos, el periodo con el cual se toman es $\tau_d$ = 1/sample_rate=2.08e-5 seg, en la sigueinte grafica

doblez3 y la FT como funcion de H

import meep as mp
import numpy as  np
import math
from matplotlib import pyplot as plt

############## -> H
NH = 220
Hi = 0.0
Hf = 11
dH = (Hf-Hi)/(NH)
HV =  np.zeros(NH)
maxV = np.zeros(NH)
QQxV = np.zeros(NH)
maxV = np.zeros(NH)
############# <- H

############# -> FT
i_0=500
i_I=800
I=i_I-i_0
Qx_i=0.01
Qx_f=1.0
NQx=100
dQx=(Qx_f-Qx_i)/NQx
EQx=np.zeros(NQx)
QxV=np.zeros(NQx)
mapaFT=np.zeros((NH,NQx))
############ <- FT


for iH in range(NH):
 H = Hi+iH*dH
 HV[iH] = H
 print(iH,H)
    
############################################### -> meep 
 L = 40
 t_until=400
 cell = mp.Vector3(2*L, L, 0)
 geometry = [mp.Block(mp.Vector3(L,  1, mp.inf),center=mp.Vector3(-L/2,-H/2),material=mp.Medium(epsilon=3.6*3.6),e1=[1.0, 0.0],e2=[0.0, 1.0]),
             mp.Block(mp.Vector3(1,H+1, mp.inf),center=mp.Vector3(   0,   0),material=mp.Medium(epsilon=3.6*3.6),e1=[1.0, 0.0],e2=[0.0, 1.0]),
             mp.Block(mp.Vector3(L,  1, mp.inf),center=mp.Vector3(+L/2,+H/2),material=mp.Medium(epsilon=3.6*3.6),e1=[1.0, 0.0],e2=[0.0, 1.0])]
 sources = [mp.Source(mp.ContinuousSource(frequency=0.1), component=mp.Ez, center=mp.Vector3(-L+2, -H/2))]
 pml_layers = [mp.PML(1.0)]
 resolution = 10
 sim = mp.Simulation(cell_size=cell,boundary_layers=pml_layers,geometry=geometry,sources=sources,resolution=resolution,)
 print('----------------- H = ',H)
 plt.figure(dpi=100)
 sim.plot2D()
 plt.savefig('a'+str(H)+'.png')
 plt.show()
 sim.run(until=t_until)
 plt.figure(dpi=100)
 sim.plot2D(fields=mp.Ez)
 plt.savefig('b'+str(H)+'.png')
 plt.show()
 eps_data = sim.get_array(center=mp.Vector3(), size=cell, component=mp.Dielectric)
 ez_data = sim.get_array(center=mp.Vector3(), size=cell, component=mp.Ez)
 plt.imshow(eps_data)
 plt.savefig('c'+str(H)+'.png')
 plt.show()
 plt.imshow(ez_data)
 plt.savefig('d'+str(H)+'.png')
 plt.show()
 plt.plot(ez_data[:,int((L/2+H/2)*10)])
 plt.savefig('e'+str(H)+'.png')
 plt.show()
################################################## <- meep


################################################## -> FT
 for iQx in range(NQx):
  Qx = Qx_i+iQx*dQx
  QxV[iQx]=Qx
  Exr=0
  Exi=0
  for i in range(I):
   Exr=Exr+ez_data[i_0+i,int((L/2+H/2)*10)]*math.cos(2.0*math.pi*Qx*((i_0+i)/resolution))
   Exi=Exi+ez_data[i_0+i,int((L/2+H/2)*10)]*math.sin(2.0*math.pi*Qx*((i_0+i)/resolution))
  EQx[iQx] = math.sqrt(Exr*Exr+Exi*Exi)  
  mapaFT[iH,iQx]=EQx[iQx]
 plt.plot(QxV,EQx)
 plt.savefig('f'+str(H)+'.png')
 plt.show()
 maxV[iH] =  np.max(EQx)
 el_indice = np.argmax(EQx)
 QQxV[iH] = QxV[el_indice]
#################################################### <- FT

plt.plot(HV,maxV,'-o')


plt.xlabel('H', fontsize=14)
plt.ylabel('max[E(Q_x=2.57]', fontsize=14)


# Mostrar la gráfica
plt.grid()
plt.show()

Entendiendo las tripas del meep: doblez2 y la FT

import meep as mp

L = 40
H = 1
t_until=400
cell = mp.Vector3(2*L, L, 0)
geometry = [mp.Block(mp.Vector3(L,  1, mp.inf),center=mp.Vector3(-L/2,-H/2),material=mp.Medium(epsilon=3.6*3.6),e1=[1.0, 0.0],e2=[0.0, 1.0]),
            mp.Block(mp.Vector3(1,H+1, mp.inf),center=mp.Vector3(   0,   0),material=mp.Medium(epsilon=3.6*3.6),e1=[1.0, 0.0],e2=[0.0, 1.0]),
            mp.Block(mp.Vector3(L,  1, mp.inf),center=mp.Vector3(+L/2,+H/2),material=mp.Medium(epsilon=3.6*3.6),e1=[1.0, 0.0],e2=[0.0, 1.0])]
sources = [mp.Source(mp.ContinuousSource(frequency=0.1), component=mp.Ez, center=mp.Vector3(-L+2, -H/2))]
pml_layers = [mp.PML(1.0)]
resolution = 10
sim = mp.Simulation(cell_size=cell,boundary_layers=pml_layers,geometry=geometry,sources=sources,resolution=resolution,)


from matplotlib import pyplot as plt
plt.figure(dpi=100)
sim.plot2D()
plt.show()

sim.run(until=t_until)
plt.figure(dpi=100)
sim.plot2D(fields=mp.Ez)
plt.show()

eps_data = sim.get_array(center=mp.Vector3(), size=cell, component=mp.Dielectric)
ez_data = sim.get_array(center=mp.Vector3(), size=cell, component=mp.Ez)

print(ez_data.shape)
plt.imshow(ez_data)
plt.plot(ez_data[:,200],'-k')
plt.plot(ez_data[:,210],'-b')
plt.plot(ez_data[:,220],'-g')
plt.plot(ez_data[:,230],'-r')
import numpy as np
import math

############ FT
i_0=500
i_I=800
I=i_I-i_0
Qx_i=0.0
Qx_f=1.0
NQx=1000
dQx=(Qx_f-Qx_i)/NQx
EQx=np.zeros(NQx)
QxV=np.zeros(NQx)
############# FT
for ic in range(NQx):
  Qx = Qx_i+ic*dQx
  QxV[ic]=Qx
  Exr=0
  Exi=0
  for i in range(I):
   Exr=Exr+ez_data[i_0+i,int((L/2+H/2)*10)]*math.cos(2.0*math.pi*Qx*(i_0+i)/resolution)
   Exi=Exi+ez_data[i_0+i,int((L/2+H/2)*10)]*math.sin(2.0*math.pi*Qx*(i_0+i)/resolution)
  EQx[ic] = math.sqrt(Exr*Exr+Exi*Exi)
plt.plot(QxV,EQx)



# Graficar


# Personalización de la gráfica
plt.title('Gráfica de la función y = x', fontsize=16)
plt.xlabel('x', fontsize=14)
plt.ylabel('y', fontsize=14)
plt.axhline(0, color='black',linewidth=0.5, ls='--')  # Línea horizontal en y=0
plt.axvline(0, color='black',linewidth=0.5, ls='--')  # Línea vertical en x=0
plt.grid(color = 'gray', linestyle = '--', linewidth = 0.5)
# Mostrar la gráfica
plt.show()