Procesamiento Adaptativo de Señales

Dr. Ing. Marcelo Risk, Ing. Federico Milano

Índice1. Introducción 2

2. Filtrado adaptativo 22.1. Estructura general de los filtros adaptativos . . . . . . . . . . . . . . 22.2. Identificación de los parámetros . . . . . . . . . . . . . . . . . . . . . 32.3. Corrección adaptativa de señales . . . . . . . . . . . . . . . . . . . . . 32.4. Técnicas de estimación óptima . . . . . . . . . . . . . . . . . . . . . . 42.5. Método de resolución directa . . . . . . . . . . . . . . . . . . . . . . . 62.6. Métodos iterativos . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2.6.1. Algoritmo de gradiente mediante el método de Newton . . . . 72.6.2. Método de gradiente buscado por descenso escalonado . . . . . 82.6.3. Algoritmo LMS de Widrow . . . . . . . . . . . . . . . . . . . 8

3. Aplicaciones del filtrado adaptativo 93.1. Modelización adaptativa en la síntesis de filtros FIR . . . . . . . . . . 9

3.1.1. Cancelación adaptativa de ruido . . . . . . . . . . . . . . . . . 93.1.2. Modelización adaptativa en la síntesis de filtros IIR . . . . . . 10

3.2. Implementación de filtros adaptativos . . . . . . . . . . . . . . . . . . 103.2.1. Estructura básica de un filtro adaptativo . . . . . . . . . . . . 10

3.3. Implementación del combinador lineal adaptativo . . . . . . . . . . . 143.3.1. Código fuente de FILTADAPT.H . . . . . . . . . . . . . . . . 163.3.2. Código fuente de FILTADAPT.CPP . . . . . . . . . . . . . . 163.3.3. Código fuente de FiltroAdaptativo.py . . . . . . . . . . . . . . 19

4. Ejemplos de aplicaciones del filtrado adaptativo 204.1. Ejemplo de identificación de un sistema . . . . . . . . . . . . . . . . . 20

4.1.1. Código fuente de FILTRO.CPP . . . . . . . . . . . . . . . . . 214.1.2. Código fuente de IdentificacionSistema.py . . . . . . . . . . . 26

4.2. Ejemplo de cancelación de interferencias . . . . . . . . . . . . . . . . 284.2.1. Código fuente de FILTRO2.CPP . . . . . . . . . . . . . . . . 294.2.2. Código fuente de CancelacionRuido.py . . . . . . . . . . . . . 34

1

1. IntroducciónEn el presente apunte se puede apreciar una breve introducción al procesamiento

adaptativo de señales, y su relación con las redes neuronales (RN), es especial conlas arquitecturas de RN basadas el modelo de neurona denominado perceptrón.

La segunda sección de este documento describe los fundamentos del filtradoadaptativo, y luego una implementación de los mismos a partir de un combinadorlineal adaptativo. Finalmente se estudiarán dos aplicaciones, uno a la identificaciónde un sistema y otro para la cancelación de interferencias.

2. Filtrado adaptativoLos filtros adaptativos son aquellos filtros autoprogramables cuya respuesta en

frecuencia o función transferencia está alterada o adaptada para dejar pasar sindegradación las componentes deseadas y atenuar las señales no deseadas, o quegeneran interferencias, o reducir cualquier distorsión de senal de entrada.

Los primeros filtros adaptativos o de autoaprendizaje fueron descriptos por Lucky[3], quien los utilizó en el diseño de un ecualizador para compensar la distorsión ensistemas de transmisión de datos.

Como generalmente no existe información disponible a priori, los filtros adaptati-vos requieren un período inicial de aprendizaje, también denominado de adaptación.

Los algoritmos utiizados en el procesamiento digital de señales han tenido suprecedente en la teoría de los sistemas adaptativos de control, si bien existen conrespecto de éstos ciertas diferencias.

2.1. Estructura general de los filtros adaptativos

Un filtro adaptativo está compuesto principalmente por tres partes:

a) El índice de performance, el cual debe ser optimizado.

b) El algoritmo que recalcula los parámetros del filtro.

c) La estructura del filtro, que realiza las operaciones requeridas sobre la señal.

El índice de performance depende de la aplicación, es decir, está determinadopor el tipo de trabajo del sistema adaptativo. En general se toma como evaluadorla minimización del error cuadrático medio e el cual generalmente es satisfactoria,como lo muestra la ecuación 9, donde y es la salida del sistema, ymodelo es la salidadel modelo.

e =∑

(y − ymodelo)2 (1)

El algoritmo de optimización es el mecanismo mediante el cual los parámetros,optimizando el criterio, son calculados. Existen dos tipos de algoritmos, el primero

2

es no recursivo, y requiere la colección de todos los datos en una ventana temporal,para luego resolver un sistema de ecuaciones.

El segundo tipo de algoritmo es el algoritmo del gradiente, el cual requiere lasolución de un conjunto de ecuaciones lineales por inversión de una matriz y de estaforma los resultados no estarán disponibles en tiempo real.

En rigor los algoritmos más utilizados desde el punto de vista del procesamientodigital son del tipo recursivo, el cual se actualiza a partir de cada entrada de la señalde entrada o con un pequeño grupo de muestras.

Los resultados están disponibles inmediatamente y es posible el seguimiento deseñales no estacionarias. Un ejemplo de ello es el algoritmo LMS (least mean squa-re). Como es obvio la estructura del filtro depende tanto del algoritmo como de laaplicación.

2.2. Identificación de los parámetros

La identificación de parámetros de un sistema es un procedimiento muy utilizadopara el análisis en sistemas de control. Para controlar el comportamiento óptimo deun sistema es necesario conocer su conducta dinámica.

Esta está usualmente dada por ecuaciones diferenciales que relacionan la entraday la salida del sistema. Si sólo se reconoce la estructura de las ecuaciones, pero nolos parámetros, algún tipo de algoritmo de computación podría ser aplicado paraestimarlos.

La figura 1 muestra el diagrama en bloques de un sistema de identificación adap-tativo.

Las entradas y salidas ruidosas son medidas tanto durante la aplicación normal odurante la prueba de especificaciones. Estas son ingresadas a un filtro con coeficientesvariables, los cuales están ajustados por un algoritmo optimizador.

Luego de la adaptación, el filtro representa el mejor modelo para el sistema, ysus coeficientes son los parámetros identificados del sistema.

2.3. Corrección adaptativa de señales

Si para una entrada del sistema innaccesible se debe actuar sobre la señal desalida de forma tal de realizar algún tipo de corrección de la señal, como por ejemplola eliminación de interferencia en la señal de electrocardiograma (ECG), o para laseparación del ECG materno del fetal, entonces en este caso la información de laseñal y de la corrección requerida es de extrema importancia. Un diagrama de estetipo de filtro correctivo se ve en la figura 2.

Otro tipo de ejemplo muy utilizado de filtro adaptativo que se conoce como filtrode predicción lineal, es aquel que provee a su salida las salidas futuras del proceso.

La misión específica de un filtro de predicción lineal es la de estimar la salidafutura del proceso como una combinación lineal de las salidas presente y pasada.

Originariamente los filtros de predicción se utilizaron para predecir las trayecto-rias de un blanco en movimiento en aplicaciones militares, y hoy en día se utilizan

3

e(t) s(t)

ruido

ruido

Sistema

Modelo

Ajuste del

Modelo

Criterio

de Error

+

+

Figura 1: Diagrama en bloques de un sistema de identificación adaptativo.

para predecir la evolución de procesos muy diversos. La figura 3 muestra el dDia-grama en bloques de un filtro de prediccion lineal.

La señal de entrada s, y versión retardada de la misma son enviadas al procesadoradaptativo, el cual debe intentar predecir la entrada para hacer que con y y d secalcule el error, buscando que este último se minimize. Los filtros de predicción seutilizan además en decodificación de señales y reducción de ruido.

2.4. Técnicas de estimación óptima

Como se dijo anteriormente el error es la resta (según la aplicación) entre lasalida del sistema y la del modelo:

e = y − ymodelo (2)

En el caso de señales discretas este será:

en = yn − wn (3)

si para simplificar la nomenclatura hacemos ymodelo = w:

en = yn − wn (4)

4

FiltroAdaptativo

AlgoritmoAdaptativo

Ruido

Señal corregida

Señal a ser corregida

w

Conocimientoapriori de

señal y ruido

Criterio

Σ

Figura 2: Diagrama en bloques de un sistema de correccion adaptativa de señales.

Delay ProcesadorAdaptativo

s d

error

y+

- Σ

Figura 3: Diagrama en bloques de un filtro de prediccion lineal.

5

si tomamos para identificar un modelo tipo ARMA:

W (z)

X(z)=

∑mj=0 βjz

−j

1−∑p

i=0 αiz−j(5)

de donde:

wn =

p∑i=1

α′

iwn−i +m∑j=0

β′

jxn−j (6)

Entonces para el error entre la salida del sistema físico y modelo se utiliza laecuación 3, luego:

yn =

p∑i=1

αiyn−i +m∑j=0

βjxn−j (7)

en =

p∑i=1

(αi − α′

i)yn−i +m∑j=0

(βj − β′

j)xn−j (8)

El error cuadrático medio (ECM ) es una magnitud se calcual como:

ξ = E(e2n) (9)

La expresión anterior representa en el espacio generado por los coeficientes αi yβj un hiperparaboloide de dimensiones p+m+1, que posee un mínimo absoluto enel punto dado por:

αj = α′

j 1 ≤ i ≤ p

βj = β′

j 0 ≤ j ≤ m(10)

El problema de estimar los coeficientes se reduce a calcular los coeficientes αi yβj que minimizan el error cuadrático ξ (ver ecuación 9).

2.5. Método de resolución directa

Los métodos de resolución directa consisten en calcular el punto en que se anulantodas las derivadas del error cuadrático medio, para las siguientes condiciones:

δξ

δαi

= 0 1 ≤ i ≤ p

δξ

δβj= 0 0 ≤ j ≤ m

(11)

Dado que la función ξ no posee máximos para los valores finitos de αi y βj lasolución de las ecuaciones conducen al mínimo deseado.

6

Si se trata de un proceso cuasiestacionario, debe emplearse la esperanza mate-mática de ECM como la suma temporal dentro del período de estacionareidad delproceso, luego el cálculo debe repetirse periódicamente, de manera de cubrir todo eltiempo en estudio.

2.6. Métodos iterativos

Estos métodos involucran un menor número de operaciones matemáticas y pro-veen además una estimación de α′

i y β′j adaptable a las fluctuaciones de los coe-

ficientes verdaderos αi y βj. De los diversos métodos iterativos sólo estudiaremosel algoritmo del gradiente mediante el método de Newton-Raphson, el método delgradiente buscado por descenso escalonado, y el algoritmo LMS de Widrow.

2.6.1. Algoritmo de gradiente mediante el método de Newton

Este método consiste en descender sobre el hiperparaboloide de ECM siguiendola dirección impuesta por el método de Newton-Raphson para la determinación deuna raíz y por ende con una sencilla transformación, el mínimo de:

f′(x0) =

f(x0)

x0 − x1(12)

entonces:

x1 = x0 −f(x0)

f ′(x0)(13)

luego:

xk+1 = xk − f(xk)

f ′(xk)k = 0, 1, ... (14)

La convergencia depende del valor inicial y de la naturaleza de f(x). Si ahoratomamos:

f(x) = ξ′(x) (15)

en vez de un cero se tendrá un mínimo, entonces:

αk+1i = αk

i − µξ′(x)

ξ′′(x)1 ≤ i ≤ p

βk+1i = βk

i − µξ′(x)

ξ′′(x)0 ≤ j ≤ m

(16)

El coeficiente µ es una constante que determina la dimensión del peso, y gobiernala velocidad de convergencia, es decir un valor de µ demasiado pequeño requiere unelevado número de iteraciones para alcanzar el mínimo.

7

2.6.2. Método de gradiente buscado por descenso escalonado

En este caso los coeficientes o pesos son ajustados según el gradiente en cadapaso, lo cual simplifica notablemente el cálculo.

Este método está gobernado por las ecuaciones en la cual µ es una constanteque regula el tamaño del paso. Al igual que el mtodo anterior, un valor demasiadopequeño de µ requerirá un elevado número de iteraciones para alcanzar un mínimo,y un valor excesivo puede ocasionar inestabilidades en el algoritmo.

Tanto en este caso como en el anterior se parte utilizando todo el conocimientoprevio sobre los valores a estimar y luego es desplazada la estimación inicial descen-diendo sobre la superficie de ECM en la dirección del gradiente (perpendicular a lascurvas de nivel) hasta alcanzar un mínimo.

La diferencia entre este método y el de Newton-Raphson radica en la practicidaddel algoritmo puesto que en el primero es necesario resolver una matriz inversa.

El cálculo de los coeficientes en este método se realiza con las siguientes ecuacio-nes:

αk+1i = αk

i − µδξ

δαi

1 ≤ i ≤ p

βk+1i = βk

i − µδξ

δβj0 ≤ j ≤ m

(17)

Donde el supraíndice indica el orden de iteración.

2.6.3. Algoritmo LMS de Widrow

En el algoritmo anterior se requiere la determinación del gradiente del ECM encada iteración, es decir se estimo el gradiente de ξ.

Bernard Widrow analizó las condiciones de estabilidad, la constante de tiempode convergencia, y el error de estimación del algoritmo anterior, reemplazando elvalor del gradiente por una estimación del mismo empleando una sola muestra [7].

Es decir, aproximando la esperanza matemática del ECM por el mismo pero encada iteración. En otras palabras se toma a e2k como una estimación directa de ξ,entonces:

αk+1i = αk

i − µδ(e2k)

δαi

1 ≤ i ≤ p

βk+1i = βk

i − µδ(e2k)

δβj0 ≤ j ≤ m

(18)

Si tomamos la definición del ECM y elevamos al cuadrado ambos miembros,calculamos la derivadas parciales y reemplazamos en la ecuación anterior se obtiene:

αk+1i = αk

i + 2µekyk−1 (19)

βk+1i = βk

i + 2µekxk−1 (20)

8

3. Aplicaciones del filtrado adaptativo

3.1. Modelización adaptativa en la síntesis de filtros FIR

La idea básica es asociar a un pseudofiltro las especificaciones ideales de un filtroque generalmente no será físicamente realizable.

El pseudofiltro es únicamente conceptual; el esquema utilizado puede verse en lafigura 4.

Pseudofiltro

FiltroAdaptativo

d

error

y+

-

x

f 1

f 2

f n

~~~

ΣΣ

Figura 4: Diagrama en bloques de un sistema de filtrado adaptativo a través de unpseudofiltro.

3.1.1. Cancelación adaptativa de ruido

La situación básica de cancelación de ruido está ilustrada en la siguiente figura5.

Una señal es transmitida de un canal a un sensor que recibe la señal más ruido nocorrelacionado n0. Esta combinación aditiva forma la primera entrada. Un segundosensor recibe ruido n1 el cual no está correlacionado con la señal pero sí lo está dealguna manera con el ruido n0. Esta constituye la entrada de referencia al cancelador.el ruido n1 es filtrado de forma tal de obtener una réplica n0.

La salida se resta de la entrada primaria para producir una salida del sistemaigual a s + n0 + y.

Ejemplo de ello son la cancelación de ruido en señales de audio (n1 es la in-terferencia), cancelación de la interferencia en el ECG del corazón transplantado y

9

FiltroAdaptativo

Fuentede señal

EntradaPrimaria

EntradaReferencia

Fuentede ruido

error

Salida

Cancelador Adaptativo de Ruido

y

-

+

nO

s+nO

n1

Σ

Figura 5: Diagrama en bloques de un sistema de cancelación de ruido con filtradoadaptativo.

separación del ECG materno del fetal como así también la cancelación de ecos encircuitos telefónicos de larga distancia.

3.1.2. Modelización adaptativa en la síntesis de filtros IIR

Como se ha visto anteriormente la transferencia de un filtro de respuesta infinitaal impulso (IIR) a sintetizar será:

H(z) =A(z)

1−B(z)(21)

El diagrama en bloques de tal filtro se puede apreciar en la figura 6.

3.2. Implementación de filtros adaptativos

3.2.1. Estructura básica de un filtro adaptativo

Los filtros adaptativos se pueden implementar tanto como filtros de respuestafinita al impulso (FIR) o de respuesta infinita al impulso (IIR), pero en general seprefiere el esquema FIR debido a una mayor simplicidad y estabilidad, teniendo encuenta que el algoritmo de adaptación es el encargado de calcular los coeficientesdel filtro, por lo cual un esquema FIR asegura la estabilidad del filtro a pesar de lascombinaciones de coeficientes calculadas.

La estructura de un filtro FIR se puede desarrollar con un combinador lineal,cuyo diagrama se muestra en la figura 7.

10

+A(z)

x

y

B(z)

Figura 6: Diagrama en bloques de un filtro adaptativo IIR.

X0

X1

X2

X3

W0

W1

YW2

W3

Σ

Figura 7: Diagrama en bloques de un combinador lineal.

11

La expresión del combinador es:

Y =N−1∑k=0

WkXk (22)

donde N es la cantidad de coeficientes u orden del filtro, X es el vector deentrada, W es el vector de pesos, e Y es la salida.

La figura anterior muestra el combinador lineal en la forma de entrada paralelo,pero en la mayoría de las aplicaciones es necesario ingresar los datos (muestras) enforma serie, especialmente cuando el filtro se implementa en tiempo real y por lotanto las muestras van llegando una a una, en este caso se debe utilizar la formaserie del combinador lineal, como lo muestra la figura 7.

X(n)

X(n-1)

X(n-2)

X(n-3)

W0

W1

Y(n)W2

W3

T

T

T

Σ

Figura 8: Diagrama en bloques de un combinador lineal serie.

En el combinador lineal con entrada serie podemos apreciar cómo se forma elvector de entrada X con la muestra actual X(n) y las muestras anteriores X(n− 1)a X(n − 3), para lograr el retardo entre dichas muestras se intercala un bloque deretardo T . La expresión del combinador lineal con entrada serie es la siguiente:

Y (n) =N−1∑k=0

WkXk(n− k) (23)

Como se describió en la introducción, los sistemas de lazo cerrado utilizan ladiferencia entre la salida deseada d(n) y salida actual Y (n) para calcular en error:

e(n) = d(n)− Y (n) (24)

12

por lo tanto:

e(n) = d(n)−N−1∑k=0

WkXk(n− k) (25)

Elevando el error al cuadrado obtenemos el denominado error cuadrático instan-táneo:

e2(n) = d2(n)− 2N−1∑k=0

WkXk(n− k) +

(N−1∑k=0

WkXk(n− k)

)2

(26)

Por otro lado se define a error cuadrático medio como:

ξ = E[e2(n)

]=

M−1∑m=0

e2(n) (27)

quedando:

ξ =M−1∑m=0

d2(n)− 2N−1∑n=0

Wkrdx(n) +N−1∑n=0

N−1∑l=0

WkWlrxx(n− l) (28)

donde las ecuaciones rdx y rxx las correlaciones cruzadas entre la salida deseaday la entrada y la autocorrelación da la entrada respectivamente.

El error cuadrático medio ξ describe una superficie de performance, cuando selo grafica con respecto a los pesos W , en la figura siguiente podemos observar unasuperficie de performance para el caso de dos pesos, donde también apreciamoslos pesos óptimos W1 y W2 que corresponden al mínimo de ξ de la superficie deperformance.

El proceso de adaptación busca minimizar el error cuadrático medio, para lograreste objetivo se propone el método de descenso a pasos a través de la superficie deperformance.

Con el método de descenso a pasos es posible recalcular cada peso teniendoinformación acerca del peso anterior y del gradiente en el punto de peso anteriorsobre la superficie de performance, multiplicado por un coeficiente de velocidad deconvergencia, como lo indica la ecuación 29.

Wk(n+ 1) = Wk(n) + µ(−5k) (29)

En el método de adaptación basado en el error cuadrático medio, el gradiente sedefine como:

5k =δ[e2(n)

]δWk(n)

(30)

entonces sustituyendo la ecuación 30 en la 29:

Wk(n+ 1) = Wk(n)− µδ[e2(n)

]δWk(n)

(31)

13

Se utiliza el error instantáneo porque está disponible muestra a muestra du-rante la optimización de los pesos del filtro. Entonces desarrolando la derivada yreemplazando por la definición de la ecuación 25, tenemos:

Wk(n+ 1) = Wk(n)− 2µe(n)δ[e(n)

]δWk(n)

(32)

y reemplazando con 25, quedando:

Wk(n+ 1) = Wk(n)− 2µe(n)δ[d(n)−

∑N−1k=0 WkXk(n− k)

]δWk(n)

(33)

Note que d(n) y Xk(n) son independientes de Wk, por lo tanto:

Wk(n+ 1) = Wk(n)− 2µe(n)Xk(n) (34)

La ecuación 34 es la expresión final del algoritmo del método de los cuadradosmínimos.

El coeficiente de velocidad de convergencia µ debe elegirse de forma tal de llegaral valor óptimo lo suficientemente rápido, pero sin perder estabilidad, por lo tantolos valores recomendados para µ deben estar en el rango descripto en la ecuación??:

0 < µ <1

20NPx

(35)

donde N es la cantidad de coeficientes del filtro, y Px es la potencia de los datosdurante el proceso de optimización:

Px =1

M + 1

M−1∑n=0

X2(n) (36)

En la ecuación 36 M es la cantidad de muestras durante la optimización de lospesos.

El combinador lineal cuando utiliza el método de los cuadrados mínimos se puededescribir con el diagrama en bloques de la figura 9, de esta forma se denominacombinador lineal adaptativo.

3.3. Implementación del combinador lineal adaptativo

El combinador lineal adaptativo se puede implementar como una clase, en len-guaje C++ [4], y en Python [5, 2]. Dentro de la clase FiltroAdaptativo están encap-suladas todas las propiedades necesarias para el funcionamiento de un objeto quemaneja un vector de datos Xk, un vector de pesos Wk, y otras variables tales elerror Ek, la cantidad de elementos Cantidad y el valor de coeficiente de velocidad deoptimización Alfa, demás de las funciones necesarias para la carga de los vectores,el cálculo de la salida y la optimización de los coeficientes del filtro.

14

X(n)

X(n-1)

X(n-2)

X(n-3)

W0

W1

d(n)

e(n)

Y(n)W2

W3

T

T

T

Algoritmo decuadrados mínimos

-+

Σ

Σ

Figura 9: Diagrama en bloques de un combinador lineal serie adaptativo.

15

3.3.1. Código fuente de FILTADAPT.H

// FILTADAPT.H

#include <math.h>#include <stdlib.h>#include <iostream.h>#include <fstream.h>

class FiltroAdaptativo {double *Xk;double *Wk;double Ek;double Sk;int Cantidad;

public:double Alfa; // en Alfa a mu por 2

FiltroAdaptativo( int n = 20);~FiltroAdaptativo();

double Calcula( void);void LlenaPesos( double *Array, int cantidad);void LlenaEntradas( double *Array, int cantidad);void LlenaPesos( double *Array);void LlenaEntradas( double *Array);double MuestraSalida( void); void MuestraSet( void);void CalculaNuevoSet( double e);double put ( double d); double operator ()( double d){return put (d);}int GuardaPesos( char *name); int LeePesos( char *name);};

3.3.2. Código fuente de FILTADAPT.CPP

// FILTADAPT.CPP:

#include "filtadapt.h"#include <stdlib.h>

FiltroAdaptativo::FiltroAdaptativo( int n){

16

Alfa = 1.0 / (100.0*double(n));Ek = 1.0;Sk = 0.0;Cantidad = n;Xk = new double[Cantidad];Wk = new double[Cantidad];for( int j = 0; j < Cantidad; j++){Xk[j] = 0.0; Wk[j] = 0.0;}}

FiltroAdaptativo::~FiltroAdaptativo(){delete Xk;delete Wk;}

void FiltroAdaptativo::LlenaPesos( double *Array, int cantidad){for( int j = 0; j < cantidad; j++)Wk[j] = Array[j];}

void FiltroAdaptativo::LlenaPesos( double *Array){for( int j = 0; j < Cantidad; j++)Wk[j] = Array[j];}

void FiltroAdaptativo::CalculaNuevoSet( double e){Ek = e;for( int j = 0; j < Cantidad; j++){Wk[j] = Wk[j] + Alfa*Ek*Xk[j];}}

void FiltroAdaptativo::MuestraSet( void){cout << endl;for( int j = 0; j < Cantidad; j++)cout << Wk[j] << ", ";}

17

double FiltroAdaptativo::Calcula( void){double Suma = 0.0;for( int j= 0; j < Cantidad; j++)Suma += Wk[j] * Xk[j];Sk = Suma;return Sk;}

double FiltroAdaptativo::put( double d){for( int j = Cantidad-1; j > 0; j--)Xk[j] = Xk[j-1];Xk[0] = d;return Calcula(); }

void FiltroAdaptativo::LlenaEntradas( double *Array, int cantidad){for( int j = 0; j < Cantidad; j++)Xk[j] = Array[j];Calcula();}

void FiltroAdaptativo::LlenaEntradas( double *Array){for( int j = 0; j < Cantidad; j++)Xk[j] = Array[j];Calcula();}

double FiltroAdaptativo::MuestraSalida( void){return Sk;}

int FiltroAdaptativo::GuardaPesos( char *name){ofstream A;A.open(name);if( !A)return -1; A << Cantidad << endl;for( int j = 0; j < Cantidad; j++)A << Wk[j] << endl;

18

A.close();return 0;}

int FiltroAdaptativo::LeePesos( char *name){ifstream A;A.open(name);if( !A)return -1;int j;A >> j;if( j != Cantidad)return -2;for( int j = 0; j < Cantidad; j++)A >> Wk[j];A.close();return 0;}

3.3.3. Código fuente de FiltroAdaptativo.py

import sysimport mathfrom random import *from pylab import *

class FiltroAdaptativo:

def __init__(self,n):self.Ek = 1self.Sk = 0self.Cantidad = nself.Alfa = 1.0 / (100.0*float(self.Cantidad)) # en Alfa va mu por 2self.Xk = zeros(self.Cantidad)self.Wk = zeros(self.Cantidad)

def Calcula(self):Suma = 0.0for j in range(self.Cantidad):

Suma += self.Wk[j] * self.Xk[j]self.Sk = Suma

def LlenaPesos(self,Array):

19

for j in range(self.Cantidad):self.Wk[j] = Array[j]

def LlenaEntradas(self,Array):for j in range(self.Cantidad):

self.Xk[j] = Array[j]

def MuestraSalida(self):print "Sk = " + self.Sk

def MuestraSet(self):print self.Wk

def MuestraEntradas(self):print self.Xk

def CalculaNuevoSet(self,ee):self.Ek = eefor j in range(self.Cantidad):

self.Wk[j] = self.Wk[j] + self.Alfa*self.Ek*self.Xk[j]

def put(self,dd):for j in range(1,self.Cantidad):

self.Xk[self.Cantidad-j] = self.Xk[self.Cantidad-1-j]self.Xk[0] = ddself.Calcula()return self.Sk

4. Ejemplos de aplicaciones del filtrado adaptativo

4.1. Ejemplo de identificación de un sistema

Los filtros adaptativos son capaces de identificar un sistema desconocido, en otraspalabras de desarrollar una modelización, el objetivo es encontrar un filtro FIR quemodelice el sistema desconocido [1].

El diagrama en bloques del proceso de modelización se muestra en la figura3b. En los párrafos siguientes se describirá una implementación en C++ para laidentificación de un sistema, en este ejemplo el sistema desconocido es un filtro IIRcon la siguiente ecuación de recurrencia:

y(n) = ax(n) + by(n− 1) (37)

donde a y b son los coeficientes del filtro, cuya condición de estabilidad es que lasuma de los coeficientes sea menor o igual a uno.

20

A continuación se describen los pasos a seguir para lograr dicho proceso:

1) Ambos sistemas, el filtro y el sistema desconocido, deben ser excitados por ungenerador de ruido aleatorio de valor medio nulo.

2) El filtro adaptativo comienza un proceso de optimización, comenzando los pesoscon valores iniciales aleatorios pequeños; utilizando el método de los cuadradosmínimos se optimizan los pesos (coeficientes del filtro) minimizando el error.

3) La cantidad de iteraciones necesarias para la optimización se pueden calcularbuscando un error mínimo y una cantidad máxima de iteraciones, en caso de nolograr el error mínimo.

En el ejemplo presentado el sistema desconocido (filtro IIR) y el filtro tienen unarespuesta al impulso y al escalón como lo muestran las figuras 10 y 11, respectiva-mente, mientras que la figura 12 muestra el error2 para las 1000 primeras iteracionesde la adaptación.

Figura 10: Respuesta del filtro IIR y del filtro adaptativo a un impulso.

4.1.1. Código fuente de FILTRO.CPP

// FILTRO.CPP: Identificacion de un sistema con un filtro

21

Figura 11: Respuesta del filtro IIR y del filtro adaptativo a un escalón.

22

Figura 12: error2 para las 1000 primeras iteraciones de la adaptación.

23

// adaptativo

#include <iostream.h>#include <fstream.h>#include <math.h>#include "filtadapt.h"

void main( void) {int CantidadCoeficientes = 20;

cout << "Orden del filtro adaptativo = ";cin >> CantidadCoeficientes;

FiltroAdaptativo F(CantidadCoeficientes);

int Number = 5000;

ofstream Salida;Salida.open( c:\neural\\filtro.txt);

int j;double Px = 0.0;double Mean = 0.0;double *noise;noise = new double[Number];// Calcula ruido blnco sin valor mediofor( j = 0; j < Number; j++){

noise[j] = double(random(1000))/1000.0;mean += noise[j];}

Mean /= double(Number);for( j = 0; j < Number; j++){

noise[j] -= Mean;Px += (noise[j]*noise[j]);}

// Calcula potencia del ruidoPx /= double(Number+1);

double x, y, d, dAnt, e;// coeficientes del sistema a identificardouble a = 0.3;double b = 0.7;cout << "Coeficientes del sistema a identificar: " << endl;cout << "y[n] = a*x[n] + b*y[n-1]" << endl;cout << "a = ";

24

cin >> a;cout << "b = ";cin >> b;

dAnt = 0.0;

// ajuste valor Alfa (coeficiente de velocidad deconvergencia)// del algoritmo LMSF.Alfa = 1.0 / (20.0*double(CantidadCoeficientes)*Px);cout << "Busca la solucion optima por el metodo LMS ..." << endl;for( j = 0; j < Number; j++){

x = noise[j];d = a*x + b*dAnt;y = F(x);e = d - y;F.CalculaNuevoSet( e);dAnt = d;if( j == 0)

cout << "Error inicial = " << e << endl;}

cout << "Error final = " << e << endl;F.GuardaPesos( "c:\\neural\\pesos.txt");

// Respuesta al impulso de ambos filtrosSalida << "Respuesta al impulso" << endl;cout << "Respuesta al impulso..." << endl;Salida << "x\td\ty" << endl;for( j = 0; j < 50; j++){

x = 0.0;if( j == 25)

x = 1.0;d = a*x + b*dAnt;y = F(x);Salida << x << "\t" << d << "\t" << y << endl;dAnt = d;}

Salida << endl;Salida << "Respuesta al escalon" << endl;cout << "Respuesta al escalon..." << endl;Salida << "x\td\ty" << endl;for( j = 0; j < 50; j++){

x = 0.0;if( j == 25)

25

x = 1.0;d = a*x + b*dAnt;y = F(x);Salida << x << "\t" << d << "\t" << y << endl;dAnt = d;}

delete noise;cout << endl << "FIN";}

4.1.2. Código fuente de IdentificacionSistema.py

import syssys.path.append(’....’) # poner el path de FiltroAdaptativo.py

import mathfrom random import *from pylab import *from FiltroAdaptativo import *

F = FiltroAdaptativo(20)Number = 5000Px = 0.0Mean = 0.0noise = zeros(Number)

seed()for j in range(Number):

noise[j] = normalvariate(0.0,1.0)Px += noise[j] * noise[j]

Px /= float(Number+1)

x = 0.0y = 0.0dAnt = 0.0e = 0.0a = 0.3b = 0.7

F.MuestraSet()

F.Alfa = 1.0 / (20.0*float(F.Cantidad)*Px)

26

print F.Alfa

for j in range(Number):x = noise[j]d = a*x + b*dAnty = F.put(x)e = d - yF.CalculaNuevoSet(e)dAnt = d

F.MuestraSet()

N = 100xx = arange(N)

yy = zeros(N)zz = zeros(N)

yy[N/2] = 1.0

F.LlenaEntradas(zeros(20))

for j in range(1,N):zz[j] = F.put(yy[j])yy[j] = a*yy[j] + b*yy[j-1]

subplot(211)plot(xx,yy,’b’)ylim(-0.1,1.1)title(’Filtrado con $y(n)=ax(n)+by(n-1)$’)

subplot(212)plot(xx,zz,’g’)ylim(-0.1,1.1)title(’Filtrado con $FiltroAdaptativo$’)

savefig("EjercicioIdentificacionSistemaRespImpulso.eps")close()

yy = zeros(N)zz = zeros(N)

for j in range(N/2,N):yy[j] = 1.0

27

F.LlenaEntradas(zeros(20))

for j in range(1,N):zz[j] = F.put(yy[j])yy[j] = a*yy[j] + b*yy[j-1]

subplot(211)plot(xx,yy,’b’)ylim(-0.1,1.1)title(’Filtrado con $y(n)=ax(n)+by(n-1)$’)

subplot(212)plot(xx,zz,’g’)ylim(-0.1,1.1)title(’Filtrado con $FiltroAdaptativo$’)

savefig("EjercicioIdentificacionSistemaRespEscalon.eps")close()

4.2. Ejemplo de cancelación de interferencias

En muchas ocasiones una señal de gran ancho de banda w(n) se encuentra inter-ferida por la suma de ruido de ancho de banda estrecho s(n), conformando la señalsuma x(n); las dos señales que forman x(n) no están correlacionadas:

x(n) = w(n) + s(n) (38)

La solución a este problema es diseñar un filtro de ancho de banda estrecho(notch), pero en ciertas ocasiones la banda de ruido no es conocida o varía lentamenteen el tiempo, para lo cual los filtros adaptativos pueden utiizarse, pues en cualquierade las dos situaciones son capaces de optimizar su funcionamiento [6].

El filtro adaptativo debe ser capaz de sustraer s(n) de x(n), teniendo en cuentaque debe predecir la componente en cada muestra a través de las muestras anteriores,dicha predicción se logra gracias al ancho de banda estrecho de s(n), y debido a quecomo el ancho de banda de s(n) es mucho menor que el de w(n), las muestras des(n) tienen tienen una alta correlación; por otro lado las muestras de w(n) tienenuna muy baja correlación.

En la figura 13 podemos apreciar el diagrama en bloques del sistema de can-celación de interferencias, donde el decorrelacionador introduce un retardo de Tmuestras, con el propósito de decorrelacionar las componentes s(n) de la señal x(n).

28

retardoT

FiltroAdaptativo

x(n) salida

decorrelacionador e(n) = x(n) - s’(n)

s’(n)

+

-

Σ

Figura 13: Diagrama en bloques de un filtro adaptativo para cancelación de interfe-rencias.

La salida del filtro adaptativo tiene la siguiente expresión:

s′(n) =

N−1∑k=0

Wk(n)x(n− k − T ) (39)

La señal de error e(n) se utiliza en el proceso de optimización del filtro adaptativo,con el método de los cuadrados mínimos. En dicho proceso de optimización x(n) esla suma de una secuencia aleatoria de valor medio nulo, sobre la cual se suma elruido de ancho de banda estrecha, por ejemplo un tono senoidal, o bien una suma dedos o tres tonos puros. El retardo del decorrelacionador en la mayoría de los casoses suficiente con una muestra de retardo.

Una vez adaptado el filtro, la señal a filtrar debe ser ingresada directamenteal filtro; eventualmente el filtro, puede reoptimizarse para seguir las variaciones deseñal interferente.

El ejemplo presentado trata de suprimir una interferencia producida por un tonopuro sobre una señal de electrocardiograma (ECG). El programa perteneciente a esteejemplo es FILTRO2.CPP. y la figura iguiente muestra los resultados del filtradosobre una porción de la señal de ECG.

La figura 14 muestra el error2 para las primeras 500 iteraciones durante la adap-tación del filtro, mientras que la figura 15 muestra en la señal superior el ECGoriginal sin ruido ni interferencias, en la señal central el mismo ECG pero con lainterferencia sumada, y finalmente en la señal inferior la salida del filtro adaptativo,note en esta última el comportamiento para las primeras 20 muestras, en las cualesel filtro adaptativo todavia no se enganchó.

4.2.1. Código fuente de FILTRO2.CPP

// FILTRO2.CPP: cancelacion de interferencias con un filtro adaptativo.

#include <iostream.h>

29

Figura 14: error2 para las primeras 500 iteraciones durante la adaptación del filtro.

30

Figura 15: La señal superior el ECG original sin ruido ni interferencias, en la señalcentral el mismo ECG pero con la interferencia sumada, y finalmente en la señalinferior la salida del filtro adaptativo, note en esta última el comportamiento paralas primeras 20 muestras, en las cuales el filtro adaptativo todavia no se enganchó.

31

#include <fstream.h>#include <math.h>#include "filtadapt.h"

void main( void) {

int CantidadCoeficientes = 20;

cout << "Orden del filtro adaptativo = ";cin >> CantidadCoeficientes;FiltroAdaptativo F(CantidadCoeficientes);

int Number = 4000;

ofstream Salida;Salida.open("filtro2.txt");

int j;double Px = 0.0;double Mean = 0.0;noise = new double[Number];double *ecg;ecg = new double[Number];double *ecgl;ecg1 = new double[Number];

double Pi = 3.14;// Calcula ruido blanco sin valor medio mas ruido// de senoidalfor( j = 0; j < Number; j++){

noise[j] = double(random(1000))/10000.0;noise[j] += sin( 1000.0*2.0*Pi*double(j)/double(number));

Mean += noise[j];}

Mean /= double(Number);for( j = 0; j < Number; j++){

noise[j] -= Mean;Px += noise[j]*noise[j]);}

// Calcula potencia del ruidoPx /= double(Number+1);

double x, y, d, dAnt, dAnt2, e;

char NameInput[80] = "apb.txt";

32

ifstream In;In.open( NameInput);if( !In){

cout << "I can’t open the file !!!!!";return}

double data;for( j = 0; j < Number; j++){

In >> data;data /= 10.0;

ecg[j] = data;ecg[j] += sin(

1000.0*2.0*Pi*double(j)/double(Number));ecg1[j] = data;In >> data;}

// ajuste valor Alfa (coeficiente de velocidad deconvergencia)// del algoritmo LMSF.Alfa = 1.0 / (20.0*double(CantidadCoeficientes)*Px);cout << "Alfa = " << F.Alfa << endl;// Busca la solucion optima por el metodo LMScout << "Busca la solucion optima por el metodo LMS ..." << endl;dAnt = 0.0;dAnt2 = 0.0;for( j = 0; j < Number; j++){

x = noise[j];y = F(dAnt2);dAnt = x;dAnt2 = dAnt;e = x - y;// cout << e << endl;F.CalculaNuevoSet( e);if( j == 0)

cout << "Error inicial = " << e << endl;}

cout << "Error final = " << e << endl;

Salida << endl;Salida << "Filtrado de ruido" << endl;cout << "Filtrado de ruido..." << endl;Salida << "ecg\tx\ty" << endl;dAnt = 0.0;

33

dAnt2 = 0.0;for( j = 0; j < Number; j++){

x = ecg[j];y = F(dAnt2);dAnt = x;dAnt2 = dAnt;d = x - y;Salida << ecgl[j] << "\t" << x << "\t" << d << endl;dAnt = d;}

delete noise;delete ecg;delete ecgl;cout << endl << "FIN";}

4.2.2. Código fuente de CancelacionRuido.py

import syssys.path.append(’...’)

import mathfrom random import *from pylab import *from FiltroAdaptativo import *

F = FiltroAdaptativo(20)

Number = 4000ecg = zeros(Number)ecgOrig = zeros(Number)ecg1 = zeros(Number)noise = zeros(Number)Px = 0.0Mean = 0.0

# calcula ruido blanco + senoidalfor j in range(Number):

noise[j] = normalvariate(0.0,1.0)/10noise[j] += math.sin(1000.0*2.0*math.pi*float(j)/float(Number))Mean += noise[j]

Mean /= float(Number)

34

for j in range(Number):noise -= MeanPx += noise[j]*noise[j]

Px /= float(Number+1)

fh = open(’apb.txt’)for j in range(Number):

linea = fh.readline()w = linea.split()ecgOrig[j] = (float(w[0])-128.0)/10.0ecg[j] = (float(w[0])-128.0)/10.0 +math.sin(1000.0*2.0*math.pi*float(j)/float(Number))

ecg1[j] = (float(w[1])-128.0)/10.0

x = zeros(Number)y = zeros(Number)d = zeros(Number)dAnt = 0.0dAnt2 = 0.0e = zeros(Number)

F.Alfa = 1.0 / (20.0*float(F.Cantidad)*Px)

# adaptacion del filtro adaptativofor j in range(Number):

x[j] = noise[j]y[j] = F.put(dAnt2)dAnt = x[j]dAnt2 = dAnte[j] = x[j] - y[j]F.CalculaNuevoSet(e[j])

# filtrado del ECGdAnt = 0.0dAnt2 = 0.0for j in range(Number):

x[j] = ecg[j]y[j] = F.put(dAnt2)dAnt = x[j]dAnt2 = dAntd[j] = x[j] - y[j]dAnt = d[j];

plot(ecgOrig[:250]+20,’g’)

35

plot(x[:250]+10,’b’)plot(d[:250],’r’)savefig("EjercicioCancelacionRuidoECG.eps")close()

e2 = e*eplot(e2[:500])xlabel(’iteraciones’)ylabel(’$error^2$’)savefig("ErrorCancelacionRuidoECG.eps")close()

Referencias[1] Ricardo Armentano, Javier Fochesatto, and Marcelo Risk. Análisis de Señales

y Sistemas. Editorial Rocamora, 1996.

[2] Perry Greenfield and Robert Jedrzejewski. Using Python for interactive dataanalysis. Association of Universities for Research in Astronomy, 2007.

[3] R. W. Lucky. Automatic equalization for digital communication. Bell SystemsTechnology Journal, 44:547–588, 1965.

[4] Al Stevens and Clayton Walnut. Standard C++ Bible. IDG Books, 2000.

[5] Guido van Rossum and Fred L. Drake. Python tutorial, Release 2.5. PythonSoftware Foundation, docs@python.org, 2006.

[6] Bernard Widrow, John Glover, John MacCool, John Kaunitz, Charles Williams,Robert Hearn, James Zeidler, Eugene Dong, and Robert Goodlin. Adaptive noisecancelling: principles and applications. Proceedings of the IEEE, 63(12):1692–1716, 1975.

[7] Bernard Widrow and Michael Lehr. 30 years of adaptive neural networks: percep-tron, madaline, and backpropagation. Proceedings of the IEEE, 78(9):1415–1442,1990.

36