Leer, graficar y reproducir archivos de sonido

La siguiente rutina lee un archivo de audio de nombre diapazon.wav, realiza una gráfica teniendo en cuenta la frecuencia de muestreo Fs  con la que fue registrada y luego reproduce el sonido.

Descargue el archivo de sonido aquí. Lógicamente, usted puede crear su propio archivo de sonido.


%---INICIO
close all
clc
clear all
%--(1). Lectura
[FILE,RUTA] = uigetfile('*.wav'); %para escoger la ruta del archivo a leer
NOM=[RUTA,'/',FILE];
[y, Fs] = audioread(NOM); %lee el archivo
% Fs-->frecuencia de muestreo del microfono utilizado (muestras por segundo)
%  y-->datos registrados
%--(2). Gráfica del sonido (oscilación temporal)
dt= 1/Fs; %--tiempo entre cada muestra
L=length(y); %--numero de muestras registradas
t=dt*linspace(0,L-1,L); %---dominio del tiempo
f=figure;
plot(t,y)
grid on
%--(3). Reproducción
h = uicontrol('Position',[300 100 200 40],'String','Reproducir',...
              'Callback','uiresume(gcbf)');%--crea el botón para reproducir
uiwait(gcf); % espera hasta que se presione el botón
delete(gco)
sound(y,Fs) % reproduce el sonido
%-------FIN

RESULTADO

Figura 1. Oscilación temporal registrada de la onda de sonido

 

Figura 2. Zoom de la figura 1 en los primeros 10 milisegundos

Métodos para ecuaciones diferenciales de primer orden: Metodo de Euler, Euler mejorado Y Runge-Kutta 4° orden

Con esta entrada comenzaremos con el uso de métodos numéricos para resolver ecuaciones diferenciales. Comenzaremos con los sencillos métodos de Euler y Euler mejorado para la solución de ecuaciones diferenciales de primer orden con valor inicial y finalizaremos con el método de Runge-Kutta de cuarto orden.

 Ecuación diferencial de primer orden con valor inicial:

 y' = f(x,y),         y(xo) = yo.

Utilizaremos los métodos de Euler para resolver la ecuación diferencial
   
y' = y + x,      y(0) = 0.

Esta ecuación tiene solución analítica    y  =  -1 + exp(x) - x.

A modo de comparación graficaremos la solución encontrada con el método numérico y la solución analítica.
 --------------------------------------
1). METODO DE EULER
 ---------------------------------------

yn+1 = yn + hf(xn, yn);              n = 0,1,2,3...,N 

El programa realizado en MatLab es el siguiente:

%-----INICIO-----
clc
close all
clear all
h = 0.2; %---paso
N=10;    %---numero de valores a calcular
%----valores iniciales
xo = 0;
yo = 0;
%---Definición del vector "y" de longitud N e inicialmente con ceros.
y = zeros(1,N);
y(1)=yo; %--valor inicial asignado al primer elemento de y.
yn = yo;
for n=1:N-1;
    xn = xo+(n-1)*h;
    f = (yn+xn); %--funcion f que depende de la ODE a resolver
    y(n+1) = yn + h*f;
    yn = y(n+1);
end
%---graficas
x = 0:h:h*(N-1);
ya = -1 + exp(x)-x; %--solución analítica
plot(x,ya,'k',x,y,'r')
hold on
plot(x,ya,'k*',x,y,'ro')
xlabel('x')
ylabel('y')
legend('analytical solution','Euler Method')
grid on
%---FIN

Resultado

 
 ----------------------------------------------------------
2). METODO DE EULER MEJORADO
 ----------------------------------------------------------
 y*n+1 = yn + hf(xn + yn)

yn+1 = yn + (1/2)h[f(xn, yn) + f(xn+1 , y*n+1)];              n = 0,1,2,3...,N 

El programa realizado en MatLab es el siguiente:

%-----INICIO-----
clc
close all
clear all
h = 0.2;
N=10;  
xo = 0;
yo = 0;
y = zeros(1,N);
y(1)=yo;
y1 = yo;
for m=1:N-1;
    x1 = xo+(m-1)*h;
    x2 =x1 +h;
    f1 = y1+x1;
    y2 = y1+h*f1;
    f2 = y2+x2;
    y(m+1) = y1 + 0.5*h*(f1+f2);
    y1 = y(m+1);
end

 x = 0:h:h*(N-1);
ya = -1 + exp(x)-x;

figure
plot(x,ya,'k',x,y,'r')
hold on
plot(x,ya,x,y,'*r')
legend('Analytical solution', 'Improved Euler Method')
xlabel('x')
ylabel('y')
grid on
%---END

Resultado



 ---------------------------------------------------
2). MÉTODO DE RUNGE-KUTTA
 ---------------------------------------------------

 k1 = hf(xn, yn)
 k2 = hf(xn + 0.5h, yn + 0.5k1)
 k3 = hf(xn + 0.5h, yn + 0.5k2)
 k4 = hf(xn + h, yn + k3)
 xn+1 = xn+h
 yn+1 = yn + (1/6)(k1 + 2k2 + 2k3 + k4)
 n = 0,1,2,3...,N 

El programa realizado en MatLab es el siguiente:

yr = zeros(1,N);
yr(1)=yo;
y1 = yo;
for r=1:N-1;
    x1 = xo+(r-1)*h;
    k1 = h*(x1+y1);
    k2 = h*(x1+0.5*h+y1+0.5*k1);
    k3 = h*(x1+0.5*h+y1+0.5*k2);
    k4 = h*(x1+h+y1+k3);
    yr(r+1) = y1 + (1/6)*(k1+2*k2+2*k3+k4);
    y1 = yr(r+1);
end 
x = 0:h:h*(N-1);
ya = -1 + exp(x)-x;
figure
plot(x,ya,'k',x,yr,'r')
hold on
plot(x,ya,'k.',x,yr,'ro')
legend('Analytical solution', 'Runge-Kutta Method')
grid on
xlabel('x')
ylabel('y')
Resultado


Comparación gráfica:  Solución analítica, método de Euler, Método de Euler mejorado Y método de Runge-Kutta.

Solución análítica de ecuaciones diferenciales ordinarias

En esta entrada usaremos MatLab para encontrar la solución analítica de una ecuación diferencial ordinaria. Para ello utilizaremos la función dsolve.
--------------------------------------------------------------------------------------------------------------------------------
dsolve
Soluciona ecuaciones diferenciales ordinarias y sistemas de ecuaciones diferenciales ordinarias
Sintaxis
S = dsolve(eqn)
S = dsolve(eqn,cond)
S = dsolve(eqn,cond,Name,Value)
Y = dsolve(eqns)
Y = dsolve(eqns,conds)
Y = dsolve(eqns,conds,Name,Value)
[y1,...,yN] = dsolve(eqns)
[y1,...,yN] = dsolve(eqns,conds)
[y1,...,yN] = dsolve(eqns,conds,Name,Value)
Descripción 
En esta entrada solo describiremos las dos primeras.
1). S = dsolve(eqn): Soluciona la ecuación diferencial ordiaria eqn. Aquí eqn es una ecuación simbólica que contiene diff para indicar derivadas.
 2). S = dsolve(eqn,cond): Soluciona la ecuación diferencial ordinaria eqn con la condición inicial o condición de contorno cond.

Las variables simbólicas se pueden definir con la función syms. Alternativamente, se puede usar un string con la letra D para indicar derivadas.
----------------------------------------------------------------------------------------------------------------------------- 
Ejemplo 1: Obtengamos la solución particular de la ecuación diferencial ordinaria de primer orden:

                            y'= sin(2*x)-y*tan(x);   y' = dy/dx,
con la condición inicial
                            y(0)=1.
----------------INICIO------------------------------------------------------------------------------------------------------

%---creamos las variables simbólicas con syms:
syms y(x)
%----Construimos la ecuación diferencial y la colocamos en la función dsolve:
y(x) = dsolve(diff(y,x) == sin(2*x)-y*tan(x),y(0)==1);
%----Evaluamos y para valores de x entre 0 y 4 a pasos de 0.1 y los valores lo asignamos a la variable y1: 
x = 0:0.1:4;
y1=eval(y);
%----graficamos y1 usando la función plot:
figure
plot(x,y1)
grid on
%---Evaluamos y para el valor x = 1 usando eval: 
 x=1;
yo=eval(y);
%---graficamos el punto  (x,yo) sobre la gráfica de y(x), de color rojo 'r' y marca 'o' :
hold on
plot(x,yo,'ro')
xlabel('x')
ylabel('y')
grid on

-------------------FIN-----------------------------------------------------------------------------------------------------------

Los resultados se muestran a continuación:
La solución es:
y(x) = 3*cos(x)-cos(2x)-1
Evaluada para x=0 resulta:
yo = 1.0371

La gráfica resultante es:


------------------------------------------------------------------------------------------------------------------------------
Ejemplo 2: Obtengamos la solución particular de la ecuación diferencial ordinaria de segundo orden:
                          h'' + 2h'+2h = 0;  h' = dh/dt,
con la condición inicial
                           h(0)=0  ;  h''(0)=15 .
----------INICIO-------------------------------------------------------------------------------------------------------
%---creamos las variables simbólicas con syms:
syms h(t)
%----Construimos la ecuación diferencial y la colocamos en la función dsolve. Usamos D en lugar de diff para definir la ODE colocándola como string (entre comillas simples):
h(t) = dsolve('D2h == -2*Dh-2*h', h(0)==0, 'Dh(0)==15')
%----Evaluamos h para 100 valores de t entre 0 y 10 y lo asignamos a la variable h1:
N=100;
t = linspace(0,10,N);
h1=eval(h);
%----graficamos h1 usando la función plot:
figure
plot(t,h1)
grid on
%---Evaluamos h para el valor t = 2:
t=2;
ho=eval(h);
%---graficamos el punto (0,ho) sobre la gráfica de h(t), de color rojo 'r' y marca '*' :
hold on
plot(t,ho,'r*')
xlabel('t')
ylabel('h')
grid on
-------------FIN-------------------------------------------------------------------------------------------------------------

Los resultados se muestran a continuación:
La solución es:
h(t) = 15*exp(-t)*sin(t)
Evaluada para t=2 resulta:
ho = 1.8459
La grafica resultante es:

Gráficas de campos direccionales con MatLab

Esta rutina en Matlab realiza la gráfica de un campo direccional a partir de la ecuación diferencial :

y'=f(x,y)=-5*x.^4.*y.^2.

Sobre el campo de direcciones se realiza la gráfica de la solución particular:

y = 1./(x.^5+1).

%=====INICIO========
%-1)---Región--------
x1=0;
x2=2;
y1=0;
y2=2;
N = 20;
x = linspace(x1,x2,N);
y = linspace(y1,y2,N);
[X,Y]= meshgrid(x,y);

%-2)---calculo del angulo de inclinación de cada vector
m = -5*X.^4.*Y.^2;
T = atan(m);

%-3)--componentes de cada vector
U = cos(T);
V = sin(T);

%-4)--Dibujo de los vectores tangentes
h = quiver(X,Y,U,V);

%-5)--superposición en el campo de direcciones de la solución particular
hold on
yp = 1./(x.^5+1);
plot(x,yp)
grid on
%===FIN========

El resultado es el siguiente:

Fig.1. Campo de direcciones (azul). Solución particular (rojo)