Víctor

Simulación de Cuerda preparada y percutida utilizando Scilab y Matlab.

 

  • Eficiencia del software matemático libre en la simulación de instrumentos musicales.
  • En este proyecto se pretende dar un análisis del por que podría ser viable el uso del software libre sobre el comercial (además del ámbito monetario), en la simulación de sistemas dinámicos, como un instrumento musical.

Resultados.

Se ha realizado la simulación de una cuerda percudida y preparada (un objeto metálico interviene sus vibraciones). Está simulación se ha realizado a partir de la solución numérica de la ecuación diferencial parcial que describe el fenómeno descrito. Para realizar los cálculos numéricos se uso el software libre Scilab y el software de pago Matlab.

El corazón del programa y el punto donde la velocidad del calculo se pone a prueba entre estos dos programas , es un loop que realiza 220498 operaciones. Esta cantidad de operaciones es la que se requiere para tener la posibilidad de generar un sonido de 5 segundos de duración.

Como se podrá ver, la sintaxis entre cada uno de los programas no cambia en absoluto. La única diferencia es la inclusión de la linea stacksize. Esta linea se tiene que colocar únicamente en Silab. Con la cual se incrementa el número de elementos que se pueden insertar en una estructura de datos. En este caso las estructuras son vectores.

Los dos algoritmos realizados, no generan sonido, solamente envían los datos que potencialmente al ser leídos a una velocidad de sampleo de audio, pueden producir el sonido simulado. En este punto entra la primera diferencia entre los dos programas; en Scilab no se pudo correr la orden para generar el sonido, mientras que en Matlab dicha orden se implementa de forma muy fácil. Es importante señalar que el comando para generar el sonido es similar en ambos programas.

La gran diferencia entre estos dos programas (además del precio) es la velocidad para realizar los cálculos. Mientras que Matlab tarda 12 segundos en generar 5 segundos de sonido, Scilab tarda 30 segundos.

Velocidad contra precio, es una decisión que no pareciera tan complicada. Pero se hace evidente la necesidad de la velocidad cuando estás probando tus códigos y quieres saber razonablemente rápido si todo va bien.

 

 

 

Programa:Scilab

Tiempo de ejecución: 30 s

 

 

clear all
Fs=44100;
dur=5;
NF = floor(Fs*dur);

L=0.83;
ms=3.93/1000;
ten=59;
b1=0.5;
b3=0.0006;
epsilon=1.6;

Kh=4.5*10^9;
mh=2.3/1000;
p=2.5;
velh= 2.5;
//parametros de tornillo
w=0;
w1=340;
sP=7.0;

C=sqrt((ten/(ms/L)));
k= 1/Fs;
frac=5;
alpha=.12;
betha= (%pi/L);
f0=(sqrt(C^2*betha^2+epsilon^2*betha^4-(b1+b3*betha^2)^2))/(2*%pi);

hu =sqrt(((C^2*k^2)/2+2*b3*k+sqrt((((C^2*k^2)/2+2*b3*k)^2+4*epsilon^2*k^2))));
M=floor(L/hu);
h = L/(M-(1+hu));

mu= (epsilon*k)/(h^2);;
lambda=C*k/h;
v=(2*b3*k)/h^2;

xp=.6;
Xp=floor(xp*M);
iotap=(xp/h-Xp);
X0=round((1-iotap)*Xp+iotap*Xp+1);
xe=floor(alpha*M);
iota=alpha/h-xe;
M0=round((1-iota)*xe+iota*xe+1);

J=1/h;

a3=(2-2*lambda^2-6*mu^2-2*v)/(1+b1*k);
a2= (lambda^2+4*mu^2+v)/(1+b1*k);
a1= (-mu^2)/(1+b1*k);
a4=(-1+2*v+b1*k)/(1+b1*k);
a5 = -v/(1+b1*k);
af=(k^2/(ms/L))/(1+b1*k);

stacksize('max');

y=zeros(M,NF);
yh=zeros(1,NF);
F=zeros(1,NF);

a3p=(2-2*lambda^2-6*mu^2-2*v)/(1+b1*k+J*(k^2*w^2/2+k^2*w1^4*y(X0,1)^2/2+sP*k));
a2p= (lambda^2+4*mu^2+v)/(1+b1*k+J*(k^2*w^2/2+k^2*w1^4*y(X0,1)^3/2+sP*k));
a1p= (-mu^2)/(1+b1*k+J*(k^2*w^2/2+k^2*w1^4*y(X0,1)^2/2+sP*k));
a4p=(-1+2*v+b1*k-J*(k^2*w^2/2+k^2*w1^4*y(X0,1)^2/2-sP*k))/(1+b1*k+J*(k^2*w^2/2+k^2*w1^4*y(X0,1)^2/2+sP*k));
a5p= -v/(1+b1*k+J*(k^2*w^2/2+k^2*w1^4*y(X0,1)^2/2+sP*k));

y(:,1)=0;
y(1,2)=0;
y(M,2)=0; 
yh(2)=velh*k;

y(2:M-1,2)=(y(3:M,1)+y(1:M-2,1))/2;
F(2)=Kh*abs(yh(2)-y(M0,2))^p;

y(1,3)=0;
y(M,3)=0;
 
y(2:M-1,3)=y(1:M-2,2)+y(3:M,2)-y(2:M-1,1);
y(M0,3)=y(M0-1,2)+y(M0+1,2)-y(M0,1)+af*(J*F(2));
y(X0,3)=(a3p*y(X0,2)+a2p*(y(X0+1,2)+y(X0-1,2))+a1p*(y(X0+2,2)+y(X0-2,2))+a4p*y(X0,1)+a5p*(y(X0+1,1)+y(X0-1,1)));
yh(3)=2*yh(2)-yh(1)-((k^2)*F(2)/mh);

F(3)=Kh*abs(yh(3)-y(M0,3))^p;

for n=3:NF

a3p=(2-2*lambda^2-6*mu^2-2*v)/(1+b1*k+J*(k^2*w^2/2+k^2*w1^4*(y(X0,n))^2/2+sP*k));
a2p= (lambda^2+4*mu^2+v)/(1+b1*k+J*(k^2*w^2/2+k^2*w1^4*(y(X0,n))^2/2+sP*k));
a1p= (-mu^2)/(1+b1*k+J*(k^2*w^2/2+k^2*w1^4*(y(X0,n))^2/2+sP*k));
a4p=(-1+2*v+b1*k-J*(k^2*w^2/2+k^2*w1^4*(y(X0,n))^2/2-sP*k))/(1+b1*k+J*(k^2*w^2/2+k^2*w1^4*(y(X0,n))^2/2+sP*k));
a5p= -v/(1+b1*k+J*(k^2*w^2/2+k^2*w1^4*(y(X0,n))^2/2+sP*k));

y(1,n+1)=0;
 y(M,n+1)=0;    
 y(2,n+1)=a1*(y(4,n)-y(2,n))+a2*(y(3,n)+y(1,n))+a3*y(2,n)+a4*y(2,n-1)+a5*(y(3,n-1)+y(1,n-1));
 y(M-1,n+1)=a1*(-y(M-1,n)+y(M-3,n))+a2*(y(M,n)+y(M-2,n))+a3*y(M-1,n)+a4*y(M-1,n-1)+a5*(y(M,n-1)+y(M-2,n-1));
  
 y(3:M-2,n+1)=a3*y(3:M-2,n)+a2*(y(4:M-1,n)+y(2:M-3,n))+a1*(y(5:M,n)+y(1:M-4,n))+a4*y(3:M-2,n-1)+a5*(y(4:M-1,n-1)+y(2:M-3,n-1));
         
 
y(X0,n+1)=(a3p*y(X0,n)+a2p*(y(X0+1,n)+y(X0-1,n))+a1p*(y(X0+2,n)+y(X0-2,n))+a4p*y(X0,n-1)+a5p*(y(X0+1,n-1)+y(X0-1,n-1)));
y(M0,n+1)=a3*y(M0,n)+a2*(y(M0+1,n)+y(M0-1,n))+a1*(y(M0+2,n)+y(M0-2,n))+a4*y(M0,n-1)+a5*(y(M0+1,n-1)+y(M0-1,n-1))+af*(J*F(n)); 
     yh(n+1)=2*yh(n)-yh(n-1)-(k^2)*F(n)/mh;
    
    if (yh(n+1)-y(M0,n+1))>0
        F(n+1)=Kh*abs(yh(n+1)-y(M0,n+1))^p;
    else 
        F(n+1)=0;
    end
    
   
end
 Y=y(M0,:)';

Programa:Matlab

Tiempo de ejecución: 12 s

clear all
Fs=44100;
dur=5;
L=.83 ;
ms=3.93/1000;
ten=59;
b1=0.5;
b3=0.0006;
epsilon= 1.6;

Kh=4.5*10^9;
mh=2.3/1000;
p=2.5;
w=0;
w1=340;
sP=7.0;

C=sqrt((ten/(ms/L)));
k= 1/Fs;
frac=5;
alpha=.12;
betha= (pi/L);
f0=(sqrt(C^2*betha^2+epsilon^2*betha^4-(b1+b3*betha^2)^2))/(2*pi);
gamma=44100/(2*f0);
gamma2 = C;
NF = floor(Fs*dur);
X=sqrt((gamma2^2*k^2+4*b3*k+sqrt(((gamma2^4*k^4+4*b3*k)+16*epsilon^2*k^2)))/(2));

M=floor(L/X);

h = L/(M-(1+X));
R0=sqrt(ten*ms/L);
mu= (epsilon*k)/(h^2);%;
lambda=gamma2*k/h;
v=(2*b3*k)/h^2;
xp=.6;
Xp=floor(xp*M);
iotap=(xp/h-Xp);
X0=round((1-iotap)*Xp+iotap*Xp+1);
xe=floor(alpha*M);
iota=alpha/h-xe;
M0=round((1-iota)*xe+iota*xe+1);
velh= 2.5;
neu=h/C;
J=1/h;
a3=(2-2*lambda^2-6*mu^2-2*v)/(1+b1*k);
a2= (lambda^2+4*mu^2+v)/(1+b1*k);
a1= (-mu^2)/(1+b1*k);
a4=(-1+2*v+b1*k)/(1+b1*k);
a5 = -v/(1+b1*k);
af=(k^2/(ms/L))/(1+b1*k);


lambda^2+4*mu^2

y=zeros(M,NF);
yh=zeros(1,NF);
F=zeros(1,NF);
yP=zeros(1,NF);
FP=zeros(1,NF);
a3p=(2-2*lambda^2-6*mu^2-2*v)/(1+b1*k+J*(k^2*w^2/2+k^2*w1^4*y(X0,1)^2/2+sP*k));
a2p= (lambda^2+4*mu^2+v)/(1+b1*k+J*(k^2*w^2/2+k^2*w1^4*y(X0,1)^3/2+sP*k));
a1p= (-mu^2)/(1+b1*k+J*(k^2*w^2/2+k^2*w1^4*y(X0,1)^2/2+sP*k));
a4p=(-1+2*v+b1*k-J*(k^2*w^2/2+k^2*w1^4*y(X0,1)^2/2-sP*k))/(1+b1*k+J*(k^2*w^2/2+k^2*w1^4*y(X0,1)^2/2+sP*k));
a5p= -v/(1+b1*k+J*(k^2*w^2/2+k^2*w1^4*y(X0,1)^2/2+sP*k));

y(:,1)=0;
y(1,2)=0;
y(M,2)=0; 
yh(2)=velh*k;

y(2:M-1,2)=(y(3:M,1)+y(1:M-2,1))/2;
F(2)=Kh*abs(yh(2)-y(M0,2))^p;

y(1,3)=0;
y(M,3)=0;
 
y(2:M-1,3)=y(1:M-2,2)+y(3:M,2)-y(2:M-1,1);
y(M0,3)=y(M0-1,2)+y(M0+1,2)-y(M0,1)+af*(J*F(2));
y(X0,3)=(a3p*y(X0,2)+a2p*(y(X0+1,2)+y(X0-1,2))+a1p*(y(X0+2,2)+y(X0-2,2))+a4p*y(X0,1)+a5p*(y(X0+1,1)+y(X0-1,1)));
yh(3)=2*yh(2)-yh(1)-((k^2)*F(2)/mh);

F(3)=Kh*abs(yh(3)-y(M0,3))^p;

for n=3:NF
%
a3p=(2-2*lambda^2-6*mu^2-2*v)/(1+b1*k+J*(k^2*w^2/2+k^2*w1^4*(y(X0,n))^2/2+sP*k));
a2p= (lambda^2+4*mu^2+v)/(1+b1*k+J*(k^2*w^2/2+k^2*w1^4*(y(X0,n))^2/2+sP*k));
a1p= (-mu^2)/(1+b1*k+J*(k^2*w^2/2+k^2*w1^4*(y(X0,n))^2/2+sP*k));
a4p=(-1+2*v+b1*k-J*(k^2*w^2/2+k^2*w1^4*(y(X0,n))^2/2-sP*k))/(1+b1*k+J*(k^2*w^2/2+k^2*w1^4*(y(X0,n))^2/2+sP*k));
a5p= -v/(1+b1*k+J*(k^2*w^2/2+k^2*w1^4*(y(X0,n))^2/2+sP*k));

y(1,n+1)=0;
 y(M,n+1)=0; 
 y(2,n+1)=a1*(y(4,n)-y(2,n))+a2*(y(3,n)+y(1,n))+a3*y(2,n)+a4*y(2,n-1)+a5*(y(3,n-1)+y(1,n-1));
 y(M-1,n+1)=a1*(-y(M-1,n)+y(M-3,n))+a2*(y(M,n)+y(M-2,n))+a3*y(M-1,n)+a4*y(M-1,n-1)+a5*(y(M,n-1)+y(M-2,n-1));
 
 y(3:M-2,n+1)=a3*y(3:M-2,n)+a2*(y(4:M-1,n)+y(2:M-3,n))+a1*(y(5:M,n)+y(1:M-4,n))+a4*y(3:M-2,n-1)+a5*(y(4:M-1,n-1)+y(2:M-3,n-1));
 
 
y(X0,n+1)=(a3p*y(X0,n)+a2p*(y(X0+1,n)+y(X0-1,n))+a1p*(y(X0+2,n)+y(X0-2,n))+a4p*y(X0,n-1)+a5p*(y(X0+1,n-1)+y(X0-1,n-1)));
y(M0,n+1)=a3*y(M0,n)+a2*(y(M0+1,n)+y(M0-1,n))+a1*(y(M0+2,n)+y(M0-2,n))+a4*y(M0,n-1)+a5*(y(M0+1,n-1)+y(M0-1,n-1))+af*(J*F(n)); 
 yh(n+1)=2*yh(n)-yh(n-1)-(k^2)*F(n)/mh;
 
 if (yh(n+1)-y(M0,n+1))>0
 F(n+1)=Kh*abs(yh(n+1)-y(M0,n+1))^p;
 else 
 F(n+1)=0;
 end
 
 
end
 Y=y(M0,:)';

Responder

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s