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

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.

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

Roco

Gráficos en GNU Octave, la alternativa a Matlab

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,

Roco

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

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