Código de trayectoria de partículas en movimiento browniano

Simulación de la trayectoria de una partícula en  movimiento browniano

En el siguiente ejemplo se simularemos la trayectoria de partículas brownianas por medio de un método aleatorio. El valor del método consiste en utilizar los números aleatorios para generar un arreglo de ángulos aleatorios para el movimiento de la partícula en el espacio.

La forma en la que la partícula se va a mover es simple: se propone un paso unitario en una dirección aleatoria y en esa nueva ubicación se hará otro desplazamiento aleatorio sin que tenga relación con la anterior. El código implementado esta abajo de estos parrafos.

%Rutina para simular la trayectoria de una partícula browniana
% Por Vicente Torres Zúñiga
clear % limpiamos memoria
close all
clc
par=1; % número de particulas,
%soporta hasta 21 particulas para diferenciar el color
% espacio para los valores finales de posicion de las partículas
X=zeros(par,1000);
Y=zeros(par,1000);
Z=zeros(par,1000);
for n=1:par %número de particulas en cada interación
x=zeros(1,1000); %vector de posicion de la particula
y=zeros(1,1000);
z=zeros(1,1000);
%se define una dirección en 3D aleatoria.
phi=2*pi*rand(1,1000);%crea un ángulo aleatorio para la particula
theta=2*pi*rand(1,1000)-pi; %se crea el otro ángulo.
% bastan con dos angulos para cubrir el espacio
for i=1:999 %se define la interacción de la partícula
x(i+1)=x(i)+cos(phi(i)).*sin(theta(i));
y(i+1)=y(i)+sin(phi(i)).*sin(theta(i));
z(i+1)=z(i)+cos(theta(i));
end
%se guarda la trayectoria instantanea de cada partícul
X(n,:)=x;
Y(n,:)=y;
Z(n,:)=z;
end
% vector para cambiar de color en cada trayectoria
C=['r' 'g' 'b' 'w' 'y' 'm' 'c' 'r' 'g' 'b' 'w' 'y' 'm' 'c' 'r' 'g' 'b' 'w' 'y' 'm' 'c'];
%graficación de la trayectoria de las partículas
hold on
for n=1:par;
plot3(X(n,:),Y(n,:),Z(n,:),C(n));
% en lugar de graficar las trayectorias de cada particula,
% se grafica en esa posición el lugar que ocupa la partícula
end
view(3);
axis tight % Ajustes en la visualización de la ventana
box on
hold off

Ejercicios:

1) Muestre que este tipo de trayectorias obedece a una distribución estadística de Gauss.
2) Ejecute el programa con hasta 21 partículas. Comente lo que observa.
3) Este script se implento en base al ángulo, escriba un script que emplee desplazamientos en el plano para simular una trayectoria aleatoria.

Muchas gracias a todos por ser parte de este curso de Matlab


Espero que el curso les gustara, que les sea útil en sus próximos trabajos y empresas.

Seguiré actualizando y colocando entradas en el blog para responder algunas dudas. Hasta enero del 2013 tendrán acceso al blog, pero si tienen dudas pueden contactarme en la Internet.

Gracias :D

Ejemplo del uso de la técnica Montecarlo en una función y su cambio de variable

Este es un ejemplo concreto del uso de técnica Montecarlo con una función y su cambio de variable.  la Integral original es:



El cambio de variable es:
 

De modo que la integral después de realizar el cambio de variable es:



Cuando graficamos la función original y la del cambio de variable, ver la Fig 1, encontramos que el cambio de variable permite obtener una función mucho más suave en el intervalo de integración.

Figura 1. Gráficas de funciones: original y con cambio de variable. La función original presenta más estructura, mientras que el cambio de variable muestra una función casi constante en el intervalo de trabajo

Ahora bien, cuando aplicamos la técnica Montecarlo varias veces a las dos funciones, observamos que las funciones convergen a un valor de integral, ver Fig. 2. Es decir, muestran estabilidad para presentar una respuesta única en un intervalo de confianza. Adicionalmente, encontramos que la función con cambio de variable converge mucho más rápido que la función original, lo cual evidencia la importancia de realizar un adecuado cambio de variable.

Fig. 2.  Graficas de los cálculos de integral por método Montecarlo (con 100 interacciones) de la función original y la de cambio de variable. Note que alcanza más rápido la estabilidad la integral con el cambio de variable.


Valor de la integral original 0.4520, con un error porcentual asociado 0.0143%. En contraste, la integral con el cambio de variable presenta un valor promedio de -0.2948 con un error porcentual = 0.0077% (Mucho menor es el error porque el proceso alcanza rápidamente la estabilidad).

% incia rutina, begins script
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function montecarloevalua
%Rutina para graficar una integral por el metodo de MonteCarlo,
%calcular el tiempo de ejecución, y contar el con la estabilidad
clc % parte donde se limpia la pantalla y la memoria
close all
clear
a=exp(-pi^2); % límite inferior de la integral
b=1; % límite superior de la integral
x1=linspace(a,b,1000); % vector donde ser realiza la integral
y1=1/2*(sin(exp(-x1.^2)))./exp(-x1.^2); % función a integrar
plot(x1,y1) % gráfica de la función original a integrar
figure
avg=[0 0]; %donde se va a contener y
%inicia método montecarlo
tic % inicia contador del reloj
iteracion=100;% numero de veces que se repite la interacción MonteCarlo,
for i=1:10 % Número de graficas para saber la estabilidad del cálculo
for N = 1:iteracion; %iteracion dentro del montecarolo
x1 = a + (b-a)*(rand(1000,1));
term= sum([(1/2*(sin(exp(-x1.^2)))./exp(-x1.^2)) (1/2*(sin(exp(-x1.^2)))./exp(-x1.^2)).^2])/1000;
avg= ((N-1)*avg + term)/N;
z(N)=(b-a)*avg(1);
end
%fin estricto del proceso MonteCarlo.
prom(i) = (b-a)*avg(1); % último punto calculado del proceso MonteCarlo
N = 1000*(1:iteracion); % vector para graficar
V = [z'];
plot(N,V);
hold on;
xlabel('N')
end
toc % finaliza el contador del reloj
promedio_int = mean(prom) %promedio de las integrales calculadas
error = std(prom); %error asociado al promedio de la integral
error_porcentual = 100*error/promedio_int % error porcentual
% termina la rutina, end script
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%inicia método montecarlo tic % inicia contador del reloj iteracion=100;% numero de veces que se repite la interacción MonteCarlo,
for i=1:10 % Número de graficas para saber la estabilidad del cálculo
for N = 1:iteracion; %iteracion dentro del montecarolo x1 = a + (b-a)*(rand(1000,1)); term= sum([(1/2*(sin(exp(-x1.^2)))./exp(-x1.^2)) (1/2*(sin(exp(-x1.^2)))./exp(-x1.^2)).^2])/1000; avg= ((N-1)*avg + term)/N; z(N)=(b-a)*avg(1); end %fin estricto del proceso MonteCarlo. prom(i) = (b-a)*avg(1); % último punto calculado del proceso MonteCarlo N = 1000*(1:iteracion); % vector para graficar
V = [z']; plot(N,V); hold on; xlabel('N') end toc % finaliza el contador del reloj promedio_int = mean(prom) %promedio de las integrales calculadas error = std(prom); %error asociado al promedio de la integral error_porcentual = 100*error/promedio_int % error porcentual % termina la rutina, end script %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Sustituyendo en la rutina las siguientes líneas
a=exp(-exp(-2*pi^2));
b=1/exp(1);
y1=-1/4*(sin(x1)./((x1.^2).*sqrt(-log(x1))));
Ejercicios:

1) Proponga otras funciones y su cambio de variable para repetir este ejercicio
2) Cómo se puede hacer más rápido esta algoritmo
3) Modifique el algoritmo para contar con una función que se alimente con la función a integrar y los limites de integración

Presentación: Introducción a la técnica Montecarlo con 2 ejemplos en Matlab


De acuerdo con la Wikipedia, el uso de los métodos de Montecarlo como herramienta de investigación, proviene del trabajo realizado en el desarrollo de la bomba atómica durante la segunda guerra mundial en el Laboratorio Nacional de Los Álamos en EE.UU.

Los métodos Montecarlo son una herramienta computacional que utilizan las propiedades deterministas de un sistema en conjunto con variables aleatorias para aproximar un resultado, es decir, se basa en un muestreo repetitivo aleatorio. En general se utiliza en problemas donde la evaluación analítica no es posible, o resulta muy complicada; en el estudio de sistemas con un gran número de grados de libertad acoplados, como son los fluidos, estructuras celulares, etc.

Más en general, los métodos de Montecarlo pueden ser útiles para modelar fenómenos con entradas de incertidumbre significantes, como en el cálculo de riesgo en los negocios.

¿Dónde más haz escuchado usar el método Montecarlo?

20 puntos extra: por mejorar interface y programa de captura de imagenes

Interface de video y analisis de una imagen capturada 

Este es un script  para obtener el video (con su archivo fig), la imagen congelada y estadística de esa imagen, funciona con la interface de una cámara web. Funciona en base a la captura de video, cuyo código se encuentra en nuestra práctica de grabación y manipulación de video por medio de Matlab.

En primer equipo que entregue este código con una función en la interface para variar y fijar la sensibilidad de la cámara web. Obtiene hasta 20 puntos (para cada uno de los miembros del equipo) en las evaluaciones de la clase.

El máximo de integrantes por equipo es 3.

La fecha límite para hacer el código es el 30 de noviembre.

Preguntas sobre este proyecto sólo serán contestadas en este post y en la clase presencial.


Links útilies:


Ayuda de Matlab Adq. Toolbox

Cómo simular en Matlab un patrón de difracción de una apertura circular

Apertura circular  y su patrón de difracción simulados con Matlab 

El siguiente guion es para usarse en Matlab y realizar una simulación de un patrón de difracción. Puede ser útil para motivarlos a estudiar más este tema:

%inicia script
%% script para simular un patrón de difración de una apertura circular
%% mediante la transfotmada discreta y rapida de Fourier
%
%% por Vicente Torres (18/8/2011), en base al script de Luis Mex

%limpieza de memoria
clear;
clc;
close all
%acondicionamiento
n=2^10;
M=zeros(n); % matriz base donde, que representa la pantalla, oscura
I=1:n;
x=I-n/2;
y=n/2-I;
[x,y]=meshgrid(x,y);
R=10;
a=(x.^2+y.^2<=R.^2); % acondicionando el área de la apertura,
% estos seran indices de la matriz
M(a)=1; % la apertura deja pasar la luz, que se representa con 1
figure
subplot(1,2,1)
imagesc(M),colormap( [0 0 0;1 1 1]),title('\bf abertura circular'),axis image;
D=fft2(M); % Calculo del patron de difraccion
D1=fftshift(D);
subplot(1,2,2)
imagesc(abs(D1)),colormap(hot),axis image, title('\bf patrón de difracción');
%fin del script

En Matlab central pueden encontrar otros ejemplos de patrones de difracción, basado en la ayuda de Matlab.


Preguntas para pensar:
La forma del patrón de difracción en campo lejano  parece concordar con el experimento, como se puede comprobar que su perfil de intensidad corresponde adecuadamente con la teoría o con el experimento. 

Ejercicios:
1) En base a este ejemplo obtenga el patrón de difracción de una rendija simple.
2) Obtenga la simulación de un patrón de difracción  de una rendija cuadrada

Cómo hacer un diagrama de barras estadístico con Matlab

Existen varias posibilidades para representar diagramas de barras. Por ejemplo, los datos se escriben como un vector y son introducidos a Matlab:

>> x = [10 2 3 4 18 20 15];

Entonces, podemos usar los comandos bar, barh, bar3 y bar3h para generar los gráficos. Usamos como apoyo el comando subplot para ver en una sola imagen todos estos resultados.

>> subplot(2,2,1), bar(x), title('Barras verticales')
>> subplot(2,2,2), barh(x), title('Barras horizontales')
>> subplot(2,2,3), bar3(x), title('Barras verticales 3D')
>> subplot(2,2,4), bar3h(x), title('Barras horizontales 3D')

El resultado lo podemos ver en la siguiente figura:

Diagrama de barras

Por supuesto,  las graficas 3D se pueden modificar empleando el comando rotate3d.

Además, los datos pueden estar agrupados. Como ejemplo, consideremos los siguientes datos

x = [1 2 3; 4 3 6; 10 9 8; 4 2 7; 12 10 7 ];

Utilizando los mismos comandos que antes, con modificaciones como: bar3(x,'group') y bar3(x,'steack'). Obtenemos la siguiente figura.

Diagrama de barras apiladas y agrupados 3D

Preguntas para pensar.

En que casos prácticos son útiles matrices de MxN mostradas como diagrama de barras.

Ejercicios:

Empleando la función rand genere un vector de 1000 datos,  agrupe los valores entre 0-10, 10-20, 20-30 y así.  ahora haga el diagrama de barraras con la agrupación, por lo que solo aparecernos 10 barras.

Código Matlab para acondicionar y filtrar señales de electrocardiogramas.


Sergey Chernenko nos muestra, explica y regala el código Matlab para detectar una característica importante de las señales provenientes de un electrocardiógrafo: el pico R de un electrocardiograma.

Para estudiar este excelente ejemplo de procesamiento de señales en ciencias medicas, lo mejor es bajar el archivo zip, con el código original.

Hay que estudiarlo para comentarlo en la clase.

Como nos comenta el mismo Sergey en su sitio Web, la tarea básica en el procesamiento de electro-cardiogramas (ECG) es la detección de los picos R. Sin embargo, debido a la respiración del paciente, la detección es complicada; pues se presenta en los ECGs artefactos que complican la detección de los picos R. Por ejemplo, las distancias irregulares entre los picos, la forma irregular en los picos, la presencia de componentes de baja frecuencia en los ECGs.

Es por ello que se deben contar con varias etapas de procesamiento que reduzcan la influencia de estos factores. Sergey nos muestra un caso donde una señal ECG original se le eliminan las frecuencias bajas por medio de FFT (Fast-Fourier-Transform). En la segunda etapa, se usa un filtro pasa banda, el cual localiza los máximos significativos en la señal. En la etapa final, se eliminan los datos más alejados del promedio. Este proceso permite la detección de los picos R.

Código para obtener la primera derivada por coeficientes de una función

Como se puede consultar en el libro: transformada de Fourier para peatones, podemos definir la derivada Yk como la “primera diferencia central” de la función F:


La cual tiene una función de transferencia, que no presenta problemas de corrimientos de fase:


Sin embargo, presenta un efecto de filtro muy pronunciado a frecuencias altas.

El siguiente código matlab aplicado a la función sin^2(x) con un paso de 0.1 funciona muy bien. En las figuras de salida del script se observan la función original, la derivada teórica y la derivada obtenida. Observe que presenta un comportamiento muy similar las dos derivadas. Efectivamente, cuando comparamos la derivada obtenida sobre la deriva teórica obtenemos en la gráfica una recta identidad. Por todo lo anterior, afirmamos que es un buen algoritmo para operar la derivada.

%%Script que calcúla numéricamente la primera derivada de la función:
%%sin(x)^2.
clear all
close all
clc
%Definiendo intervalos
x=0:0.01:1;
%función original:
funcion=(sin(x)).^2;
%preparando parámetros:
s=length(funcion); %longitud de la funcion
Y=zeros(1,s); % espacio de memoria para la derivada numerica
% Algoritmo de primera derivada
for k=2:s-1
Y(k)=(funcion(k+1)-funcion(k-1))/(2*0.01);
end
% Graficación de funciones teoricas y su compraración con la estimación
% numerica de la primer derivada
figure
subplot(2,1,1)
hold on
plot(x,funcion,'r') % función original
plot(x,2*sin(x).*cos(x), '-b') % derivada teorica de la funcion original
plot(x,Y,'ok') % derivada numerica obtenida
axis([0,1.04,0,1.1]);

xlabel('Vector <>')

ylabel('Funciones utilizadas')
text(0.02,0.8,'Funciones Sin^2(x) (rojo), derivada teorica (azul), derivada numerica (negro)')
box on
hold off 
subplot(2,1,2)
plot(2*sin(x).*cos(x),Y,'-b')% comparacion entre funciones;
% derivada teorica contra la derivada numerica
xlabel('Función 2*sin(x).*cos(x), derivada teorica')
ylabel('Derivada numerica obtenida')
axis([0,1.04,0,1.1]);
grid
%%% fin del script 


Ejercicios:

Matlab tiene su comando/función para hacer la derrivada: diff, su algoritmo de trabajo es:


1) ¿Cómo se compara este algoritmo con el propuesto en el post?

2) Muestre ejemplos del uso de este comando

3) ¿Cuales son las limitaciones de esta forma de definir la derivada?

Enlaces relacionados:

Capítulo 2, Derivación numérica del curso de métodos numéricos de Agustín A. Rosas.

Comparación del desempeño de un filtro convolución y un filtro pasa bajas de promedios ponderados en una misma función/señal

Ahora mostramos dos ejemplos de filtros. La funcion a trabajar es:




La cual tiene un pico extendido en t=5 y dos pequeños picos de estructura en t=2. Pues bien la entrada a un instrumento  cuenta con  filtro descrito por:




Efectivamente, la señal filtrada por  convolución será la transforma de Fourier inversa de la multiplicación de las transformadas de Fourier de la señal y del filtro; es decir:



A continuación aplicaremos el filtro en código Matlab.

Filtro 1: Convolución.

Para comparar la señal y la salida filtrada, a ambas las normalizamos. Observamos que hay un corrimiento en el tiempo. Este efecto en los indices y en la amplitud son muy comunes cuando se usa la transformada de Fourier.



Adicionalmente, una vista cercana a la estructura de los dos picos pequeños muestra que la salida del filtro cuenta con estos picos más pequeños y casi sin desernir. Con todo, esta forma de presentar los datos muestra a detalle que el filtro es bueno sin ser excelente para eliminar la estructura de la curva principal .


A continuación, mostramos el código Matlab utilizado para hacer estas curvas.


 %%% Inicio código 1: Filtro convolución
 clc
clear
close all

x = 0:0.0001:10;

y = 4*((sin(x-2)).^2).*exp(-50*((x-2)).^2) + 0.5*( ((x-5).^2 +0.25).^(-.25) );

r = ones(size(x));
for k= 1:length(x)
    if x(k) <= 0.15
        r(k)= .50;
    elseif x(k) <= 0.30
        r(k)= 0.25;
    elseif x(k) <= 0.60
        r(k)= 0.10;
    else
        r(k)= 0;
    end
end
plot(x,r,  'LineWidth',3)
axis([-.2,3,-0.2,0.6])
Y= ifft( fft(y).*fft(r));
figure
hold on
plot(x, y/max(y), '--b','LineWidth',3)
plot(x,Y/max(Y), 'k', 'LineWidth',3)
axis([0,10,0.28,1.05])
hold off
%%% Fin código 1: Filtro convolución

Filtro 2: Pasa bajas con promedio ponderado


Ahora bien, utilizaremos el filtro pasa bajas (PB), que ya mostramos en otro post, en este caso usamos 6,0000 iteraciones. Dado que la estructura en la función puede representar una oscilación rápida aplicaremos el filtro PB.   Observamos en la salida que tanto el cambio en la amplitud como el corrimiento en los indices de los vectores de las curvas es casi minimo. Además de que se ve más disminuida la estructura en la función.



Este filtro PB presenta mejores resultados que el de convolución; esto es por el tipo de filtro r(t) utilizado. Sin embargo, el filtro PB  tarda mucho más en desarrollarse. Puedes  usar los comandos de Matlab TIC yTOC antes y despues de que inicien las rutinas para observar la velocidad de ejecución de cada función o escript.


Podemos asegurar que dependiendo que se busquemos filtrar, hay un filtro que nos funcionara mejor que otros, pero hay que saberlo encontrar :)

A continuacion el segundo código Matlab:

%%% Inicio código 2: filtro pasa bajas con promedio ponderado
clc
clear
close all

tic % inicia el contador de tiempo
t = 0:0.001:5*pi;
f = 4*((sin(t-2)).^2).*exp((-50)*(t-2).^2) +0.5*((t-5).^2+0.25).^(-0.25); %señal de trabajo

hold on
plot(t, f, 'r') % para visualizar la señal con ruido

for corrida = 1:60000 %iteracion del filtro
for k = 2:(length(t)) -1 ;
f(k) = 0.25*(f(k-1) + 2*f(k) + f(k+1)); %filtro
end
end

plot(t, f, '--b') %para visualizar la señal filtrada
hold off

toc %finaliza el contador del tiempo

%%% Fin código 2: filtro pasa bajas con promedio ponderado

Ejercicios: 

1) Muestre el grado de atenuación de las protuberancias después de pasar por el filtro.

2) Diseñe otro filtro de convolución que atenue mejor a los dos picos que forman protuberancias en la curva original

3) Calcule la salida del filtro de convolución con r(w)*I(w) en lugar de  I(w)*r(w). Describa sus observaciones.

Ejemplos de código para formar señales discretas en Matlab

El objeto más básico en Matlab es una matriz numérica con la posibilidad de almacenar números complejos. Por supuesto, los datos obtenidos en el estudio de señales y sistemas son muy bien representados en forma de matrices. En este post usaremos Matlab para la generación de señales elementales: cuadrada, triangular, entre otras.

El ToolBox de procesamiento de señales de Matlab posee una gran variedad de funciones para la generación de señales, estas señales requieren de una representación vectorial de la variable tiempo, de manera continua o discreta. Para realizar una simulación de un intervalo continuo, se usa un vector de valores discretos con un intervalo de muestreo muy pequeño.

Como vimos en post anteriores, el siguiente comando genera un vector llamado t de valores que representan la variable tiempo, con un intervalo de muestreo de 1 ms entre 0 y 1 segundo.

t = 0:0.001:1;

Después de creado el vector que representa la variable tiempo, es posible iniciar el desarrollo de alguna señal de interés.

En Matlab una señal discreta en el tiempo se representa exactamente, porque los valores de la señal son representados como los elementos de un vector. Sin embargo las señales de tiempo continuo en Matlab son tan solo aproximaciones. La aproximación consiste de un vector cuyos elementos son muestras de la verdadera señal de tiempo continuo. Cuando se usa esta técnica para la representación de señales continuas es importante escoger el intervalo de muestreo lo suficientemente pequeño para asegurar que las muestras capturan todos los detalles de la señal.

EJEMPLOS DE SEÑALES EN MATLAB 

La generación de señales periódicas tales como ondas cuadradas y triangulares es una actividad muy fácil de realizar en MATLAB.

1) SEÑAL CUADRADA
Consideremos primero la generación de una onda cuadrada de amplitud A, frecuencia fundamental w (medida en radianes por segundo) y ciclo útil rho. Recordemos que el ciclo útil es la fracción de cada periodo en donde la señal es positiva.


Para generar dicha señal se pueden escribir lo siguiente en la linea de comandos:


>> A = 1;
>> w = 10 * pi;
>> rho = 0.5;
>> t = 0:0.001:1;
>> sq = A*square(w*t+rho);
>> plot(t,sq);

En la segunda línea de este ejemplo, pi es una función interna de Matlab que calcula el número más cercano a la constante PI en formato de coma flotante. El último comando es usado para vizualizar la señal generada. El comando plot dibuja líneas conectando los valores sucesivos de la señal y así da la apariencia de una señal en tiempo continuo.

2) SEÑAL TRINGULAR 

Consideremos ahora la generación de una onda triangular de amplitud A, frecuencia fundamental w y ancho Wdt . El periodo de la onda triangular será T con el máximo valor de la señal ocurriendo en t = WT . El comando básico para generar esta señal es:

A * sawtooth(w * t + Wdt)

El resultado se puede observar en la gráfica a la izquierda


3) SEÑAL ESCALÓN

En Matlab, el comando ones(M, N) genera una matriz de unos de tamaño MxN, y el comando zeros(M, N) es una matriz de ceros del mismo tamaño. Se puede hacer uso de estas dos matrices para generar dos señales comúnmente usadas: la señal escalón y la señal impulso.

Una señal paso de amplitud uno, puede ser generada con el siguiente comando.


U = [zeros(1, 10), ones(1, 11)];


Para la versión continua creamos un vector que represente el tiempo el cual tenga muestras de un intervalo separados por valores muy pequeños. Los comandos y los resultados se muestran a continuación:

>> u=[zeros(1,10),ones(1,11);
>> t=-1:0.1:1;
>> plot(t,u)


Como se menciono anteriormente, una señal generada en Matlab es inherentemente de naturaleza discreta. Para visualizar una señal en tiempo discreto se puede hacer uso del comando stem. Específicamente stem(n, x), bosqueja los datos contenidos en el vector x como una señal de tiempo discreto con los valores de tiempo definidos por el vector n. Los vectores n y x deben tener dimensiones compatibles, es decir deben tener el mismo número de elementos. Así, para este caso para obtener la representación de esta señal en tiempo discreto creamos un vector-tiempo el cual debe tener valores separados por una unidad.

>> u=[zeros(1,10), ones(1,11)];
>> n=-10:10;
>> stem(n,u)

Recuerde que para poder usar las funciones plot y stem, es requisito que los vectores (t y u) ó (n y u) tengan iguales dimensiones. Por esta razón el vector u se forma como una composición de diez ceros y 11 unos, debido a que los arreglos t y n, tienen dimensión 21 dado que incluyen un elemento central el cual es el número cero. Para probar este hecho, se puede hacer uso de la función Matlab llamada size que devuelve como resultado un vector con las dimensiones de la matriz que se le pasa como parámetro.

4) SEÑAL IMPULSO:


La versión discreta de la señal impulso se puede también generar con ayuda de las funciones zeros y ones, realizando una composición como sigue:


>> delta = [ zeros( 1 ,10 ), 1 , zeros( 1 ,10 ) ];
>> n = -10:10;
>> stem(n,delta);



5) SEÑAL RAMPA 

Para generar la señal rampa, tan solo es necesario recordar que esta función puede ser creada, como la composición de una recta Y(x) = x a partir de cero y de la recta Y(x) = 0 para valores de x menores de cero, así la versión discreta se muestra a continuación:


>> t1=0:0.1:10;
>> rampa1=t1;
>> rampa=[zeros(1,101),rampa1];
>> t2=-10:0.1:0;
>> t=[t2,t1];
>> plot(t,rampa)

Ejercicios:

1) Desarrollar un conjunto de comandos Matlab para aproximar las siguientes señales periódicas en tiempo continuo, dibujando 5 ciclos de cada una:

a) Onda Cuadrada, de amplitud 5 Volts, frecuencia fundamental 20 Hz y ciclo útil del 60%.

b) Señal diente de sierra, amplitud 5 Volts y frecuencia fundamental 20Hz

Práctica 2: Señales harmónicas moduladas y su filtrado

Ya pueden consultar su práctica 2, para lo cual les recomiendo que leen sus libros sobre transformada de Fourier.

Por cierto, esta presentación sobre función de transferencia les será útil en los siguientes semanas.









Por cierto, también deben consultar el artículo: How can Trigonometry help to understand Lock-in Amplifier operation? De Ernesto Marín, R. Ivanov. Pues la técnica lock-in también es un método de filtrado de señales, muy aplicado cuando la razón señal/ruido es muy baja.

Código y ejemplos con diferentes iteraciones de un filtro pasa bajas de promedio ponderado


Estas dos figuras se obtuvieron con pequeñas modificaciones del código que aquí se presenta. En el primer caso se realizaron 10,000 iteraciones del mismo filtro. El tiempo que tomo fue considerable y su filtradrado fue bueno. Sin embargo con 1,000 iteracciones el filtrado es malo y el tiempo no se reduce considerablemente, lo cual es lógico pues disminuyo en un factor de 10 las iteraciones, aunque siguen siendo muchas iteraciones 1,000. Estos son sólo ejemplos de cómo aplicar filtros de promedios ponderados hay filtros alternativos que son más rápidos y eficaces en la tarea asignada. El ejercicio sirve para ver las entrañas de estos filtros.


% Inicio de código
clc
clear
close all

x = 0.0:0.001:3*pi;

yy = sin(0.5*x)+sin(40*x); % la frecuencia a limpiar es la de 40, la cual es relativamente grande

subplot(3,1,1)
plot(x,yy) % primero vemos como es la señal original

% limpiar señal

% y(k) = yy(k-1) + 2*yy(k) + yy(k+1)

final = length(yy);

y = yy;
for m =1:10000 % hacemos muchas iteraciónes pues el filtro es poco eficiente

for k =2:final-1
y(k) = (0.25)*(y(k-1) + 2*y(k) + y(k+1)); % filtro pasa bajas
end
end

subplot(3,1,2)

plot(x,y) % el resultado lo vemos aqui.

subplot(3,1,3)
hold on
plot(sin(0.5*x),y) % el resultado lo vemos aqui.
plot(y,y, '--k') % el resultado lo vemos aqui.

hold off

Código para encontrar la señal principal mediante el filtro de la función max

Les dejo uno de los ejemplos que me enviaron para hacer la evaluación.


% Inicio  del código
%% ejemplo del uso de las funciones FFT y MAX para limpiar una señal

%% con ruido
%% 12 sept 2011

% limpieza de la memoria
clc
clear
close all

% Señales con ruido
x=0:1/255:1;
y= 10*sin(2*pi*29*x)+  ...
   3*sin(2*5*x)+ ...
   2.5*sin(rand*25*2*pi*x)+...
   3.5*cos(rand*155*2*pi*x)+...
   2.5*sin(300*pi*x);

% Transformada de Fourier simple
Y=fft(y);
[a,b]=max(Y);
%indiceAmplitudmax=b;
% haciendo el vector con ceros
Y1=Y*0;
% sustituyendo en el espacio del indice el numero complejo que representa
% el maximo
Y1(b)=(a);
%tomando unicamente la parte real de la parte inversa de la transformada de
%Fourier. La parte imaginaria no nos da más información en esta parte
inversa= real(ifft(Y1));

% visualizando en la misma grafica la señal con ruido y limpia
hold on
plot(x,y,'r');
plot(x,inversa, 'k')
hold off

% Fin del código

Pregunta para pensar:

¿qué implica en la señal de principal (color negro) que sea un número complejo el máximo del espectro en la Transformada de Fourier?,  ¿acaso son son dos señales?  o ¿cambios de fase?.

Códigos adicionales para gráficas y código de filtro pasabajas para mejorarlo

Les pongo el Código para la tarea, deben hacer que funcione el filtro.



%%%%%%%%%%%%%%%%%%
clc
clear
close all

%señal a trabajar

x = 0:0.0001:5*pi;
 y = sin(x) + sin(100*x);

 subplot(2,1,1)
 plot(x, y, 'c')


 for corrida = 1:6
     for k = 2:(length(x)) -1 ;
   
 y(k)  = 0.25*(y(k-1) + 2*y(k) + y(k+1));
     end
 end

 subplot(2,1,2)
  plot(x, y, 'y')
%%%%%%%%%%%%%%%%%%

Por otro lado, estos son los comandos extra para hacer otro tipo de gráficas:

X = [linspace(0,pi*3/8,20); linspace(pi*5/8, pi, 20)];
Y = sin(X); % Tanto la X como la Y son matrices
plot(X,Y) % Esto no es lo que se esperaba!!! Hay que teneren cuenta que se imprime por columnas


pie(X)
pie(X, Explode)

x = [1 1 2 3 5 8];
ex = [0 1 0 0 1 0];
pie(x, ex);

hist(X)
hist(X, nbins)
x = randn(1,10000);
hist(x);

% Generamos números aleatorios con distribución normal y
uniforme
x = [randn(10000, 1) rand(10000, 1)*6-3]
hist(x, 50);

x = 0.5:0.5:4;
y = 1./x;
barh(x, y)

Barras apiladas

x = 0.5:0.5:4;
y = 1./x;
Y = [y' fliplr(y)'];
bar(x, Y)

x = 0.5:0.5:4;
y = 1./x;
Y = [y' fliplr(y)'];
bar(x, Y, 'stacked')
x = randn(1,100);
y = randn(1,100);
scatter(x, y)

scatter(x, y, rand(1, 100)*1000, [1 0.5 0])

Práctica 1: Trasformada rápida de Fourier discreta aplicada en simulaciones de ruido

En nuestra primer práctica (documento pdf) recodaremos algunos conceptos sobre trasformada de Fourier discreta (como también el teorema de Nyquist) y los fundamentos de la representación matemática de ondas y ruido.

Actividad.
Lee con cuidado las actividades de la práctica 1. Antes de empezar la clase debes de haber repasado un poco los temas de T. de Fourier.

Espero que todos tengan instalado el Matlab para trabajar en equipo el miércoles 24 de agosto.

Para repasar los temas de esta práctica y de otras sesiones, te recomiendo mucho el libro: Trasformada de Fourier para peatones.

Por favor, si tienen un comentario o duda, pueden decirme por este medio electrónico o en la clase presencial.
Creditos:

Imagen del poster de la película White Noise (2005)

¡Bienvenidos a su blog didáctico de programación en Matlab!

Este blog tiene la intención de apoyarte en la clase de programación en Matlab.
Por ello es importante que realices las actividades online que aquí se incluyen. Además del material adicional que te ayudara en la clase: videos, notas, calificaciones, hojas de cálculo para hacer trabajo colaborativo. Adicionalemte, en este blog encontraras fechas importantes para el curso, notas y otra información importante.
Ciertamente, el blog es un medio de comunicación, por lo cual esperamos que tanto en la clase presencial como en este espacio participes con tus recomendaciones a material didáctico extra que puede ayudarnos a dar un mejor curso.

Recuerda el la clase la hacemos todos =) .

Bienvenido al curso y esperamos que aprendas muchas cosas nuevas, útiles y divertidas ^_^ .

Actividad 1: Dinos quien eres y sabremos cómo hacer la clase mejor



Recuerda enviar la forma para que tengamos la información.

fecha limite para la actividad: 24 de agosto del 2011

Cómo escribir operaciones básicas en Matlab (0001)




Resumen.
Se muestra como se abre el programa Matlab, se hacen operaciones básicas: suma, resta, multiplicación y división de números. También se dan ejemplos de la función del punto y coma, el empleo de básico de variables, dónde se guarda temporalmente la información de una variable. Finalmente, se muestra el uso de comandos clc y exit, los que hay que evitar usar como variables.

Preguntas abiertas.
1) ¿El comando clc borra permanentemente la información que escribí en la línea de comandos del programa?

2) ¿Cómo puedo investigar si una variable es usada previamente por el programa Matlab?

Ejercicios.
Realiza algunas operaciones básicas y después combínalas entre ellas. Comprueba la salida de datos con tus propias cuentas a mano.

Enlaces relacionados.

Cómo escribir fórmulas básicas en Matlab (0002)



Resumen.
Se muestra el uso de los botones “arriba” y “abajo” como atajo para recuperar expresiones pre-escritas en la línea de comando. Ejemplo de variables largas. También se muestra como hacer fórmulas sencillas en una sola línea de comando. Finalmente, se muestra el comando “help”

Preguntas abiertas.

1) ¿Dónde puedo encontrar la información sobre el límite de características que puede tener el nombre de una variable larga?

2) ¿Cuáles son las ventajas de usar nombres de variables largas?, ¿y cuales son sus desventajas?

Ejercicios.

Realiza la fórmula en una sola línea de la ecuación que describe el alcance en un tiro parabólico

Enlaces relacionados.

Cómo escribir matrices en Matlab (0003)



Resumen.
Se muestra cómo escribir en Matlab vectores y matrices, con números complejos. Además, mostramos ejemplos del uso de  índices. Ponemos énfasis que Matlab diferencia entre mayúsculas y minúsculas.

Preguntas abiertas.
1) ¿Se pueden escribir en Matlab matrices de orde superior a tres, es decir matrices del tipo NxMx...xZ?
2) ¿Qué sentido práctico puede tener una matriz de hiperdimensionalidad?

Ejercicios.
Escribe una matriz de 3x3 y utiliza subíndices para construir 3 vectores horizontales y 3 verticales, de modo que todos sean vectores diferentes

Enlaces relacionados.

Matrices como variables en Matlab (0004)



Resumen.
Mostramos como las matrices pueden ser variables en Matlab: ponemos operaciones en matrices, usamos matrices pequeñas para construir matrices más grandes. También hacemos matrices pequeñas derivadas de otras más grandes (e.g. por medio de dos puntos). Finalmente, empleamos CTRL+C para abortar una operación atascada en Matlab

Preguntas abiertas.
1) ¿Cómo se puede ahorrar memoria en nuestras instrucciones de Matlab? Mencione ejemplos.

Ejercicios.
Utiliza los dos puntos para construir una matriz formada de dos matrices originalmente independientes. Utiliza dimensiones de 2x3 y 3x2.

Enlaces relacionados.

Suma de vectores y matrices en Matlab (0005)



Resumen.
Se revisa la suma y resta de matrices y submatrices (con números complejos o reales puros), conforme a las reglas del álgebra lineal y otras propias de Matlab sobre la combinación de dimensiones.

Preguntas abiertas.
1) ¿Se puede considerar que los números añaden una dimensión extra a cualquier matriz?
2) ¿Cómo puedo construir una matriz de orden superior a 3?

Ejercicios.
Construye dos matrices de dimensiones diferentes. Utilizando subíndices y dos puntos, suma elementos de las anteriores matrices para construir nuevas.

Enlaces relacionados.


Multiplicación de matrices (0006)




Resumen.
Hacemos ejemplos de multiplicación de matrices: simple (siguiendo las reglas del algebra lineal) y con punto (elemento por elemento).

Preguntas abiertas.
1) ¿En qué casos tienen una aplicación práctica la multiplicación de matrices elemento por elemento?

2) Las proyecciones de vectores son un ejemplo del uso práctico de la multiplicación de matrices (como dice el algebra lineal). Mencione ejemplos de uso en otras aéreas del conocimiento o de aplicadas.

Ejercicios.
Construye dos matrices (2x2) que puedan conmutar bajo la multiplicación simple. Demuestra su propiedad.

Enlaces relacionados.

División entrada a entrada y división de matrices en Matlab (0007)



Resumen.

Programamos unas cuantas matrices para mostrar la función del comando “/” para hacer divisiones.  Hemos mostrado  cómo hacer divisiones entrada por entrada, una matriz dividida por un escalar, luego el uso de “./” para dividir una matriz o escalar entre una matriz.

Después, aquí mostramos como dividir dos matrices. En algebra lineal NO EXISTE LA DIVISION ENTRE MATRICES. Sin embargo, en Matlab la división tiene la interpretación como la división entre una matriz entre la inversa de la matriz dividida. Es decir, para dos matrices A, B “A/B”, equivale a “A*inv(B)”, donde inv es un comando de Matlab. Por supuesto, cuando las matrices deben tener las mismas dimensiones y la matriz “B” debe presentar inversa; de otro modo, los resultados son erróneos.

Preguntas abiertas.

¿Por qué hay matrices que no tienen inversa?, ¿Cuáles son las características de una matriz para carecer de inversa?

Ejercicios

1) Compruebe, con diferentes ejemplos, que para dos matrices A, B en Matlab “A/B”, equivale a “A*inv(B)”.

2) Escriba dos matrices no cuadradas, y utilice índices para reescribir sub-matrices que si se puedan dividir

Enlaces relacionados

Calendario del curso: sus fechas importantes para estar al día en la clase

Algunos tutoriales extra, online y en español para introducirse a Matlab

En la red encontramos muchas lecturas complementarias para el curso. Las recomendamos pues te pueden dar otra perspectiva para la correcta programación en Matlab, además de que puedes sentirte a gusto consultando estos materiales

1) Aprenda Matlab 7.0 como si estuviera en primero, esta es una guía muy completa

2) Tutorial de Manuel Vargas, un tutorial corto, básico y muy práctico, ideal para quien inicia

Actividad: dejamos un comentario con otra guía de Matlab en español: link y el porqué te grado. Recuerda poner tu nombre y no repetir links.
Con este material complementariemos la información de tutoriales.

Fecha límite para dejar comentario: 1 octubre 2011
Related Posts Plugin for WordPress, Blogger...