Señales aleatorias y filtrado de imágenes

Trabajo realizado por:

• David Cabrero Sañudo
• Ainara Contreras Echebarria
• Roberto Rodríguez Gallego

1 Procesamiento de imágenes (Oversampling)

El proceso de Oversampling consiste en aumentar el número de pixeles de la imagen original.

1.1 Generar un filtro Low pass que contenga 2N + 1 elementos

Filtro Low-Pass

En la imagen de la izquierda vemos 15 puntos ya que el valor de N introducido a sido de 7 (2N + 1 = 15) y en la imagen de la derecha podemos ver el filtro en tres dimensiones.

Filtro Low-Pass

1.2 Aplicar una ventana de Hamming al filtro

Ventana Hamming

Igual que para el filtro de Low-Pass vemos que en la imagen de la izquierda nos aparecen 15 puntos y en la imagen de la derecha vemos la representación en tres dimensiones. Al aplicarle la ventana de Hamming lo que conseguimos es amortiguar las altas frecuencias y amoldarse mejor a las bajas frecuencias.

Ventana Hamming

1.3 Realizar un filtro Low-Pass en dos dimensiones tomando variables separadas para (x) e (y)

La aplicación de un filtro de frecuencia baja tiene el efecto de eliminar frecuencias altas y medias dando como resultado una imagen que tiene un menor contraste, una apariencia más suave. Es por esto que este proceso es también denominado “suavización de imágenes” (image smoothing) y al filtro de frecuencia baja se le llama filtro de suavizado o de homogeneización (smoothing filter).

%N es el número de elementos del filtro
%M es la proporción de frecuencias que queremos eliminar
Nx=7;
Mx=4;
Ny=7;
My=4;
% Filtros para x e y 
 h1=filtrolp(Nx,Mx);
 h2=filtrolp(Ny,My);
% Filtros para x e y (hamming)
 h1m=filtrolph(Nx,Mx);
 h2m=filtrolph(Ny,My);
% filtro dos dimensiones low pass
 subplot(1,2,1) 
 h=h1'*h2;
 surf(h)
 title('low pass')
 % filtro dos dimensiones Hamming
 subplot(1,2,2) 
 hm=h1m'*h2m;
 surf(hm)
 title('Hamming')
 % dtft con low pass
 DTFTh=fftshift(abs(fft2(h,2^6,2^6)));
 n1=2^6; n2=2^6; p1=1/n1; p2=1/n2;
 tau1=[-1/2:p1:1/2-p1];
 tau2=[-1/2:p2:1/2-p2];
 [tt1,tt2]=meshgrid(tau1,tau2)
 figure(2)
 subplot(1,2,1)
 surf(tt1,tt2,DTFTh)
 title('DTFT low pass')
 %dtft con hamming
 DTFThm=fftshift(abs(fft2(hm,2^6,2^6)));
 subplot(1,2,2)
 surf(tt1,tt2,DTFThm)
 title('DTFT Hamming')

DTFT Low-Pass

DTFT ventana Hamming

1.4 Aplicar un factor de ampliación a una imagen

El primer paso será añadir M-1 ceros entre cada dos valores de la muestra inicial y luego aplicarle el filtro de Low-Pass.

function z = interM(imagen,N,M)

[a,b]=size(imagen);
%matriz de ceros
imagen2=zeros(a*M,b*M);
imagen2(1:M:a*M,1:M:b*M)=imagen(1:a,1:b);

% Filtros para x e y 
 h1=filtrolp(N,M);
 h=h1'*h1;
 
imagenfil=filter2(h,imagen2);
imagesc(imagenfil);
colormap('gray'); 
axis('image'); 
title('filtrada  ')


figure(2)
imagesc(imagen2);
colormap('gray'); 
axis('image'); 
title('con ceros  ')

En la imagen vemos como cada dos pixeles de la imagen original a intercalado un valor nulo, ya que hemos utilizado M=2.

Imagen con ceros

Después de aplicar el filtro de Low Pass el resultado es el siguiente:

Imagen con ceros con filtro Low-Pass

Vemos que en los pixeles donde añadió un cero esta más oscuro, pero queda mucho más suavizado que antes de aplicarle el filtro.

1.5 Aplicar un factor de ampliación a una imagen, aplicando un filtro con la ventana de Hamming

function z=interMh(imagen,N,M)
[a,b]=size(imagen);

imagen2=zeros(a*M,b*M);
imagen2(1:M:a*M,1:M:b*M)=imagen(1:a,1:b);

% Filtros para x e y 
 h1=filtrolph(N,M);
 h=h1'*h1;
 
imagenfil=filter2(h,imagen2);
imagesc(imagenfil);
colormap('gray'); 
axis('image'); 
title('filtrada  ')


figure(2)
imagesc(imagen2);
colormap('gray'); 
axis('image'); 
title('con ceros  ')

En la imagen vemos como cada dos pixeles de la imagen original a intercalado un valor nulo, ya que hemos utilizado M=2.

Imagen con ceros

Después de aplicar el filtro de Low Pass con la ventana de Hamming el resultado es el siguiente:

Imagen con ceros con ventana de Hamming

Para comparar los resultados obtenidos mediante la ventana Hamming y sin ella:

Vemos que mediante la ventana Hamming los resultados son peores, se parecen más a la imagen original.

Comparación

1.6 Multiplicar el tamaño de la imagen lena por cuatro y aplicar un filtro con 2N + 1 valores siendo N=10

%leer la imagen
lena=imread('lena.jpg','jpeg');
%multiplicar por 4 la imagen
N=10;
M=4; 
figure (1)
h=1/M;
lena4=lena(1:h:512,1:h:512);
imagesc(lena4);
colormap('gray');  
axis('image'); 
title('imagen x 4')
%aplicar un filtro low pass
theta=pi*[1:N];
h=sin(2*theta/(2*M))./theta; 
h=[fliplr(h) 1/M h];
h=h/(sum(h))*M; 
hb2=h'*h
lenafil=filter2(hb2,lena4);
%anadir ventana Hamming
theta=pi*[1:N];
h=sin(2*theta/(2*M))./theta; 
hh=h.*(0.54+0.46*cos(pi*theta/N));
%simetrizamos
hh=[fliplr(hh) 1/M hh];
hh=hh/(sum(hh))*M; % evitamos que el filtro amplifique la senal
z=hh;
hb3=hh'*hh;
lenafil2=filter2(hb3,lena4);
%dibujar
figure(2)
subplot(1,2,1)
imagesc(lenafil);
colormap('gray'); 
axis('image');
title('low pass')
subplot(1,2,2)
imagesc(lenafil2);
colormap('gray'); 
axis('image'); 
title('hamming')

En la imagen de la izquierda vemos la imagen ampliada, es decir, que paso de tener (512 x 512) pixeles a tener (2048 x 2048) pixeles. En la imagen de la derecha vemos primero la imagen aplicándole el filtro de Low-Pass y en la segunda imagen aplicándole la ventana de Haming.

Imagen ampliada

En la siguiente imagen hemos ampliado una zona para poder apreciar mejor cuales son los resultados del filtro, y vemos que aplicándole el filtro Low-Pass la imagen se ve mejor y al añadirle la ventana de Hamming que lo que hace es suavizar la imagen.

Imagen ampliada

1.7 Realizar lo mismo que en el ejercicio anterior pero con N=1

Vemos que al ser menor el número de elementos del filtro, es decir, el valor de N, la imagen filtrada se parce más a imagen original.

Imagen ampliada

1.8 Ahora ampliamos la imagen por ocho y el número de elementos del filtro lo igualamos a diez

Al aumentar más el tamaño de la imagen, se vera más pixelada pero vemos que manteniendo el número de elementos del filtro el resultado es similar al que obtuvimos al aumentar la imagen por 4.

Imagen ampliada

3 Modelos autoregresivos y de media móvil

3.1 Proceso MA-4

Un proceso aleatorio MA-q se define por: donde w es un ruido blanco de varianza δ² y b₁,...,b_q son los coeficientes.

En nuestro caso vamos a representar un proceso MA-4 con 4 coeficientes con un valor comprendido entre -2 y 2, con lo que este proceso quedaria definido como:

x(k)=w(k)+f₁*w(k-1)+f₂*w(k-2)+f₃*w(k-3)+f₄*w(k-4)

Vamos a realizar este experimento con un diferentes numero de valores (200,300,500 y 1000)

% Proceso MA-4
clear all
% Parametros aleatorios entre -2 y 2
a1=randn(1);
if a1<0
    a1=-1;
else
    a1=1;
end
f1=a1*(2*rand(1))

a2=randn(1);
if a2<0
    a2=-1;
else
    a2=1;
end
f2=a2*(2*rand(1))

a3=randn(1);
if a3<0
    a3=-1;
else
    a3=1;
end
f3=a3*(2*rand(1))

a4=randn(1);
if a4<0
    a4=-1;
else
    a4=1;
end
f4=a4*(2*rand(1))

% Ruido blanco
N=300;
a=rand(1,N);
% Proceso MA4
x(1)=randn(1);
x(2)=randn(1);
x(3)=randn(1);
x(4)=randn(1);

for n=5:N
    x(n)=a(n)+f1*a(n-1)+f2*a(n-2)+f3*a(n-3)+f4*a(n-4);
end

 % Dibujamos el proceso
 figure (1)
 plot(a)
 hold on
 plot(x,'x-r')
 title('N=300')

 % Reproducir el proceso usando filtros
 h=[1 f1 f2 f3 f4];
 y=filter(h,1,a);
 plot(y,'go')        
 legend('Ruido Blanco','Proceso MA-4','Filtro','Location','South')
 
 
  % Calculo de la PSD 
  PSD2=fftshift(abs(fft(y)));
  figure (2)
  z2=PSD2.^2/N;
  tau3=-1/2:1/N:1/2-1/N;
  plot(tau3,log(z2),'rx-')
    
  % Segundo metodo  
  H1=abs(fftshift(fft(h,N))).^2;
  tau2=-1/2:1/N:1/2-1/N;
  hold on
  plot(tau2,log(H1),'b')
  title('N=300')

Proceso de Media Móvil de 4 coeficientes con diferentes valores

Representación de las PSD de los procesos MA-4

Si nos fijamos en el principio de las gráficas, se puede apreciar como el filtro tarda 4 puntos en empezar a actuar. A partir de ahí coincide exactamente con el proceso MA-4

centro

3.2 Proceso Autoregresivo

Un proceso aleatorio AR-q se define por: donde w es un ruido blanco de varianza δ² y a₁,...,a_q son los coeficientes.

Construiremos un proceso x(t) con parámetros ρ=0.95, θ=π/4 y N=100, 200, 500 y 1000. También dibujaremos el periodograma y la aproximación de la PSD de estos procesos.

% Construir un proceso AR-2
clear all
ro=0.95;
th=pi/4;
phi1=-2*ro*cos(th);
phi2=ro^2;
N=100;
a=randn(1,N);

H=[1 phi1 phi2];
x=filter(1,H,a);
% Dibujo del filtro
figure(1)
plot(x,'x-')
title('N=100')

% Periodograma
figure(2),
title('N=100')
hold on
PSD1=fftshift(abs(fft(x)));
x1=PSD1.^2/N;
tau1=-1/2:1/length(x):1/2-1/length(x);
plot(tau1,log(x1),'xr-')

% Aproximación de la PSD
PSD=fftshift(abs(fft(H,N))).^2;
PSD=1./PSD;
tau2=-1/2:1/N:1/2-1/N;
plot(tau2,log(PSD))
legend('Periodograma','Aproximación PSD',0)

Procesos autoregresivos para 100, 200, 500 y 1000 valores

Representación de las PSD del proceso x(t) y sus aproximadas

Repetimos el proceso anterior pero variando los valores de ρ y θ pero manteniendo el mismo N

Proceso AR-2 con valores de ρ=0.5, θ=π/4 y N=1000

Proceso AR-2 con valores de ρ=0.5, θ=π/4 y N=1000

Proceso AR-2 con valores de ρ=0.1, θ=0 y N=1000

Proceso AR-2 con valores de ρ=0.1, θ=0 y N=1000