RocoBlog

El Motor paso paso temperado a ±0.5[Hz] ( o Música con motores eléctricos)

2010-04-21T01:21:00.008-04:00

El título hace referencia al "Clavecín bien temperado". Como nunca se pueden temperar o afinar bien las cosas, se indicó el máximo error por resolución en la afinación, que corresponde a medio [Hz].

Como ejemplo, a continuación se trata de reproducir, a dos voces - que es lo que permite el controlador- el Himno Nacional de Chile, para acentuar el nacionalismo post terremoto, la Badinerie de J.S. Bach, pues a el se le debe el título y el Himno de la Universidad de Concepción mi alma mater y dueña de los motores.

Un intento de explicación

Cuando se puede controlar la frecuencia con la resolución de ±1[Hz] en una máquina y se puede controlar el tiempo de funcionamiento a esa frecuencia, se puede hacer música con ella.

Este es el caso de los motores paso paso o "stepper motors". En el caso particular del ejemplo de más abajo, usando un controlador ISEL basta con enviarle los movimientos vía puerto serial.

Usando la fórmula para convertir notas musicales en frecuencia y modificando el programa para hacer música en GNU-Octave/Matlab, de manera que entregue la sintaxis de los movimientos de los motores en vez de la salida del audio, es una opción para hacer música con motores paso paso.

Es importante decir que la sintaxis típica del motor es:

[np motor 1],[vel motor 1],[np motor 2],[vel motor 2]

por lo que el programa debe hacer la conversión de tiempo a numero de pasos basado en la frecuencia de giro. Existe un problema en este punto, porque para bajas frecuencias el error en tiempo crece, y se desea que para cada línea el tiempo de ejecución sea el mismo, porque por cada línea se da la instrucción de mover 2 motores simultáneamente.

El programa para la conversión, está probado en GNU-Octave y requiere de la función notacnc_2.m, similar a la del programa anterior, para funcionar. En la sintaxis del programa se adjuntan las mismas canciones, expuestas anteriormente, como el Badinerie.

El Código

El programa

La función

En resumen, si puedes controlar la frecuencia, puedes hacer música.

¡Fuerza Chile!

2010-03-05T13:22:00.006-03:00

Tras el Terremoto de 8,8 grados en la escala Richter que ha afectado a la zona centro sur de Chile se necesita ayuda

Manual de Cómo ayudar

Alternativas de Ayuda

Chile Ayuda a Chile

Cuenta Única 2702 Banco de Chile/Banco Santander

La foto que representa la campaña Chile ayuda a Chile, que muestra a un compatriota levantando una bandera chilena embarrada y rota entre escombros. Ese simbolismo contagioso me ha llevado a utilizar el parlante interno de mi PC para levantar simbólicamente nuestra bandera musical, el Himno nacional de Chile.

El himno está arreglado para dos voces, en el video aparece una voz grabada previamente y superpuesta, para que suene mejor.

El código del programa, el ejecutable y las notas de la primera y segunda voz del himno están para la descarga vía megaupload y vía rapidshare

A modo de comentario es importante decir que los terremotos no son un acontecimiento aislado en Chile. De hecho la misma zona ha sentido terremotos en el periodo histórico de la colonia. En 1751 uno destruyó Concepción y Chillán, fue en ese entonces, que la ciudad de Concepción se re-fundó en el lugar que está hoy -antes estaba en la localidad de Penco (Por eso el gentilicio de Concepción es penquista)- y la Ciudad de Chillán se reconstruyó donde hoy está Chillán Viejo. El Chillán actual, que está al lado de Chillán Viejo, fue construido por otro terremoto, en 1835. Y sufriría otros terremotos en 1939, 1953 y 1960.

Por eso ánimo a todos, fuerza Chile, te levantarás.

Francisco Roco

Cómo hacer música con el parlante interno del PC y C++

2009-12-14T04:41:00.005-03:00

Al pensar en el parlante o altavoz interno del PC, en inglés, PC speaker, se nos viene a la cabeza el bip, en inglés beep que suena como aviso.

Usando la consola en cualquier distribución basada en GNU/Linux se puede controlar el parlante, de modo de hacer frecuencias y duraciones a voluntad y poder hacer música. A diferencia del uso de GNU-Octave y Matlab para hacer música, noi se necesita construir la onda, solo se debe controlar la frecuencia y duración.

Tomando como referencia la entrada de notas usada para hacer música en Matlab, en la que se ingresan notas de forma parecida a como se hace en los programas compositores en los celulares Nokia, especificando la duración en el primer caracter, la nota en clave americana, y la octava tomando como octava central la 3ra, se puede hacer una función teniendo como entrada las notas en dicho formato y como salida las frecuencias y las duraciones de esas notas.

La salida de la función -las frecuencias y las duraciones pueden ser ingresadas a otros dispositivos, en el caso del Matlab era a otra función que creaba las ondas, a motores controlados electronicamente, etcétera.

En este caso nos comunicamos con el parlante enviandole frecuencias y notas:

//  La Función que hace los sonidos
int suena(int freq,int len){
  
 int fd = open("/dev/console", O_WRONLY);
      ioctl(fd, KIOCSOUND, (int)(1193180/freq));
      //usleep(1000*len);

 usleep(1000*len);
      ioctl(fd, KIOCSOUND, 0);
      close(fd);
}

Esta parte puede hacerse de forma distinta instalando gnubeep. Con el que se puede escribir

gnubeep -f 2600 -l 1000000

y sonarán 600 [Hz] por 1 segundo.

El código

El código escrito en C++ funcionó en mi PC con Mandriva 2009.1, y probablemente funcione en cualquier distribución GNU/Linux.

Puede ser compilado para Windows descomentando las partes indicadas con //win .

Está a disposición de vustra merced en el enlace: musica.cpp.

Algunos ejemplos para ingresarle al programa:

Badinerie BWV 1067 de J.S. Bach. La transcripción es para la primera voz y la parte del bajo
- Badinerie primera voz
- Badinerie segunda voz

Suite para Cello N°1, BWV 1007 J.S. Bach
- Suite de Cello

Invention 13 (BWV 784) de J.S. Bach
- Invención 13 primera voz
- Invención 13 segunda voz

Ejemplos

Usando dos computadores y sincronizando el instante en que empiezan las canciones se probaron los ejemplos anteriores.

Badinerie BWV 1067

Invención 13, BWV 784

Suite para Cello 1007

Comentario

Si quiere desactivar el pitido, escriba en la consola como superusuario:

# modprobe -r pcspkr

Francisco Roco

Como hacer musica en GNU Octave o Música en GNU Octave o bien migrando "Musica en Matlab" a Octave

2009-09-10T23:53:00.017-04:00

En la una de las entradas anteriores se describe un programa para hacer música en Matlab, programa que fue probado en Windows.

La idea es si existe un computador con alguna distribución de GNU/Linux, y en vez de Matlab ocupamos el siempre útil y nunca bien ponderado GNU Octave para ejecutar este programa.

Muchos comandos de Matlab se ejecutan de igual manera en GNU Octave.
Excepto por algunos comandos importantes :

menu('title','opc 1','opc 2')

Este comando no está en GNU Octave, por lo menos en la instalación del paquete Octave-forge, así que el menú con botones de Matlab, pasa a ser un sencillo menú por consola.

Otro comando que existe para Matlab y funciona en Octave es :

sound(y,fm)

las primeras veces que ejecuté el programa arrojaba un error así:


octave:2> Fs=44100;
octave:3> t=0:1/Fs:1;
octave:4> y=sin(2*pi*440*t);
octave:5> sound(y,Fs)
sh: ofsndplay: command not found
warning: warn_fortran_indexing is no longer a built-in variable; please read the NEWS file or type `news' for details
warning: warn_fortran_indexing is no longer a built-in variable; please read the NEWS file or type `news' for details
warning: broken pipe -- some output may be lost

La solución es simplemente un

clear all

antes de ejecutar el comando. Existe un error al hacer esto cuando un instala el Octave-forge, pero la solución es borrar unos paquetes que por lo menos en este caso no son prioridad, en detalle está acá.

El funcionamiento muestra el detalle de lo reproducido:



octave:7> clear all
octave:8> Fs=44100;
octave:9> t=0:1/Fs:1;
octave:10> y=sin(2*pi*440*t);
octave:11> sound(y,Fs)
warning: warn_fortran_indexing is no longer a built-in variable; please read the NEWS file or type `news' for details
warning: warn_fortran_indexing is no longer a built-in variable; please read the NEWS file or type `news' for details
play au: header size 24 is too small

-: (au)

Encoding: Signed PCM
Channels: 1 @ 16-bit
Samplerate: 44100Hz
Replaygain: off
Duration: unknown

In:0.00% 00:00:01.00 [00:00:00.00] Out:44.1k [      |      ] Hd:0.0 Clip:0
Done.

Al intentar correr el programa en Fedora 12 obtenía el siguiente error:

octave:1> Fs=44100;t=0:1/Fs:1;y=sin(2*pi*440*t);sound(y)
error: sound.m: No command line utility found for sound playing
error: called from:
error:   /usr/share/octave/packages/audio-1.1.4/sound.m at line 108, column 5

La solución para esto, que dan en los foros de ubuntu es en resumen la siguiente:

Tener instalado alsa-utils o pulseaudio-utils.
Luego dice poner en ~/.octaverc

global sound_play_utility = 'aplay';

o bien

global sound_play_utility = 'paplay';

o bien ponerlo directamente en Octave.

Vemos que funciona:

octave:2>  Fs=44100;t=0:1/Fs:1;y=sin(2*pi*440*t);sound(y,Fs)
Playing Sparc Audio 'stdin' : Signed 16 bit Big Endian, Rate 44100 Hz,Mono

Otro comando que se ejecuta de manera parecida en Octave es:


wavwrite( 'nombre',y, Fs, nbits)

la diferencia es importante es que la señal guardada en un vector y, debe ser ingresada necesariamente como un vector columna. En este caso y es un vector fila, cual debe ser ingresado traspuesto:


octave:13> wavwrite( 'Roco_audio.wav',y', Fs, 16)

Con esas consideraciones, el programa esta listo para ejecutarse.
Básicamente, los cambios que hay que hacerle al código Matlab para ejecutarse en Octave son relativamente menores.

El programa fue probado en GNU Octave, version 3.0.3 sobre Mandriva Linux Free 2009 Spring.

El programa

El menú del programa se ve así, obviamente las tildes son omitidas para evitar errores de codificación de caracteres :


   MENU     
-----------------------
Seleccione Opcion

1) Ingresar notas 1ra voz
2) Ingresar notas 2da voz
3) Ingresar notas 3ra voz
4) Ingresar Percusion
5) Sonido
6) Guardar Cancion
7) Duracion de la negra [s]
8) Salir
-----------------------

La idea es la misma que en aca, salvo que en Octave.
Acá el programa y las funciones

y un ejemplo que no incluye el uso de la percusión :

Versión de un fragmento de la canción "Rosas" de la Oreja de van Gogh GNU Octave

La Badinerie de J.S. Bach BWV 1067

Chao

Francisco Roco

Gráficos en GNU Octave, la alternativa a Matlab

2009-05-30T00:11:00.006-04:00

En GNU Octave el programa para cálculos numéricos similar a Matlab, la mayoría de los comandos en la sintaxis son iguales.

Para graficar, se ocupa el comando plot(x,y) el cual entrega gráficos menos agradables a la vista que los logrados en Matlab. Octave ocupa para esto otro Software que se llama GNU Plot.

Se puede cambiar el uso del GNU Plot por otro programa, que mejora la presentación de los gráficos, usando comandos iguales en algunos casos, y otros parecidos.

Ese programa es OctPlot.

Para obtener octplot + Octave:

$ sudo apt-get install octave octplot

Luego para activar OctPlot, despues de hechar a correr Octave se pone:

octave:1> toggle_octplot

Lo que desactiva GNUplot y deja trabajando a OctPlot.

Un Ejemplo de un gráfico obtenido en Octplot, a continuación la sintaxis:

t=0:.1:10;
y=sin(pi*t);
x=cos(pi*t);
z=2*cos(pi*t+pi/4)+1;
plot(t,y,'k')
hold on
plot(t,x,'k--')
plot(t,z,'r.-')

title('Titulo del Grafico');
xlabel('Variable Independiente');
ylabel('Vatiable Dependiente');

axis([0 10 -1 3]);
legend('Signal 1','Signal 2','signal 3');
text(5,2,'Texto en el grafico');
print('MiGrafico.png', "-dpng");

Chao.

Francisco Roco

Órbita de Molniya o Molniya Orbit o молния

2008-05-19T23:07:00.014-04:00

La órbita de Molniya es una órbita inventada por los camaradas de la ex Unión Soviética para utilizarla en sus satélites de comunicación. Antes de hacer esto, pensé que era fácil comunicar vía satélite a cualquier parte del globo, pero no es así, pues sólo para lugares que estan cercanos al Ecuador se puede poner un satélite geoestacionario.

¿Que pasa con los lugares que están cercanos a los polos?

El satélite geoestacionario permenece sobre el punto que se requiere porque mantiene su velocidad igual a la de este punto. ¿Que características debe tener la órbita para cumplir con esta misión?

La velocidad de rotación de la tierra es constante, y si quiero que el satélite siga a un punto en la superfice la velocidad de rotación - o velocidad angular - del satélite debe ser la misma. La órbita con la velocidad angular constante es la circular, que usan los geoestacionarios, en ella, la distancia del satélite a la tierra permanece igual siempre -no como el la órbita elíptica-.

La órbita del satelite estará sobre un plano, este plano contendrá siempre el centro de la tierra, es decir , si quiero que mi satélite geostacionario pase por Concepción, defino el plano por el que esta esta ciudad, y me doy cuenta que al girar sobre este punto el satélite no pasará por el centro de la tierra, es decir no puedo poner un satelite geoestacionario ahí, ¿Dónde puedo? En el ecuador porque el plano de las ubicaciones que están ahi contiene el centro de la tierra.

En ese momento es cuando el DT llama para que entre a la cancha a la órbita elíptica. La característica que tiene esta órbita es que tiene una distancia máxima, el apogeo, y una mínima, el perigeo. La gracia esta en que en el perigeo la velocidad del satélite

La que Esta órbita es una buena solución para el problema de cubrir zonas de alta latitud, dicho de otra forma , de zonas cercanas a los polos. Los camarradas pensarón en esta órbita para su Madrrre RRRusia, pero también puede servir pa nuestro Chile querido porque tenemos territorio muy al sur del mundo.

El objetivo de las siguientes lineas de sintaxis es reproducir por cuales zonas sobre el globo pasa nuestro satélite y de qué zonas cubre, o sea que zonas en ese instante es capaz de dar comunicación.

Hay que decir que esta órbita tiene 63.4° de inclinación respecto del Ecuador, y tiene su apogeo cerca de Ciudad de Mejico, y algún lugar de Rusia.

Lo notable de esto es que se hizo en tiempos donde no existían metodos como Matlab u Octave para hacer todos estos gráficos y cálculos, los camarradas la sacaron con lápiz, regla de cálculo y mucho ingenio. Debo confesar que esto me costo y si hubiese sido uno de esos rusos probablemente hubiese terminado en el Gulag por no presetar cálculos a tiempo.

Sintaxis en MATLAB

clc;
clear all;
close all

Datos

rt=6400;            %Longitud del Radio de la Tierra
g0=9.8;             %acel. gravedad
rper=rt+600;        %Perigeo
betha=(rper/rt);    %Periogeo Adimensional
ii=-63.4*pi/180;    %angulo orbita -ecuador
lon_ap=-100;        %longitud del apogeo
lon_ap=lon_ap+90;
uv=((rt*g0/1000)^(0.5));  %Adimensional de velocidad
ut=rt/uv;                 %Adimanensional de Tiempo
T=12*3600/ut;             %Período de la órbita adimensional

Cálculos

Se calculan las  constantes de la ecuacion polar de la orbita y se calcula el  semi eje mayor

a=(T/(2*pi))^(2/3);


e=1-betha/a;
k=(a*(1-e^2))^0.5;


t=linspace(0,24,50*24)*3600/ut;     %vector tiempo en seg
E=zeros(1,length(t));               %anomamalía verdadera
error=1e-8;                         %Tolerancia para el cálculo de E
ini=0;

Cálculo de E en función del tiempo

Se resuelve la ecuación

Para ello se utiliza el método de Newton o bisección, para ubicar la solución de E, en la ecuación que la relaciona E con el tiempo

for i=1:length(t)
ini=0;
f=t(i)-T/(2*pi)*(E(i)-e*sin(E(i)));      %Función a buscar solución
z=1;
while abs(f)>error
E(i)=ini;
f=t(i)-T/(2*pi)*(E(i)-e*sin(E(i)));
fa=f;
while sign(f) == sign(fa)   %Búsqueda del cambio de signo
E(i)=E(i)+pi/180*10^(2-z);
fa=f;
f=t(i)-T/(2*pi)*(E(i)-e*sin(E(i)));
ini=E(i)-pi/180*10^(2-z);
end
z=z+1;
end
end

Determinación v y rho para cada tiempo

La anomalía verdadera v, se obtiene a partir de E con:

La ecuación de la Elipse, que da la posición del radio adimensional rho para cada ángulo de la anomalía verdadeda v es:

v=2*atan(((1+e)/(1-e))^0.5*tan(E/2));


r=k^2./(1+e*cos(v));
rm=r.*rt;

Paso a ejes cartesianos

y=rm.*cos(v).*cos(ii);
z=rm.*cos(v).*sin(ii);
x=sqrt(rm.^2-y.^2-z.^2).*(-sign(sin(v)));

Paso a coordenadas polares

tm=atan2(y,x)-pi/2;
fim=atan2(z,sqrt(x.^2+y.^2));

Se toma en cuenta el movimiento de la tierra y la posición del perigeo

tm=tm-2*pi/(24)*t*ut/3600+lon_ap*pi/180;

Orden del los vectores para que queden entre -180 y 180

for i=1:length(tm)
while abs(tm(i))>pi
tm(i)=tm(i)-2*pi*sign(tm(i));
end
end

Ordenamos el Vector para que no aparezcan líneas cruzadas

i=1;
while i < class="keyword">if sign(tm(i))~=sign(tm(i+1)) &&  abs(tm(i))>pi/180*100
tm =[ tm(1:i) NaN  tm(i+1:end)];
fim=[fim(1:i) NaN  fim(i+1:end)];
i=i+2;
end
i=i+1;
end

Zonas de cobertura

Se calcula el angulo de cobertura para cada posición de la órbita y se ordenana los vectores para cuando cruzan de una parte del planeta a la otra

alfa=.57006302;
g=pi/2-alfa-asin(sin(pi/2+alfa)*rt./rm);
for i=1:length(g)
xz(i,:)=tm(i)+g(i).*cos(linspace(0,2*pi,60));
yz(i,:)=fim(i)+g(i).*sin(linspace(0,2*pi,60));
end

Ordenamiento de la zona de cobertura

[l,a]=size(yz);
for i=1:l
for j=1:a
if abs(yz(i,j)) > pi/2
yz(i,j)=pi-yz(i,j);
xz(i,j)=-xz(i,j);
end
end
end

Gráfico de La Órbita

Se carga la costa del mundo, la órbita y la zona de cobretura, además se grafica en forma de animación para ver la zona de cobertura en cada proyección de la orbita del satélite

load coast;                            %Carga costas del mundo
for i=1:length(g)
plot(tm*180/pi,fim*180/pi,'r')     %Grafica Órbita
hold on
plot(long,lat,'k')                 %Carga mapa
xlim([-180 180])                   %Límites
ylim([-90 90])
plot(180/pi*xz(i,:),180/pi*yz(i,:),'.')
plot(tm(i)*180/pi,fim(i)*180/pi,'o')
pause(.01)
hold off
end

Audio - Tacometro con Matlab ( Como medir las RPM con Matlab y un microfono)

2007-07-24T15:53:00.005-04:00

Los sonidos que emite un motor de explosiones- o alternativo- son más agudos al incrementar las RPM -su velocidad de rotación del cigüeñal, aumenta su frecuencia-. Por otro lado, los sonidos son más graves cuando las RPM son más bajas. El Objetivo es usar esta idea para hacer un "Audio-Tacómetro" que funcione con la frecuencia del sonido ¿ la nota musical- que haga el motor, usando la Trasformada de Fourier -FFT- y Matlab. La idea es hacer una especie de afinador para el motor que en vez de mostrar la nota muestre la frecuencia del sonido. Conocida la frecuencia del sonido hay que saber como se relaciona esta con las RPM.

Supuestos

La fuerza de la explosión del motor es de forma cualquiera, pero es periódica, por lo que puede usarse la transformada de Fourier ¿La FFT - para descomponerla en una suma de fuerzas periódicas de la forma.
El periodo de la fuerza, es el periodo de tiempo entre cada explosión producida. La velocidad angular, las RPM, que es lo que queremos medir, dependen del movimiento sube-baja de cada pistón, y este a su vez depende de la periodicidad de las explosiones. Por lo tanto, podemos dejar la velocidad angular en función de la frecuencia de las explosiones, que es lo que vamos a medir.

Midiendo la frecuencia de las explosiones con sonido

Hay varias formas de medir el periodo de las explosiones, por ejemplo con el encendido de las bujías. En realidad lo que se busca es la frecuencia fundamental del sonido que musicalmente, es la nota que hace el motor.
La respuesta del sistema a la fuerza excitadora ¿la Fuerza de la explosión-es el movimiento de forma de vibración. Esta vibración se transmite al aire, y es el sonido, que se mide por micrófono. Se supone el sistema lineal, entonces la respuesta de este es la suma de las respuestas provocadas por cada sumando de la descomposición de la fuerza periódica del motor por la serie de Fourier. Suponiendo que la fuerza es no impulsiva, con periodo mayor al periodo natural de vibración de la estructura que produce el sonido. Esto para suponer que la respuesta tiene una parte transiente y una estacionaria. Esta última es la que nos importa porque contiene la frecuencia de la fuerza excitadora.

En el mejor de los casos, la fuerza excitadora es sólo una sinusoide, que en el espectro se ve una sola línea a la frecuencia de la fuerza, pero esto no pasa, pues hay más componentes múltiplos en el espectro, entonces lo que se que se debería obtener para una RPM constante, debería mostrar las respuestas provocadas por cada sumando de la fuerza excitadora. Como cada respuesta de cada sumando de serie de la fuerza corresponde a un múltiplo de la frecuencia fundamental -por los n de la serie-, la diferencia entre cada línea del espectro corresponde a la frecuencia fundamental. Como la frecuencia fundamental depende de la velocidad angular del cigüeñal si se conoce esta, se puede conocer la velocidad angular del cigüeñal.

Relación de la frecuencia fundamental con la velocidad angular

Para un motor de 4 tiempos de 4 cilindros, tenemos 2 explosiones por cada vuelta o revolución de cigüeñal.
la frecuencia de las explosiones es:
Así las RPM en función de la frecuencia de las explosiones es
Para un motor de dos tiempos y un pistón, como el del avioncito del GIANT, se calcula de la misma manera y la relación da
El régimen de funcionamiento para este motor esta aproximadamente entre las 2000 rpm correspondientes al ralentí, y un poco mas de 9000 que es el máximo.

Esto es 33,3[Hz] y 150[Hz] por lo que las notas musicales que se podrían hacer con el motor - poco más de una octava- se ven en el gráfico.

Las Mediciones

Para medir se usa un micrófono de PC que en las especificaciones dice que tiene una respuesta de frecuencia de 100 Hz a 10 kHz, es decir, que deja fuera las frecuencias que queremos medir.
Curiosamente cerca del motor de mi auto, el micrófono registra señales de bajas frecuencias, ¿Cómo ocurre esto? No sé, sólo sé que se puso cerca del motor del auto y mostraba una señal y su frecuencia correspondía con las RPM. Para el caso del avioncito del GIANT, el micrófono no es capaz de registrar las frecuencias bajas, pero empieza a registrar una señal periódica concordante con las RPM del motor a las 5500-6000 rpm unos 100[Hz] Aproximadamente.

Para medir se hace un programa en Matlab que grabe un tiempo pequeño, le haga la FFT, y muestre el gráfico. La grabación se hace con el comando wavread, donde TM es el tiempo de la grabación y FS es la frecuencia de muestreo.

Y=wavrecord(TM*FS,FS);

Para hacer el FFT hay que tener ciertas consideraciones. El algoritmo que hace el FFT se puede bajar de la página de The MathWorks ¿ Support , este entrega todo listo, yo la hice una funcion que se llama mifft (aunque no es mía) . Usamos como frecuencia de muestro el doble de la máxima frecuencia que puede alcanzar el motor, 600[Hz] (Criterio de Nyquist). Graficamos para cada instante el espectro, para tratar de hacerlo en tiempo real y ver los resultados. A continuación el FFT instantáneo.

clear all
FS=1000; %input(' Frecuencia de Muestreo: ');
TM=.1;
while 1
    Y=wavrecord(TM*FS,FS); 
    [x,y]=mifft(Y,FS);
    hold off   
    plot(0,0)
    [yo,i]=max(y);
    xo=x(i);
    frecf=num2str(xo);  
    frecf=['La fundamental :',frecf];
    text(3*FS/8,3/4*max(y),frecf)
    text(3*FS/8,2/4*max(y),frecf)
    hold on
    
    stem(x ,y ) 
    pause(TM)
    clc
    %fprintf(1,'RPM         %.0f \n',30*xo)
    %fprintf(1,'frecuencia  %.0f  ',xo)
    if max(Y) >=.99 | min(Y)<=-.99
        fprintf(1,'\n\n   ¡ONDA SATUDADA!',xo)
    end 
end

Un problema que puede ocurrir es que la señal medida puede estar saturada y el FFT mostrará muchos armónicos que no están presentes. Por eso graficamos para cada instante la señal que graba el computador, de forma parecida a un osciloscopio, esto para asegurar que no existe la saturación, y para ver la forma de onda y asegurarnos de que sea periódica. En este caso ocupamos una frecuencia de muestreo alta para ver la forma de la onda, 44100 [Hz] La idea es tener ambos gráficos. A continuación el "Osciloscopio Hechizo"

clear all
FS=44100/2;
TM=.1;
amx=.01;
while 2<3
    Y=wavrecord(TM*FS,FS); 
    hold off
    X=linspace(0,length(Y)/FS,length(Y));
    plot(0,-amx)
    hold on
    plot(TM,amx)
    plot(X,Y)
    pause(.5*TM)
end

Aquí se ve la diferencia de un La440 Hz con la señal saturada y no saturada
Espectro de de la onda saturada

Vista de la onda saturada

Espectro de onda no saturada
Onda no saturada

La prueba

Poniendo el micrófono en un lugar seguro y no caliente, cerca del motor de un Suzuki Aerio® para las 1500 rpm medidas en el tacómetro del vehículo que tiene una sensibilidad de 500 rpm (se mantuvo lo más cerca de la rayita de los 1500) y asegurando que no esta saturada la señal.

Aca esta la grabación del motor del Suzuki Aerio® para las 1500 rpm.

El fft de la grabación es:

La línea mde mayor amplitud esta cercana a los 50[Hz]

La onda muestreada a 44100[Hz]es esta, notar que no está saturada.

La onda muestreada a 600[Hz]es esta.

El espectro es como en el mejor caso, es decir, solo muestra una línea de frecuencia, que corresponde a los 50[Hz] que son las 1500[rpm]=30*50[Hz]. Esto pasa porque se ocupa una frecuencia de muestreo baja, y no se muestran los armónicos mayores a 300[Hz], Los armónicos menores a 300 [Hz] existen, pero son muy pequeños. La onda muestreada a 44100[Hz] tiene forma periódica pero no de sinusoide, por esto debe tener armónicos.

El problema de estas mediciones ocurre cuando se acelera de forma rápida, pues la muestra que se graba y se analiza no es una sinusoide de periodo constante, como lo exige el FFT. El periodo de la señal cambia, y el programa entrega una frecuencia que no corresponde necesariamente al sonido del motor.

Caso del avión del GIANT

El sistema funciona con el avioncito midiendo RPM mas altas que 5000~6000 RPM. Hay dos problemas para medir las RPM con este método para el avioncito...

1) La saturación de la señal medida por la gran amplitud del sonido del motor a altas RPM, esto se soluciona cubriendo el micrófono con géneros o esponjas aislantes, pero ocurre que cuando esta funcionando al ralentí el sonido es muy bajo y el micrófono no lo registra de buena forma.

2) Los sonidos producidos por el motor cuando este funciona a 5000-6000 RPM son cercanos a 100 Hz, limite inferior del rango de frecuencias que capta el motor. A diferencia de lo que ocurre en el auto, el micrófono no es capaz de captar estas ondas de baja frecuencia.

Cualquier ayuda, sugerencia comentario es bienvenida

Francisco Roco

Musica en Matlab (Como hacer musica en matlab)

2007-05-12T18:47:00.007-04:00

Imagínese usted la versión de famosas canciones interpretadas por un Nintendo.
Eso es lo que he intentado hacer con Matlab. El resultado a continuación:

Escuche:

Versión de la cortina del Tetris, se llama Korobushka, Tetris

Versión de la introducción de la "Reina del Tamarugal", Tirana

Versión de la cortina de "Informe Especial"

Versión de la canción de los Kjarkas "Llorando se fue" más conocida como "Lambada" (La versión del plagio)

Este es una humilde propuesta, pues hay maneras mas eficientes de hacer música con otros programas, como Noteworthy Composer o el Guitar Pro , etc. O simplemente su merced puede tocar música con instrumento. Recomiendo no tocar harmónica porque pone nerviosa a la gente, sobre todo a tus compañeros de Universidad cuando hay certamenes.

La ventaja de la propuesta de hacer música en MATLAB es que acá empiezas de cero, y puedes crear tus propios tipos de sonidos, además de entender como funciona la música. Antes he escrito un poco sobre estos temas. Más abajo se hace referencia al sonido y a sus características, además de el cálculo de una frecuencia para una nota en particular.

Para explicar mejor el programa veamos las partes.
El programa lee música escrita en formato parecido al de los celulares NOKIA, arma las señales de sonido , las reproduce, y las guarda.

Para ello se ocupan 3 funciones, aca se incluye el vínculo a las sintaxis:

a) música , Programa principal que una vez que tiene las señales las suma las reproduce y las guarda como .wav

b) notas Función que con las notas ingresadas en formato de escritura parecido al del celular NOKIA, crea señales en un vector con una frecuencia de muestreo dada.

c) perc Crea la percusión variando el tiempo de sonido de un sonido parecido a un ruido de interferencia electríca, o ruido blanco.

Los detalles y observaciones importantes se discuten a continuación:

1) Ingreso de notas, lectura de notas, creación de la señal
Se introducen como en los celulares NOKIA de la forma

[Duración][Nota en clave americana][Octava del piano]

-La duración de la redonda será 1, y las figuras siguientes serán nombradas por el inverso del múltiplo que le corresponden. Así:
redonda 1, blanca 2, negra 3, corchea 8, semicorchea 16.
- La nota en clave americana A hasta la G, si es sostenido sera acompañada por una s, los bemoles se hacen con los sostenidos.
-Las octavas son 5 , la central del piano es la 3.

Ejemplo: La 440, corchea de la octava central del piano ..................4a3
Ejemplo: Do# blanca de la "cuarta octava" (la del co central se toma como la tercera) ..................2cs4

Hice una misma función que leyera la nota y guardara en una matriz, la duración y la nota
asi [duración , notas], si tenemos n notas con sus respectivas duraciones será una matriz de n x 2 filas.

La creación de la señal se hace con la matriz anterior, se lee la frecuencia y la nota, y luego se arma una señal con la frecuencia y duración dada.

Para armar el sonido se hace la sinusoide y se le suman los armónicos que queremos que tengan. Yo personalmente intentaba imitar el sonido de un nintendo, asi que le puse unos armonicos que encontre por ahí en internet, la señal me queda como una suma de esta manera:

sin(fr*2*pi.*x)+1/4*sin(3*fr*2*pi.*x)+1/6*sin(5*fr*2*pi.*x)+1/8*sin(7*fr*2*pi.*x)

Para armar el trozo de sonido para el par frecuencia/duración dejo un espacio libre con silencio , para que no suenen pegadas:

for i=1:A
fr=N1(i,1);
t=N1(i,2);
x=(0:(1/fm):t*8/10);
y1=[y1  sin(fr*2*pi.*x)+1/4*sin(3*fr*2*pi.*x)+1/6*sin(5*fr*2*pi.*x)+1/8*sin(7*fr*2*pi.*x) zeros(size((0:(1/fm):t*2/10)))];
end

Las notas/frecuencias las armamos asi


% ==FRECUENCIAS==
r=2^(1/12);
fr=440;%frecuencia del LA
fr=r^-9*fr;%el Do
c=r^0*fr*[1/4 1/2 1 2 4];
cs=r^1*fr*[1/4 1/2 1 2 4];
d=r^2*fr*[1/4 1/2 1 2 4];
ds=r^3*fr*[1/4 1/2 1 2 4];
e=r^4*fr*[1/4 1/2 1 2 4];
f=r^5*fr*[1/4 1/2 1 2 4];
fs=r^6*fr*[1/4 1/2 1 2 4];
g=r^7*fr*[1/4 1/2 1 2 4];
gs=r^8*fr*[1/4 1/2 1 2 4];
a=r^9*fr*[1/4 1/2 1 2 4];
as=r^10*fr*[1/4 1/2 1 2 4];
b=r^11*fr*[1/4 1/2 1 2 4];
s=0;
i=1;
ii=1;

2)LA PERCUSIÓN
Se hace de forma similar a la anterior pero solo se ingresa la duración. 1 2 5 8 16 . La señal que se crea es algo parecido ala interferencia , que suena como ch chc chhhhhchchchchch (asi).
Es absolutamente igal , lo que varía es como armamos la señal.

y1=[y1,o*sin(100000*2*pi.*x.*100000000000000.*sin(x)) zeros(size((0:(1/fm):t*6/10)))];

ese seno es la interferencia que suena cchhchchchchhchcchchc.

3)LOS INSTRUMENTOS
Ahora que se tienen las señales hay que hacerlas sonar.
Como ya creamos las señales de sonido, en el programa hay 3 de melodía más la percusión, lo que hay que hacer , es sumarlas

y=1/4*(y1+y2+y3+y4);

y reproducirlas, vamos, que suene...

sound(y,fm)

Como tenemos almacenada la señal en el vector y y tenemos su frecuencia de muestreo, con el comando

wavwrite(y,fm,8,Nom)

guardamos en un archivo formato wav la canción.
Al inicio de la sintaxis están los códigos de algunas canciones. Están los 3 instrumentos más la percusión.
Cualquier comentario duda o sugerencia, sólo hazlo :D

Francisco Roco

Como hacer un gif con Matlab

2007-03-13T19:01:00.004-04:00

Para algunos trabajos o presentaciones, mejor que un esquema , es una presentación. En nuestro caso teniendo ecuaciones de movimiento podemos hacer una en formato gif , que se pueden pegar en un Pagüer-point.

Más abajo se ve como hacer animaciones en matlab, con un sistema mecánico. (Aquí)

Esas animaciones las podemos guardar como videos *.avi, pero como ese formato pesa mucho, personalmente prefiero los gifs.

La gran pregunta es ¿Cómo hago los gifs?

Buscando en internet encontré la página de un gringo que dice como hacerlo, y en realidad voy a decir como él lo hace, hay que agradecerle a él. La página es esta: http://uwcem.ece.wisc.edu/creatingAnimatedGifs.htm
Del poco inglés que sé pude rescatar lo siguiente: (Corrígeme si estoy mal)

1) Necesitamos Matlab y un programita que esta en este vínculo
http://uwcem.ece.wisc.edu/downloadable/GaLE10.exe

2) Para cada gráfico hecho, que se supone que será un cuadro en tu animación, debe estar en un ciclo "for" como muestra el ejemplo , despues deberás agregar unas lineas de sintaxis que crearán imágenes bmp.

3) se crearan tantas imágenes como tu ciclo for, esas imágenes , las metes al programa para los gifs y listo tendrás tus gif hechos en matlab

Ejemplo

El ejemplo son puntos de una párabola puse las ecuaciones de movimiento parametrizadas...

t=-2:.1:2;
x=100*t;
y=-9.8/2*t.^2+20;
for i=1:length(t)
plot(-200,0)
hold on
plot(200,25)
plot(x(i),y(i),'o')
plot(x(1:i),y(1:i),'.')
pause(.05)

% ===para Hacer GIF===
Image = getframe;
P = frame2im(Image);
number = num2str(i);
extension = '.bmp';
filename = [number,extension];
imwrite(P,eval('filename'), 'bmp');
hold off
end

Lo que se debería ver como el gif de arriba
Cualquier duda comentario aporte crítica es bienvenida

Francisco Roco

Como calcular la frecuencia de una nota musical y hacer música con Matlab / Gnu-Octave

2007-01-23T11:23:00.008-03:00

Bueno como más abajo dice como hacer un sonido, ahora podemos hacer música, pues la música no es más que un conjunto ordenado de sonidos. Veamos como ordenarlos...

Si disponemos de un lenguaje de programación, que pueda actuar sobre algún dispositivo de sonido, o de cualquier cosa que pueda vibrar o producir sonidos a una frecuencia audible, sólo falta saber a que nota musical corresponde cierta frecuencia en la variación de presión del aire.

Los sonidos que se usan son los que corresponden a las notas musicales. A cada nota musical le corresponde una frecuencia. (Ver Características del Sonido)

Introducción al Tema

(Si quieres ver la fórmula para el cálculo sáltate esta parte)

Muy grosso modo la música posee, entre otras características importantes:

Altura, las distintas notas que se representan en los espacios y líneas en el pentagrama

Duración que se representa con las figuras, redonda, blanca, negra, corchea , etc.

Timbre el tipo de sonido si es de un piano, un motor sonando o de una soprano cantando

Para hacer música, del tipo occidental, (Otras culturas ocupan distintas escalas) disponemos de un número discreto de sonidos limitados por el rango audible de un ser humano, por las frecuencias 20[Hz] como cota inferior y por arriba de 20.000[Hz], lo que varía según la persona. Pues para algunas personas da lo mismo escuchar Reguetón por los altavoces de un teléfono móvil o escuchar una grabación del deutsche grammophon en formato alta fidelidad, y para otras no da lo mismo, pero eso es para otro tema.

Los sonidos usados por los músicos, las notas musicales, do re mi fa sol la si y sus respectivos bemoles y sostenidos, tienen frecuencias conocidas que se pueden calcular, y se pueden hacer en algún lenguaje de programación -como C, C++, Java, Matlab, Pascal, Basic, etc - para que suenen. La historia de las distintas escalas utilizadas -tipos de frecuencias y sus cálculos para las notas- es un poco larga y tiene que ver con criterios matemáticos (Véase acá), es muy interesante.

Brevemente, la historia empieza con Pitágoras, que inventó la escala que lleva su nombre, pero la parte que ahora interesa ver habla sobre el periodo barroco y Johann Sebastian Bach. Cuando se usaba la escala Pitagórica y uno quería cambiar de tono las canciones (frecuentemente pasa cuando a un cantante no le alcanza el registro), se producía una no correspondencia, es decir, no sonaba correctamente ( sonaba mal o en palabras simples, sonaba feo).

Y ahí es donde aparece la escala cromática, que el famoso Bach popularizó con su obra "El clave bien temperado" o "Das wohltemperierte Klavier" , el nombre lo dice, el instrumento clave estaba bien temperado, es decir, afinado con la escala cromática. ¿Cómo se populariza un estilo de afinación? Componiendo un preludio y una fuga en todos los tonos posibles, es decir un preludio y fuga en Do mayor y menor, en Do sostenido mayor y menor en Re mayor ... y así. Hasta hace poco había leído que Bach había inventado esta afinación, lo que no es cierto, pero esto no quita mérito a semejante obra.

"El clave bien temperado" se puede ver y escuchar en esta página, ahí explican hasta las partes. Es bastante interesante verlo y sobretodo escucharlo.

¿Como se construye la escala cromática?
La respuesta es "La distancia entre semitonos debe ser la misma". Eso desde el punto de vista musical. Porque esas palabras traducidas a la práctica son así :

Para armar la escala temperada necesitamos una frecuencia de una nota conocida, llamada de referencia, en nuestro caso y por norma ISO-16 el LA de la octava central del piano tiene 440 [Hz]. (Para afinar los músicos usan los diapasones). En tiempos de Bach cada cual elegía su frecuencia de referencia, lo que no causaba problemas si se toca solo un instrumento.

¿Cual es la relación entre frecuencias-notas en la escala cromática?

Fórmula para calcular la frecuencia de una nota.

Dada una frecuencia conocida, la que llamaremos f_Ref, la nota distante a i semitonos sobre la frecuencia de referencia f_Ref tendrá la frecuencia:

La norma ISO 16:1975 dice que la frecuencia de referencia para el La de la octava central del piano es:

La norma agrega también que debe ser generado con una precisión de 0.5 [Hz].

Entonces, siguiendo la norma, la frecuencia del i-ésimo semitono a partir del La de la octava central del piano es:

Ejemplo

Calculemos la frecuencia del Do que sigue del La-440, el de la octava a la derecha de la octava central:

Veamos a que semitono corresponde. Para ello vemos en la figura y contamos:
- 0 La
- 1 La#
- 2 Si
- 3 Do
El Do corresponde al semitono 3 a partir del La.

Reemplazamos en la fórmula
La frecuencia del do es:

Hacemos música en Matlab

Como se tienen las notas, basta con darles duración para crear música, usando como guía la parte vista más abajo, Como hacer sonido, escribimos la siguiente sintaxis, creando primero con la forma vista arriba, las notas musicales, y luego creando sonidos.

%==Frecuencia de Muestreo==
fm=8120;

% ==Elaboracion de las notas musicales==
r=2^(1/12);

% ==frecuencia de referencia LA==
fr=440;

% Tomamos la frecuencia del do para definir las octavas de mejor manera
fr=r^-9*fr;
c=r^0*fr*[1/4 1/2 1 2 4];
cs=r^1*fr*[1/4 1/2 1 2 4];
d=r^2*fr*[1/4 1/2 1 2 4];
ds=r^3*fr*[1/4 1/2 1 2 4];
e=r^4*fr*[1/4 1/2 1 2 4];
f=r^5*fr*[1/4 1/2 1 2 4];
fs=r^6*fr*[1/4 1/2 1 2 4];
g=r^7*fr*[1/4 1/2 1 2 4];
gs=r^8*fr*[1/4 1/2 1 2 4];
a=r^9*fr*[1/4 1/2 1 2 4];
as=r^10*fr*[1/4 1/2 1 2 4];
b=r^11*fr*[1/4 1/2 1 2 4];

%El Silencio, su frecuencia es cero
s=0;

%==Parte para Duracion==
% Definimos la duracion en segundos de la negra.
% tomada como 1 tiempo
n=.5;

%Hacemos un vector N = [ nota , duracion ]
NOTAS=[g(3),e(3),s,e(3),f(3),g(3),e(4),e(4),c(4)];
DURACION=[n/2,n/2,3*n/2,n/2,n/2,n/2,n,n,2*n];
N=[NOTAS',DURACION'];
y=[];

%==Armamos la señal==
for i=1:length(N)
fr=N(i,1);
t=N(i,2);
x=(0:(1/fm):t);
y=[y sin(fr*2*pi.*x)];
end

% Hacemos el sonido
% en este comando
% esta toda la magia

sound(y)

Puse en NOTAS y DURACION una canción muy conocida... jajaja, buenísima.

El método para ingresar notas se puede mejorar y también se pueden añadir armónicos a la onda. Por ejemplo to hice un método para que introducir 3 instrumentos, percusión a través de notas parecidas al celular NOKIA
El programa está acá

Francisco Roco

Características del Sonido

2007-01-23T10:48:00.004-03:00

Un sonido tiene características que se notan en la forma de onda.
Estas son:

Altura

Intensidad

Timbre

La Altura del Sonido

Al hablar de altura se hace referencia a la nota musical que se hace, en otras palabras, que tan agudo o grave es el sonido. Mientras más grave es el sonido(Como el de un Contrabajo o un cantante bajo), menor será la frecuencia de mi onda y mayor su longitud de onda. Por el contrario si el sonido es agudo (Como un flautín o una soprano), la frecuencia de la onda será mayor y menor su longitud de onda.

Si numeramos las teclas del piano de izquierda a derecha, es decir, de grave a agudo (Supongamos un piano de 96 teclas, [tienen al rededor de 88 algunos ]), y a cada tecla numerada hacemos corresponder una nota y por ende a una frecuencia, lo que veriamos es eso:

Intensidad

La intensidad del sonido hace referencia a lo que en los equipos de música se le llama "volume", es la cantidad potencia que tiene el sonido, por eso que los parlantes traen calcomanías que dicen 4000[W], por ejemplo.
En la onda la intensidad se ve como la amplitud de la onda , más fuerte suena, más intenso , más amplitud tiene el sonido, por el contrario más despaco , menor su amplitud.
La unidad de amplitud frecuentemente usada es el Decibelio [dB].

El Timbre

El timbre de un sonido dice relación con la fuente que produce el sonido que escuchamos. Me explico Puedo tener dos situaciones un cantante haciendo un La, y un pianista haciendo un la. En ambos casos la altura es la misma pues hacen la misma nota, la misma frecuencia de 440[Hz]. Imaginemos que la amplitud es la misma, así tienen el mismo volume, o intensidad. Pero aún así podemos diferenciar el sonido del piano con el del cantante, por la característica llamada timbre.

En la onda esto corresponde a la cantidad de armónicos que tiene la onda.
Los armónicos de una frecuencia X, son ondas de amplitud más pequeña que se suman a la frecuencia X , que se llama frecuencia principal.

Por ejemplo, abajo vemos la frecuencia principal en color negro, su frecuencia armónica tiene una amplitud de 1/6 y una frecuencia tres veces mayor que la principal, la suma de estas dos onditas es la onda azul.

En el caso del sonido es un poco mas complicado, pues los armónicos son mucho más que dos , como en el caso anterior.

Los sintetizadores , o los teclados electrónicos , crean electrónicamente con matlab distintos armónicos para crear sonidos como trompetas, pianos , voces etc. lo difícil es ver que armónicos agregar.

A continuación pongo ejemplos de algunos "LA-440[Hz]" que he grabado y he graficado.
Nótese los armónicos y que la onda tiene ciclos, no es difícil encontrarlos:

La por Roco

Nótese la onda graficada...

La por un Diapasón (Este idealmente no tiene armónicos)

La por un piano

El programa para Matlab con los archivos de sonido de distintos "la" esta Aquí.

Chao

Francisco Roco

El Sonido y matlab

2007-01-14T03:21:00.003-03:00

¿Cómo hago sonido con el Matlab?

El computador funciona de forma discreta, no en forma continua, es decir, trabaja con números, que guardamos en vectores.
Un ejemplo de vector es [1,2,3,4,5,6,10] (Un saludo pa’ todos los locos)

Para hacer sonidos necesitamos una señal que en este caso es un vector con muchos números, de repente millones. Este vector tiene que contener el sonido.

Un sonido es una vibración, una onda, que la vemos de forma sinuidal o de forma parecida (En realidad es consecuencia de ello pero para formar una idea buena lo planteo así).

Lo que se guarda en el vector son valores de la señal. Una analogía para entender esto la visualizo así: Una cuerda que se mueve, vibrando (Fig. 1) y le sacamos una foto en un instante. Cuando la cuerda esta quieta, no vibrando, decimos que está en el valor cero (Fig. 2). Mientras empieza a moverse de izquierda a derecha va adquiriendo valores, los que medimos al medio de la cuerda desde el punto donde esta quieta -cuando esta en cero- (Fig.3), su elongación máxima será 1 a la derecha, y -1 a la izquierda (Fig.5 y 6).

Bueno, con cada foto que le saque a mi onda que representé como cuerda, tendré un valor que guardaré en el vector. Pero como el movimiento es una vibración, tengo que guardar estos valores cada cierto intervalo de tiempo, algo así como sacar fotos cada 1 segundo a mi onda, y es ahí donde aparece el concepto de “Frecuencia de Muestro”.

La frecuencia de Muestreo fm. es, la cantidad de muestras (fotos en mi analogía), que se toman en 1 segundo. Si tomo la grafica de los valores que va tomando mi onda (cuerda), en el “eje y” pongo los valores en función del tiempo, y en el “eje x” el tiempo la frecuencia de muestreo en la figura es 1[Hz], pues muestreo (tomo una foto), 1 vez por segundo.
Lo que guardo es lo que se ve arriba

Frecuencia de muestreo

Lógicamente cuando tengo una mayor frecuencia de muestreo mejor será la calidad del sonido, porque estaré representando con más puntos mi función.
Por ejemplo el teléfono funciona a 8000 [Hz], y se escucha más o menos no más, (Sobre todo cuando hacen esos contactos con las unidades móviles en la Radio); Los CD tienen una frecuencia de muestreo de 44.100[Hz] la radio de 22100[Hz], bueno como dato no más.

Hago lo comentado arriba en matlab con el comando sound que reproduce un vector con su frecuencia de muestreo

%Frecuencia de muestreo
fm=44000;
%frecuencia sonido, nota musical La
fr=440;
%Duracion en segundos
duracion=1;
t=0:1/fm:duracion;
% Señal senoidal
y=sin(fr*2*pi*t);
%La hacemos sonar
sound(y,fm)
%La guardamos como wav
wavwrite(.5*y,fm,'la440.wav')
%Graficamos un pedacito para ver la forma
t=t(round(68.5*length(t)/100):round(70*length(t)/100));
y=y(round(68.5*length(y)/100):round(70*length(y)/100));
plot(t,y)

Para los que no tengan Matlab:
1) Sonido que se debería escuchar
2) Gráfico de un pedacito de señal, nótese que los valores que da están entre 1 y -1, y que en el eje x esta graficado el tiempo en segundos. Se nota la forma de sinusoide

Su servidor...

Francisco Roco

Animación del Motor Radial en Matlab

2007-01-01T17:45:00.003-03:00

Se muestra a continuación una animación de un Motor Radial en la página se verá una animación para que quede clarísimo como es uno, la que se expone aquí es un poco más simple, pero no menos didáctica.

Esta se ve así(En versión pequeña):

La animación es simplemente una sucesión de imágenes que van cambiando, en la cual cada una permanece un período de tiempo muy pequeño.

El cálculo de la dinámica, más abajo, se plantea para varias posiciones, en este caso las mayas tiene como variable de entrada t que es el ángulo de la manivela (cigüeñal).
Luego se calcula la posición de los pasadores de la biela principal.
Posteriormente la posición de los pistones. Para ello se necesitan como entradas un t inicial y los valores tentativos respectivos.
Como son cinco pistones se cálcula la posición de cada uno de ellos en una función independiente, la idea era hacer una para todos, pero se me complicó el asunto en sun momento y al final hice 5 :D.

Para hacer todo esto hay dos partes

A) Programa principal
B) Funciones de Cálculo.

Sintaxis del Programa Principal

En la función principal aquí abajo se ve: La parte de los valores iniciales. Un ciclo For en el que se calcula para cada t la posición de todos los pistones.Obviamente va variando el t, en esta parte se llaman a las funciones que caculan toda la imagen. Y por último la parte que hace el gráfico, que son puros plot más la parte que usa el comando pause.

Notar las tres partes indicadas:


%Valores iniciales
clc
clear all
hold off

%numero de cálculos
nc=2%5*360/5+1;

%tiempo depausa para una imagen
dt=0.01;
paso=pi/4%5*pi/(180);

%geometria
lm=15; %manivela
lb=45; %biela
rb=8; %radio de la biela
e=16; %ancho piston
d=3.4; %ancho biela
np=10; %numero de puntos para el calculo de la curvatura de las biela
 
%Valor inicial
t=0;
X=[];
Y=[];
x1=lm+lb;

% valores estimativos iniciales
x1=40;
tb=80*pi/180;

%piston2
t2=pi/8;
a2=(90-72)*pi/180;
x2=50;

% Piston3
t3=45*pi/180;
a3=54*pi/180;
x3=40;

%Piston4
t4=pi/5;
a4=a3;
x4=10;

% Piston5
t5=-30*pi/180;
a5=18*pi/180;
x5=40;


%  Calculos

for i=1:nc
 [tb,x1] = pist1(t,tb,x1,lm,lb,rb);
 x=[0,lm*sin(t)];
 y=[0,lm*cos(t)];
 X0(i,:)=x;
 Y0(i,:)=y;
 [cx,cy]=cuad(0,x1,e,0);
 CX1(i,:)=cx;
 CY1(i,:)=cy;

 % Biela
 pm=lm*[sin(t),cos(t)];
 [ppx,ppy]=pen(pm,tb,rb);
 Px(i,:)=ppx';
 Py(i,:)=ppy';
 tt(i)=t;
 xx2(i)=x2;
 
 % piston1
 X1(i,:)=[ppx(1),0];
 Y1(i,:)=[ppy(1),x1];
 
 % piston2
 [t2,x2] = pist(ppx(2),ppy(2),t2,x2,lm,lb,rb);
 [cx,cy]=cuad(x2*cos(a2),x2*sin(a2),e,(90-72)*pi/180);
 CX2(i,:)=cx;
 CY2(i,:)=cy;
 X2(i,:)=[ppx(2),x2*cos(a2)];
 Y2(i,:)=[ppy(2),x2*sin(a2)];
 xx2(i)=x2;
  
 % piston3
 [t3,x3] = pist3(ppx(3),ppy(3),t3,x3,a3,lm,lb,rb);
 [cx,cy]=cuad(x3*cos(a3),-x3*sin(a3),e,pi-a3);
 CX3(i,:)=cx;
 CY3(i,:)=cy;
 X3(i,:)=[ppx(3),x3*cos(a3)];
 Y3(i,:)=[ppy(3),-x3*sin(a3)];
 xx3(i)=x3;

 % piston4
 [t4,x4] = pist4(ppx(4),ppy(4),t4,x4,a4,lm,lb,rb);
 [cx,cy]=cuad(-x4*cos(a4),-x4*sin(a4),e,a4);
 CX4(i,:)=cx;
 CY4(i,:)=cy;
 X4(i,:)=[ppx(4),-x4*cos(a4)];
 Y4(i,:)=[ppy(4),-x4*sin(a4)];
 xx4(i)=x4;

 % piston5
 [t5,x5] = pist5(ppx(5),ppy(5),t5,x5,a5,lm,lb,rb);
 [cx,cy]=cuad(-x5*cos(a5),x5*sin(a5),e,pi-a5);
 CX5(i,:)=cx;
 CY5(i,:)=cy;
 X5(i,:)=[ppx(5),-x5*cos(a5)];
 Y5(i,:)=[ppy(5),x5*sin(a5)];
 xx5(i)=x5;

 t=t+paso; 

 % cigueñal
 [bx0,by0] = biela(X0(i,:),Y0(i,:),d,np);
 
 % las bielas
 [bx2,by2] = biela(X2(i,:),Y2(i,:),d,np);
 [bx3,by3] = biela(X3(i,:),Y3(i,:),d,np);
 [bx4,by4] = biela(X4(i,:),Y4(i,:),d,np);
 [bx5,by5] = biela(X5(i,:),Y5(i,:),d,np);

 BX0(i,:)=bx0';
 BY0(i,:)=by0';

 BX2(i,:)=bx2';
 BY2(i,:)=by2';

 BX3(i,:)=bx3';
 BY3(i,:)=by3';

 BX4(i,:)=bx4';
 BY4(i,:)=by4';

 BX5(i,:)=bx5';
 BY5(i,:)=by5';

 % biela Maestra
 [ppx,ppy]=pbiela(pm,tb,[ppx(1),0],[ppy(1),x1],3,10,rb);
 BPX(i,:)=ppx';
 BPY(i,:)=ppy';
end

BX=[];
BY=[];

a=[90,18,-54,-126,-198,90]*pi/180;

for i=1:6
  P = bloc(a(i),lm,lb,rb,e);
  BX=[BX P(:,1)'];
  BY=[BY P(:,2)'];
end


%  Animación

for i=1:nc
 plot(-100,-100)
 hold on
 plot(BX,BY,'k')%Bloc
 plot(100,100)
 plot(CX1(i,:),CY1(i,:))%CUADRADO 1
 plot(CX1(i,:),CY1(i,:))%CUADRADO 1
 plot(CX2(i,:),CY2(i,:))%CUADRADO 2
 plot(CX3(i,:),CY3(i,:))%CUADRADO 3
 plot(CX4(i,:),CY4(i,:))%CUADRADO 4
 plot(CX5(i,:),CY5(i,:))%CUADRADO 5

 plot(Px(i,:),Py(i,:),'o')%PENTAGONO

 plot(X0(i,:),Y0(i,:),'o')%FIGURA 0 manivela 2 puntos
 plot(X1(i,:),Y1(i,:),'o')%FIGURA 1 pasador biela princ piston 2 puntos
 plot(X2(i,:),Y2(i,:),'o')%FIGURA 2 pasador biela princ piston 2 puntos
 plot(X3(i,:),Y3(i,:),'o')%FIGURA 3 pasador biela princ piston 2 puntos
 plot(X4(i,:),Y4(i,:),'o')%FIGURA 4 pasador biela princ piston 2 puntos
 plot(X5(i,:),Y5(i,:),'o')%FIGURA 5 pasador biela princ piston 2 puntos

 plot(BX0(i,:),BY0(i,:),'m')%Manivela

 plot(BPX(i,:),BPY(i,:),'r')%biela maestra
 plot(BX2(i,:),BY2(i,:))%Biela
 plot(BX3(i,:),BY3(i,:))%Biela
 plot(BX4(i,:),BY4(i,:))%Biela
 plot(BX5(i,:),BY5(i,:))%Biela

 hold off
 pause(dt) %Pausa dt segundos

end

Nótese las tres partes del programa, en la primera se declaran los valores iniciales para las funciones, y junto con ello se determinan para cuantos instantes se realizarán los cálulos, es decir, para cuantas posiciones de la biela principal se hacen los cálculos.
En la segunda para cada instante, es decir, para cada posición de la biela maestra o principal, se calculan las posiciones de los pistones. Esto se hce llamando a las funciones correspondientes en cada caso, se examinará mas abajo este detalle.
En la tercera parte se grafican, para cada instante las posiciones.

Funciones de Cálculo

La Primera función es "pist.m" , para calcular la posicion de los pistones, hay una función por pistón, la que se expone más abajo es la del piston 1 el que se mueve de arriba abajo, la punta de la "estrella". La idea es tener como entrada a la función el angulo de la biela principal t, y con la ecuacion de la maya sacar los valores que faltan.
Nótese que el método de calculo usado es el mismo que se describe en una entrada anterior del blog, si no se conoce ver en Animaciones con matlab, en el que se utiliza una iteración.

function     [tb,x1] = pist1(t,tb,x1,lm,lb,rb)
e=10^-10;
c=1;
f1=lm*sin(t)-(lb+rb)*cos(tb);
f2=lm*cos(t)+(lb+rb)*sin(tb)-x1;

%  Iteracion
while abs(f1)>e | abs(f2)>e
 c=c+1;
 df1x=0;
 df1tb=(lb+rb)*sin(tb);
 df2x=-1;
 df2tb=(lb+rb)*cos(tb);
 A=[df1x,df1tb;df2x,df2tb];
 B=[-f1;-f2];
 X=A\B;
 if X==[0;0]
     break
 end
 tb=tb+X(2);
 x1=x1+X(1);
 f1=lm*sin(t)-(lb+rb)*cos(tb);
 f2=lm*cos(t)+(lb+rb)*sin(tb)-x1;
end
x1=lm*cos(t)+(lb+rb)*sin(tb);

Función para graficar los pistones

Para hacer que la animación sea más realista, aunque esquemática, se usa una función que dibuja un cuadrado orientado en el angulo del bloc del motor, de manera que parezca un piston. La función cuad.m , da un cuadrado, tiene como entrada un punto, con sus dos coordenadas x e y, el lado del cuadrado que queremos e y el ángulo de inclinación de este a. La función tiene como salidalos puntos de un cuadrado (representa un pistón), listo para ser graficado.

function [X,Y]=cuad(x,y,e,a);
r=e*sqrt(2);
b=pi/4-a;
X=[x-r*sin(b),x+r*cos(b),x+r*sin(b),x-r*cos(b),x-r*sin(b)];
Y=[y-r*cos(b),y-r*sin(b),y+r*cos(b),y+r*sin(b),y-r*cos(b)];

Función para Pasadores

La función pen.m calcula la posición de los 5 pasadores de la biela principal, como entrada tiene la posicion de la manivela, con elvector pm, el ángulo tb, y el radio de la biela rb.

function [ppx,ppy]=pen(pm,tb,rb)
f=0;
for i=1:6
 pp(i,:)=pm+rb*[-cos(tb+f),sin(tb+f)];
 f=f+72*pi/180;
end

ppx=pp(:,1);
ppy=pp(:,2);
ppx(1)=0;
ppy(1)=0;
ppx(6)=0;
ppy(6)=0;

Función para dibujar una biela

La función biela.m se usa para dibujar la biela en el gráfico. Como tenemos los puntos de los pasadores de la biela, que son el pistón y la biela principal, se dibuja al rededor de estos un rectángulo con los bordes redondeados. Como entrada a la función, se tienen : los puntos de los pasadores, enlos vectores X1 e Y1, el ancho de la biela d, y el número de puntos n que dibuja la semicircunferencia de la biela, nótese en el ciclo for.

function     [X,Y] = biela(X1,Y1,d,n)
p1=[X1(1,1),Y1(1,1)];
p2=[X1(1,2),Y1(1,2)];
a=atan2(Y1(1,2)-Y1(1,1),X1(1,2)-X1(1,1));
b=pi/2-a;

for i=1:n
 p(i,:)=p2+d*[cos(-b+i*pi/n),sin(-b+i*pi/n)];
end

for i=1:n
 p(n+i,:)=p1+d*[cos(a+pi/2+i*pi/n),sin(a+pi/2+i*pi/n)];
end
p(2*n+1,:)=p(1,:);
X=p(:,1);
Y=p(:,2);

Dibujando la biela Principal o Maestra

La función pbiela,hace una figura con forma de biela principal. Tomando pm la posición de la manivela,el angulo de esta, y

los mismos argumentos que la función anterior.

function [ppx,ppy]=pbiela(pm,tb,X1,Y1,d,n,rb)
f=0;
for i=1:6
 pp(i,:)=pm+3/2*rb*[-cos(tb+f),sin(tb+f)];
 f=f+72*pi/180;
end

ppx=pp(:,1);
ppy=pp(:,2);

t=atan2(ppy(2)-ppy(5),ppx(2)-ppx(5));
ppx(6)=ppx(1)-d*cos(t);
ppy(6)=ppy(1)-d*sin(t);

ppx(1)=ppx(1)+d*cos(t);
ppy(1)=ppy(1)+d*sin(t);

p1=[X1(1,1),Y1(1,1)];
p2=[X1(1,2),Y1(1,2)];
a=atan2(Y1(1,2)-Y1(1,1),X1(1,2)-X1(1,1));
b=pi/2-a;

for i=1:n
 ppx(i+6)=p2(1,1)-d*cos(b+i*pi/n);
 ppy(i+6)=p2(1,2)+d*sin(b+i*pi/n);
end
ppx(6+n)=ppx(1);
ppy(6+n)=ppy(1);

Dibujando el Bloc del motor

Para dibujar el bloc o carcasa del motor usamos la función bloc.m, ella tiene como entradas el angulo de las estrellas a2, el largo de la manivela lm, el largo de la bielalb, el radio de la biela rb y el ancho del pistóne.

function [P] = bloc(a2,lm,lb,rb,e)
r1=-lm+rb+lb-e;
r2=lm+rb+lb+e;
r2=r2-r1;
p1=r1*[cos(a2),sin(a2)];
p2=r2*[cos(a2),sin(a2)];
i=e*[-sin(a2),cos(a2)];
d=e*[sin(a2),-cos(a2)];
P(1,:)=p1+i;
P(2,:)=p1+i+p2;
P(3,:)=p1+i+p2+2*d;
P(4,:)=p1+i+p2+2*d;
P(5,:)=p1+d;

Los archivos *.m

Están para descargarTodas las funciones más el programa. Se agradecería avisar cualquier error.

Francisco Roco

Animaciones con Matlab/GNU-Octave

2006-12-28T03:39:00.009-03:00

Esto puede resultar de ayuda cuando se quiere calcular la dinámica de un sistema mecánico, es decir su posición en función del tiempo, su velocidad o su aceleración para el cálculo de las fuerzas de trepidación. Además de el cálculo, uno puede "ver "estas magnitudes, en un animación.

Animación en Matlab

Simplemente es ir calculando para cada instante la posición de nuestro "sistema", las guardamos en matrices o vectores, y después graficamos para cada instante. Luego si se desea se puede guardar en un gif y se ve como el pistón de arriba ( para ver más Aquí para hacer un gif con matlab

Conocer la posición de un sistema Pistón Biela Manivela

Se plantea la posición del sistema en un instante cualquiera para establecer las ecuaciones que permitan calcularla.

La geometría conocida del sistema es:

x posición del pistón respecto del eje de giro de la manivela

l_b Longitud de la biela

l_m Longitud de la manivela

t₁ Ángulo entre la línea de recorrido del pistón y la manivela

t₂ Ángulo entre la línea de recorrido del pistón y la biela

Las ecuaciones se determinan sumando en cada componente los desplazamientos a partir de un punto y llegando al mismo punto. En este caso, partiendo del punto de rotación de la manivela, los desplazamientos serían:

Es sistema de ecuaciones queda:

Conocida la geometría, el largo de la biela y la manivela, y dado el ángulo de la línea de recorrido del pistón y la manivela, el sistema de ecuaciones que queda, tiene como incógnita la posición del pistón y el ángulo.

Se puede resolver despejando en la primera ecuación el valor de t₂ y luego reemplazar en la segunda ecuación para encontrar x. Pero para generalizar el método, es conveniente resolver de otra forma el sistema para casos en que hay dos incógnitas, ambas en las dos ecuaciones y como argumento de una función trigonométrica, como en el caso del mecanismo de cuatro barras.

Resolución del sistema de ecuaciones

Para ello se define una función F de dimensión 2, que tiene componentes los miembros del sistema de ecuaciones distintos de cero:

La función se iguala a cero por el segundo miembro del sistema de ecuaciones.
Para encontrar las incógnitas se utiliza la aproximación de la función por los dos primeros términos de la serie de Taylor.

Lo que desarrollado es:

Se impone que el valor de la aproximación de la serie de Taylor sea cero, coincidente con el de la función, entonces queda:

Escrito de otra manera:

Por lo que basta con invertir la matriz para conocer el valor de las incógnitas de la aproximación de la función por la serie de Taylor. Dichos valores son una aproximación a los valores que reemplazados en el sistema lo satisfacen. Por lo que se programa una iteración, de manera de partir con ciertos valores iniciales, evaluar la función, calcular las derivadas parciales, invertir la matriz, obtener valores nuevos, sumarlos a los antiguos, evaluar la función para que sea menor a un error pequeño predefinido y si es así estos valores son considerados los que satisfacen el sistema. Si no es así, se vuelve a iterar. Lo que aparece resumido en el diagrama de flujo.

La función de cálculo

Es diagrama de flujo anterior implementado en código GNU-Octave/Matlab se muestra a continuación. Nótese que como valores de entrada están el ángulo de la manivela, para el cual el sistema tiene una única posición y además se incluyen los valores de x y t₂ iniciales, los que deben ser tentativos.
La función entrega como parametro de salidas las incógnitas x y t₂ y ademas de esos dos vectores, xx e yy, en los que se guardan la posición de los puntos para cada instante en los ejes x e y.


function [xx,yy,x,t2] = blogpiston(t1,x,t2)

% Tolerancia del error
e=0.001;

% Dimensiones del sistema
lm=35;
lb=56;

% definicion de F
f1=lm*sin(t1)-lb*sin(t2);
f2=lm*cos(t1)+lb*cos(t2)-x;

while abs(f1)>e | abs(f2)>e
 df1x=0;
 df1t2=-lb*cos(t2);
 df2x=-1;
 df2t2=-lb*sin(t2);
 A=[df1x,df1t2;df2x,df2t2];
 B=[-f1;-f2];
 
 % Resolucion del sistema
 X=A\B;
 
 % Por si acaso no converge
 if X==[0;0]
  break
 end
 
 %valores nuevos de x y t2
 x=x+X(1);
 t2=t2+X(2);
 
 %definicion de F con los nuevos x y t2
 f1=lm*sin(t1)-lb*sin(t2);
 f2=lm*cos(t1)+lb*cos(t2)-x;
end

% Puntos de la maya para el t1 dado
xx=[0,lm*sin(t1),0];
yy=[0,lm*cos(t1),x];

La animación

Para animar el sistema debemos tener calculada la posición del sistema dados varios instantes, o merjor dicho, para varios angulos de la manivela. Así, si se supone una velocidad angular constante en la manivela, puede calcularse para ángulos t₁ conocidos, a un paso conocido,y luego graficar cada resultado del sistema para cada valor de angulo t₁.

Las instrucciones anterioes se indican en el siguiente código:

% animacion
clear all;

% paso
paso=6*pi/180;

% numero de puntos
nc=100;

% datos
t1=pi/4;

% estimaciones
x=70;
t2=30*pi/180;

% definicion de los vecotres para guardar los puntos
X=[0,0,0];
Y=X;
XP=[0,0,0,0,0];
YP=XP;

% calculo para cada posicion de t1
for i=1:nc
    [xx,yy,x,t2] = blogpiston(t1,x,t2);
    X(i,:)=xx;
    Y(i,:)=yy;
    t1=t1+paso;
    %cuadrado para simular piston (Es para que se vea bonito )
    [xx,yy]=cuad(xx(3),yy(3),5,0);
    XP(i,:)=xx;
    YP(i,:)=yy;
    %
end

% grafico/animacion
for i=1:nc

    %Marco para la pantalla
    plot(-100,-100)
    hold on
    plot(100,100)
    
    %malla
    plot(X(i,:),Y(i,:),'-o')
    %cuadrado piston
    plot(XP(i,:),YP(i,:))
    %tiempode pausa
    pause(.1)
    hold off

end

Lo anterior puede guardarse como imágenes y puede elaborarse un gif de un pistón agregando la sintaxis del artículo citado, y quedará una animación como la que se ve en la cabeza del artículo.

Dibujar el pistón

Como la "animación" descrita dibuja solo lineas y puntos, agregué una función que dibuja un cuadrado para representar el pistón, tomando como referencia el punto calculado para la posición de este

Notar que x e y son las coordenadas del pistón. El pistón es un cuadrado de lado e que recorre una línea horientada con la vertical de ángulo a

function [X,Y]=cuad(x,y,e,a);
r=e*sqrt(2);
b=pi/4-a;
X=[x-r*sin(b),x+r*cos(b),x+r*sin(b),x-r*cos(b),x-r*sin(b)];
Y=[y-r*cos(b),y-r*sin(b),y+r*cos(b),y+r*sin(b),y-r*cos(b)];

Archivos *.m

Dada cualquier eventualidad los archivos *.m están en los siguientes enlaces:

Función

Programa

Cuadrado para el pistón

Francisco Roco

Volar...volar siempre fue el sueño de la humanidad

2006-08-09T17:53:00.001-04:00

Volar...volar siempre fue el sueño de la humanidad, y ahora que la mente humana puede reproducir este milagro, es visto por todos como tal. Es cosa de fijarse en expresiones que se refieren a sensaciones, como por ejemplo por efecto de la marihuana, "uhh siento volar...arghh",O simplemente al ver alguien ensimismado en su labor hacemos referencia a esta actividad "uhh el loco se voló tocando el piano" .En el enamoramiento es común que se haga analogía con la sensación de vuelo, es clásica la imagen de la pareja volando sobre un campo de flores; tambien escribió Neruda

"Hay que volar a cada instante como
las águilas, las moscas y los días,
hay que vencer los ojos de Saturno
y establecer allí nuevas campanas."

A lo largo de la historia han aparecido varios dioses que poseen alas, como una forma de ilustrar la percepción del hombre sobre ello, como también los ángeles, Gokú y Supermán que está tan de moda en estos días.

Una vez un compañero me dijo "Tu no puedes comparar el volar en avión, con esa sensación de libertad que sería poder volar sin cabinas, sin parabrisas, sin nada..."

Bueno, en parte tiene razón pero aún en la conciencia colectiva las máquinas voladoras tienen una imagen de modernidad, complejidad y de progreso. Si no, amigos míos, ¿Por qué creen ustedes que en el comercial de las toallas higiénicas mencionan las alas de esta y recalcan , y vuelven a decir "con alas"?, ja ja, no eso es estar loco. Un caso más aterrizado es ver en el comercial de Ferrari "La gasolina de la escuderia Ferrari no podía ser otra" donde sale un avión cisterna abasteciendo de combustible al modelo de la Fórmula uno mientras este esta corriendo a gran velocidad.

Como que me volé ¿ah?. Eso tiene explicación en mi gusto por el rubro aeroespacial, eso de disfrutar viendo despegar un avión , o el colosal transbordador espacial, nace de una pregunta que todos nos hacemos cuando chico: ¿Cómo lo hacen?

La respuesta a esta pregunta, obviamente es dificil entenderla a corta edad - para los que están en dicha situación http://www.aviacion.cl/ninos.htm - y a algunas personas nos tiene todavía buscándola, Leonardo Da Vinci quién también gustaba del tema dijo:

"Un pájaro es una máquina que funciona conforme a unas leyes matemáticas , y entra en la capacidad del hombre reproducir esa máquina con todos sus movimientos."

Así que a estudiar nos vamos para continuar con este arte...