Tratamiento de señales digitales mediante wavelets y su uso con MATLAB - Félix Martínez Giménez

Tratamiento de Señales Digitales Mediante Wavelets y su Uso con MATLAB

Félix Martínez Martínez Giménez Alfredo Alfredo Peris Peris Manguill Manguillot ot Francisco Francisco Rodenas Rodenas Escribá Escribá

Valencia 2004

Tratamiento de señales digitales mediante wavelets y su uso con Matlab © Félix Martínez Giménez Alfredo Peris Manguillot Francisco Rodenas Escribá ISBN: 978–84–9948–628–4 e-book v.1.0

ISBN edición en Papel: 978-84-8454-387-9

Edita: Editorial Club Universitario. Telf.: 96 567 61 33 C/. Cottolengo, 25 – San Vicente (Alicante) www.ecu.fm Maqueta y diseño: Gamma. Telf.: 965 67 19 87 C/. Cottolengo, 25 – San Vicente (Alicante) www.gamma.fm [email protected]

Reservados todos los derechos. Ni la totalidad ni parte de este libro puede reproducirse o transmitirse por ningún procedimiento electrónico o mecánic mecánico, o, incluyendo fotocopia, grabación magnética o cualquier almacenamiento de información o siste ma de reproducción, sin permiso previo y por escrito de los titulares del Copyright.

Prólogo Las wavelets proporcionan un conjunto de herramientas flexible para problemas prácticos en ciencia e ingeniería. En la última década se han aplicado con éxito al análisis de señales en disciplinas tan diversas como la medicina, la ingeniería eléctrica, teledetección y muchas otras. Una de las principales virtudes de las wavelets es que permiten modelar mejor procesos que dependen fuertemente del tiempo y para los cuales su comportamiento no tiene porqué ser suave. La transformada wavelet resulta especialmente eficiente para extraer información de señales no periódicas o de vida finita. Otra de las ventajas de dicha transformada frente a otros métodos es el de poder disponer de una amplia familia de wavelets, lo cual permite tratar señales de diversa índole. La elección de la wavelet dependerá del tipo de señal que analicemos. Algunos Algunos de los principales principales problemas problemas que afectan afectan al tratamien tratamiento to de señales e imágenes digitales, y en los que las wavelets constituyen una potente herramienta para afrontarlos, son la reducción del ruido (en señales de audio y en imágenes), la compresión de señales (de vital importancia tanto en la transmisión de grandes cantidades de datos como en su almacenamiento) o la detección de determinados objetos en imágenes o irregularidades locales en ciertos tipos de señales (electrocardiogramas, vibraciones de motores, etc.). Esta moderna teoría ha experimentado un gran desarrollo en las dos últimas décadas mostrándose muy eficiente donde otras técnicas, como por ejemplo la transformada rápida de Fourier, no resultaban satisfactorias. En esta última se maneja una base de funciones bien localizada en frecuencia pero no en tiempo, mientras que la mayoría de las wavelets interesantes presentan una buena localización en tiempo y en frecuencia, disponiendo incluso de bases de wavelets con soporte compacto. Este texto se centra fundamentalmente en la transformada wavelet discreta. Esta transformada está íntimamente ligada al análisis de multirresolución, iii

iv

formulado por Mallat a finales de los 80, y facilita su computación rápida cuando la familia de wavelets es ortogonal. La transformada wavelet discreta es una transformación de la señal que la divide en dos tipos de subseñales, la tendencia y las fluctuacio fluctuaciones. nes. La tendencia viene viene a ser una copia de la señal a menor resolución y las fluctuaciones fluctuaciones almacenan almacenan información información referida referida a los cambios locales en la señal original. La tendencia y las fluctuaciones más significativas permiten una compresión de la señal a cambio de descartar información irrelevante y de la eliminación del ruido producido por los aparatos y las condiciones de medida. Según el tipo de medición realizada el ruido correspondiente se comporta matemáticamente siguiendo distribuciones de probabilidad gaussianas, uniformes... El estudio de las fluctuaciones permite detectar anomalías o disfunciones en el comportamiento esperado de la señal inicial. También permite la comparación con patrones para detectar formas en una imagen o una señal unidimensional de forma automática. El nuevo formato de JPEG2000 basa la compresión de imágenes en la transformada wavelet. La mayor parte de las familias de wavelets que utilizaremos son ortogonales, los cual nos permite una transformada inversa de fácil computación, y tan rápida como la transformada directa. No existe una transformada wavelet única, ni que resuelva todos los problemas, a partir del modelado del proceso y de un análisis a priori del tipo de señal tratada y del objetivo que se pretenda (compresión, eliminación del ruido, segmentación de la imagen,...) se busca la familia de wavelets (Haar, Daubechies, Coiflets,...) que mejor coincida con las características de la señal a estudiar. El tratamiento con wavelets discretas permite su aplicación directa a procesos computacionales. Las wavelets continuas presentan por una parte la dificultad de su manejo al tener que evaluar un gran número de integrales y tener en consecuencia una redundancia de información, pero por otra parte permiten la flexibilidad de poder adaptarse a situaciones en las que las discretas no dan un resultado satisfactorio. Concretamente en el tratamiento de imágenes digitales se precisa precisa realzar realzar detalles detalles y detectar detectar texturas analizando analizando la imagen desde distintos ángulos, lo cual es posible con las wavelets continuas al disponer de un mayor número de parámetros que posteriormente se pueden discretizar para su tratamiento computacional. Además la transformada wavelet continua, al proveer una alta resolución temporal y espectral en el espacio transformado, permite observar aspectos sutiles no estacionarios que pueden ser obviados en la transformada discreta, lo cual compensa el mayor coste computacional.

v

Este texto es el resultado de la impartición de la asignatura “Tratamiento de señales e imágenes digitales mediante wavelets”, durante los cursos 2001/2002 y 2002/2003, en el marco del Centro de Formación de Postgrado de la Universitat Politècnica de València, València, como Curso de Formación Formación Específica que versaría sobre el uso de las wavelets para el tratamiento de señales e imágenes digitales. El público público al que iba dirigido dirigido estaba fundamen fundamentalme talmente nte formado formado por ingenieros ingenieros y estudiantes de tercer ciclo de la Universitat Politècnica de València. Este curso pretende ser una introducción al tratamiento de señales mediante wavelets. Los prerrequisitos necesarios son poco exigentes y se reducen a conocimientos básicos de algebra lineal y cálculo. El objetivo del curso es claramente práctico y por esa razón se incluyen en el texto (en los apéndices finales) ejercicios y aplicaciones. Para la realización de los ejercicios prácticos se usó el software comercial “Wavelet Toolbox for MATLAB”. Las rutinas de este paquete se ejecutan en el entorno matemático de computación MATLAB y por lo tanto también necesita estar instalado en nuestro sistema. En la red se pueden encontrar otros paquetes alternativos de libre distribución (como SCILAB) que el lector puede utilizar, bajo los cuales se podrán realizar los ejercicios propuestos si lo desea pero será necesario traducir las órdenes y rutinas descritas al nuevo entorno que se esté utilizando. El texto se estructura en capítulos y apéndices. Los capítulos están dedicados a la exposición teórica de los conceptos y los apéndices son fundamentalmente prácticos. El capítulo 2 se dedica a la explicación básica de la teoría de wavelets unidimensionales tomando el caso concreto de la wavelet de Haar. Muchas de las ideas desarrolladas en este capítulo para este tipo de wavelets serán utilizadas en capítulos posteriores. También se introducen algunas nociones básicas sobre compresión de señales. En el capítulo 3 se introducen otras wavelets ortogonales como las de Daubechies y las Coiflets, estudiando sus propiedades fundamentales. En el capítulo 4 se profundiza en las técnicas de compresión de señales y de eliminación de ruido blanco. Se abordan también cuestiones importantes a tener en cuenta en compresión como la cuantización y la entropía. Los dos últimos capítulos están dedicados a wavelet packets y a la transformada wavelet continua. Como se ha indicado anteriormente el objetivo de la segunda parte del texto es la visualización práctica de las técnicas basadas en wavelets para el tratamiento de señales. Estas prácticas se realizan con el paquete “Wavelet

vi

Toolbox for MATLAB” que se ejecuta dentro del entorno MATLAB. Por esta razón el apéndice A se dedica a dar unos conceptos básicos sobre el uso del programa MATLAB y de la “Wavelet Toolbox for MATLAB”. El apéndice B contiene la mayoría de los ejercicios que se proponen para su resolución de forma paralela al estudio del texto. En algunos se incluye el código MATLAB para su resolución e incluso se aportan imágenes con los resultados. El último apéndice C pretende ofrecer al lector un pequeño repertorio de ejemplos y aplicaciones de tratamiento de señales con wavelets. La mayoría de estas aplicaciones han sido tomadas de el paquete “Wavelet Toolbox for MATLAB” y se utilizan señales que se distribuyen con el citado paquete. Los autore autoress agrade agradecen cen la financi financiaci ación ón del Proy Proyect ectoo Inter Interdis discip ciplin linar ar “La Transformada Wav Wavelet elet en el Tratamiento Tratamiento de Señales e Imágenes Digitales” dentro del programa INNOVA, por el Vicerrectorada de Investigación, Desarrollo e Innovación de la Universitat Politècnica de València (c.e. 20020629).

Valencia, julio de 2004

Los autores

Índice Abreviado Prólogo

iii

Índice Abreviado

vii

Índice General

ix

Índice de figuras

xiii

Índice de códigos para MATLAB

xv

1 Introducción a las Wavelets y prerrequisitos

1

2 Wavelets de Haar

7

3 Familias de wavelets ortogonales

23

4 Compresión de señales y reducción de ruido

41

5 Wavelet packets

51

6 La transformada wavelet continua

55

A MATLAB y la “Wavelet Toolbox”

63

B Ejercicios

79

C Ejemplos y aplicaciones

12 5

Bibliografía

13 7 vii

Índice general Prólogo

iii

Índice Abreviado

vii

Índice general

ix

Índice de figuras

xiii

Índice de códigos para MATLAB

xv

1 Introdu Introducci cción ón a las las Wa Wavelets velets y prer prerreq requisi uisitos tos de álgeb álgebra ra lineal lineal

1

2 Wavelets de Haar

7

2.1 2.2 2.3

Scaling y wavelets . . . . . . . . . . . . . . . . . . . . . . . . . Análisis de multirresolución (MRA) . . . . . . . . . . . . . . . . Transformada de Haar . . . . . . . . . . . . . . . . . . . . . . .

3 Familias de wavelets ortogonales 3.1

3.2 3.3

Wavelets de Daubechies . . . . . . . . . . . . . . . 3.1. 3.1.11 Daube Daubech chie iess db2 . . . . . . . . . . . . . . . . Conceptos básicos . . . . . . . . . . . . . . Cara Caract cter erís ísti tica ca de las las fluct fluctua uaci cion ones es peque pequeña ñass 3.1. 3.1.22 Daube Daubech chie iess db3 . . . . . . . . . . . . . . . . Conceptos básicos . . . . . . . . . . . . . . Cara Caract cter erís ísti tica ca de las las fluct fluctua uaci cion ones es peque pequeña ñass 3.1. 3.1.33 Daube Daubech chie iess dbJ . . . . . . . . . . . . . . . . Coiflets . . . . . . . . . . . . . . . . . . . . . . . . Transformada de Fourier discre creta y MRA . . . . . ix

7 12 16

23 . . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

23 23 23 27 29 29 31 31 32 34

x

4 Compresión de señales y reducción de ruido 4.1 4.2 4.3

41

Cuantización y umbralizado . . . . . . . . . . . . . . . . . . . . Codificación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Medidas de la aproximación . . . . . . . . . . . . . . . . . . . .

42 46 49

5 Wavelet packets

51

6 La transformada wavelet continua

55

A MATLAB y la “Wavelet Toolbox”

63

A.1 MATLAB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.1.1 Breve (muy breve) introdu oducción ión . . . . . . . . . . . . . . A.1.2 Representación de señales . . . . . . . . . . . . . . . . . Datos en un fichero ASCII . . . . . . . . . . . . . . . . . Datos creados a partir de una función . . . . . . . . . . A.2 La “Wavelet Toolbox” . . . . . . . . . . . . . . . . . . . . . . . A.2.1 Wavelets avelets 1-D 1-D desde desde la línea línea de comandos comandos (primer (primeros os pasos) pasos) A.2.2 A.2.2 Wavelet aveletss 1-D 1-D desde desde el entor entorno no gráfico gráfico (prime (primeros ros pasos) pasos) . A.3 Prog ogrram amaando ndo algu alguna nass fun funcion ciones es en MATLA LAB B . . . . . . . . . .

B Ejercicios B.1 B.2 B.3 B.3 B.4 B.5 B.6

Transformada de Haar de una señal . . . . . . . . . . . . . . . . Compresión ión de señales (prime imeros pasos) . . . . . . . . . . . . . Otra Otrass wavelet eletss orto ortogo gona nale les: s: Daub Daubec echi hies es y Coifl Coiflet etss . . . . . . . . Dibujando wavelets and scaling . . . . . . . . . . . . . . . . . . Compresión ión de señales y cuantización . . . . . . . . . . . . . . . Reducción del ruido en una señal . . . . . . . . . . . . . . . . . B.6.1 Comportamie Comportamient nto o del ruido ruido blanco blanco a trav través és de una transtransformada wavelet . . . . . . . . . . . . . . . . . . . . . . B.6.2 Reducción ión del ruido en una señal . . . . . . . . . . . . . B.7 Transformada de Fourier ier Disc iscreta . . . . . . . . . . . . . . . . . B.7.1 “Fast Fourier Transform” . . . . . . . . . . . . . . . . . . B.7. B.7.22 El esp espec ectr troo de de señ señaales les wa wavelet eletss y scal scalin ingg . . . . . . . . . B.7. B.7.33 Anál Anális isis is de frec frecue uenc ncia iass y tran transf sfor orma mada da wavelet elet . . . . . B.8 Detección . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.9 Wavelets packets . . . . . . . . . . . . . . . . . . . . . . . . . . B.9.1 Analizando una señal . . . . . . . . . . . . . . . . . . . . B.9.2 Comprimiendo la señal . . . . . . . . . . . . . . . . . . .

63 63 64 65 66 68 68 74 75

79 79 80 85 88 92 98 98 100 10 5 105 108 1100 11 112 120 120 120

xi

B.9.3 Reducción de ruido . . . . . . . . . . . . . . . . . . . . . 121 B.10 Transformada wavelet continua . . . . . . . . . . . . . . . . . . 121 B.10.1 Transformada wavelet continua desde la línea de comandos121 B.10.2 B.10.2 Transfo ransforma rmada da wav wavele elett contin continua ua en el entor entorno no gráfic gráficoo . . 123

C Ejemplos y aplicaciones C.1 Detección de “breakdown poi points” . . . . C.1.1 Un cambio en la frecuencia . . . C.1.2 Evolución a largo tiempo . . . . C.1.3 Discontinuida idad en derivadas . . . C.2 Identificando frecuencias puras . . . . . C.3 Análisis de un caso real . . . . . . . . . C.4 C.4 Dete Detecc cció ión n de auto auto-s -sem emej ejan anza za (est (estru ruct ctur uraa

Bibliografía

12 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . frac fracta tal) l)

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

125 1 25 1 28 1 29 13 1 13 3 1344 13

13 7

Índice de figuras 1.1

Análisis de multirr irresolución de una señal . . . . . . . . . . . . .

2.1 2.2 2.2 2.3 2.3 2.4 2.4

Tendenc encia y fluctuación de una señal a nivel ivel 1 . . . . . Tend endenci enciaa y fluct uctuaci uación ón de una seña señall a nivel iveles es 1 y 2 . Anális álisis is de mult ultirre irreso solu luci ción ón de una una señ señal a niv nivel 2 . . Trans ransfo form rmad adaa de de Haa Haarr y perfil perfiles es de ener energí gíaa acum acumul ulad adaa

. . . .

. . . .

. . . .

. . . .

. . . .

9 11 14 19

3.1 3.2 3.3 3.4 3.5 3.6 3.6

Vectore ectoress wave wavelet letss y scaling scaling para para db2 . . . . . . Vectore ectoress wave wavelet letss y scaling scaling para para coif . . . . . . Análisis de frecuencias con DFT . . . . . . . . . Espectros de wavelets y scaling . . . . . . . . . Análisis de frecuencias ias con DFT del MRA . . . Espec spectr tros os de wavelet eletss y scal scalin ingg a vario arioss niv niveles eles

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

25 34 36 37 38 39

4.1 4.2 4.2 4.3

Eliminación de ruido blanco . . . . . . . . . . . . . . . . . . . . Trans ransfo form rmad adaa wa wavelet elet de un ruid ruidoo bla blanc ncoo ggau auss ssia iano no . . . . . . . Cuantización de una señal de audio . . . . . . . . . . . . . . . .

44 45 49

6.1 6.2 6.3 6.4 6.5

Wavelet madre “sombrero mejicano” Transformada wavelet continua . . . Espectros sombrero mejicano . . . . Wavelet madre Morlet . . . . . . . . Escalograma . . . . . . . . . . . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

57 58 59 60 61

A.1 A.1 A.2 A.3 A.4

Un ejem ejempl ploo de de plot() . . . . . . . . . . . Creando una señal . . . . . . . . . . . . . Tendencia y fluctuación de una señal . . . Señales promedi edio y detalle a prime imer nivel

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

66 67 69 70

xiii

. . . . .

. . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

3

xiv

A.5 Análisis de tres niveles de una señal . . . . . . . . . . . . . . . A.6 Herra erram mien ienta gráfi gráfica ca de la “Wa “Wavelet elet Too oolb lboox” . . . . . . . . . . .

73 76

B.1 B.1 Transf ansfor orm mada ada de Haa aarr de 1 nivel ivel de una una señ señal . . . . . . . . . . B.2 Transfo ransforma rmada da de Haar de nivel nivel 2 de una una señal y sus perfiles perfiles de energía . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.3 Compresión de una señal . . . . . . . . . . . . . . . . . . . . . . B.4 Comparación Daubec bechies y Coiflets . . . . . . . . . . . . . . . . B.5 Transformadas inversas de Haar . . . . . . . . . . . . . . . . . . B.6 Compresión ión de una señal con cuantización . . . . . . . . . . . . B.7 Transfo ransforma rmada da wav wavele elett de Haar Haar de un un ruido ruido blanc blancoo gaussi gaussiano ano . . B.8 Reducci cción del ruido blan lanco en una señal . . . . . . . . . . . . . B.9 Reducc Reducción ión del del ruido blanc blanco o en una señal señal con la orden wden . . . B.10 B.10 Trans ransfo form rmad adaa dis discr cret etaa de de Fou Fouri rier er de una una seña señall . . . . . . . . . . B.11 B.11 Anál Anális isis is de frec frecue uenc ncia iass de de una una seña señall con con DFT DFT . . . . . . . . . . B.12 Espec pectro de wavelets y scaling . . . . . . . . . . . . . . . . . . . B.13 B.13 Aná Análi lisi siss de frec frecue uenc ncia iass combi combina nand ndoo DFT DFT y wa wavelets elets . . . . . . . B.14 B.14 Dete Detecc cció ión n de de señ señal ales es co cort rtas as en seña señale less gra grand ndes es . . . . . . . . . . B.15 Detección de señales cortas en señales grandes (manipuladas con wavelets) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.16 B.16 Trans ransfo form rmad adaa wa wavelet elet co con ntin tinua de una una señ señal al . . . . . . . . . . .

80

C. 1 C. 2 C. 3 C. 4

Loca ocalización de “breakdown poi points” . . Evoluc lución a largo tiempo de una señal Disc iscontinuidad en la segunda deriva ivada Identificando frecuencias puras . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

81 83 88 90 95 100 103 1 04 106 1077 10 108 1111 11 1155 11 118 1233 12 1 26 12 8 130 13 2

Índice de códigos para MATLAB Archivo ana_nivel_3.m Anális álisis is wavelet elet de una señ señal hasta asta niv nivel 3 . . . . . . . . Archivo energy.m Energía de una señal . . . . . . . . . . . . . . . . . . . . Archivo cumenergy.m Energí ergíaa acu acumula ulada (nor (norma mali liza zad da) de una una señ señal . . . . . Archivo rms.m Error RMS entre dos señales . . . . . . . . . . . . . . . . Archivo compress.m Compresión de una señal . . . . . . . . . . . . . . . . . . Archivo daub_coif.m Mejoría Coiflets respec pecto Daubec bechies . . . . . . . . . . . Archivo wav_sca.m Construcción de wavelets y scaling . . . . . . . . . . . . Archivo compress2.m Compresión señal con cuantización . . . . . . . . . . . . Archivo comp_audio.m Comp Compre resió sión n seña señall audi audioo con con cuan cuanti tiza zaci ción ón y codifi codifica caci ción ón . Archivo white_noise.m Transformada wavelet de ruido blanco . . . . . . . . . . Archivo make_noisy_signal.m Adici dición ón de ruid ruidoo bla blanc ncoo ggau auss ssia iano no a una una seña señall . . . . . . Archivo denoise.m Reducción de ruido . . . . . . . . . . . . . . . . . . . . . xv

. . . .

71

. . . .

77

. . . .

77

. . . .

78

. . . .

81

. . . .

86

. . . .

89

. . . .

93

. . . .

96

. . . .

99

. . . . 10 1011 . . . . 1 01

xvi

Archivo denoise_auto.m Reducción de ruido (con comando wden) . . . . . . . . . . . . Archivo fourier_1.m Transformada discreta de Fourier . . . . . . . . . . . . . . . . Archivo fourier_2.m Análisis Análisis de frecuencias frecuencias con Transf Transforma ormada da discreta discreta de Fourie Fourierr . Archivo spectrum.m Espec pectros de scaling y wavelets . . . . . . . . . . . . . . . . . Archivo fourier_wav.m Anál Anális isis is de frec frecue uenc ncia iass con con wavelet eletss y Fouri ourier er . . . . . . . . . Archivo corr.m Correlación entre dos señales . . . . . . . . . . . . . . . . . . Archivo detect.m Dete Detecc cció ión n de de pat patro rone ness cor corto toss en en señ señal ales es gran grande dess . . . . . . . . Archivo corrfft.m Correlació Correlación n entre entre dos dos señales señales usando usando transfor transformada mada de de Fouri Fourier er Archivo detect_wav.m Dete Detecc cció ión n de pat patro rone ness (inc (incor orpor poraa técn técnic icas as wa wavelet elets) s) . . . . . . Archivo break_anal.m Búsque Búsqueda da de discon discontin tinuid uidade adess (‘‘b (‘‘brea reakdo kdown wn points’ points’’) ’) . . . . . Archivo evol_anal.m Evolución en el tiempo de una señal . . . . . . . . . . . . . . Archivo disc_anal.m Discon continu inuidad en las derivadas . . . . . . . . . . . . . . . . .

. 1 03 . 105 . 107 . 1 09 . 11 1100 . 114 . 11 1144 . 117 . 118 . 126 . 128 . 130

Capítulo 1

Introducción a las Wavelets y prerrequisitos de álgebra lineal Trabajaremos con vectores en RN . Es decir, para cada señal discreta f f supondremos que f = (f 1 , f 2 , . . . , f N ) RN .

∈

En general N N (la dimensión del espacio, la cual coincide con el número de datos discretos que tenemos de cada señal) será una potencia de 2 (N ( N = 2k ) suficientemente grande. Recordamos también un par de conceptos sobre señales discretas que utilizaremos en el texto: La energía de una señal f es E (f ) f ) := f := f 12 + f 22 +

· · · + f N 2 .

El soporte de f se f se define como el conjunto sop(f sop(f )) := m : f m = 0 ,

{

 }

y el valor medio de f lo f lo tomamos como f 1 + f 2 + + f N N , sop(f sop(f ))

|

··· |

siendo sop(f sop(f )) el número de índices en el soporte de f . f .

|

|

1

2

Intro ducción a las Wavelets y prerrequisitos de álgebra lineal

En esta sección preliminar preliminar intentare intentaremos mos introducir, introducir, de un modo simple, simple, alguna de las nociones básicas que serán desarrolladas con rigor en los posteriores capítulos. Tomemos, por ejemplo, la siguiente señal f := := (1, (1, 7, 4, 3). 3). A partir de f generamos f generamos otra señal A1 de la siguiente manera: promediamos los dos primeros elementos de f f y los dos últimos, lo cual nos da 4 y 3.5 respectivamente. Si los dos primeros elementos son sustituidos por su promedio, al igual que los dos últimos, obtenemos la señal A1 := (4, (4, 4, 3.5, 3.5). 5). De esta forma la nueva señal viene a ser una aproximación de f de menor resolución . Es lo que llamaremos más adelante como señal de prome promedio dio de Haar a primer nivel . Examinemos ahora la diferencia entre f y f y la aproximación D1 := f := f

− A1 = (−3, 3, 0.5, −0.5). − 5).

Esta nueva señal D1 será la señal de detalle de Haar a primer nivel . Es precisamente lo que necesitamos añadir a la señal de promedio para recuperar la señal original f . f . Podemos observar que D 1 consiste en, hacer las diferencias promediadas entre los elementos 1 y 2 , 2 y 1 , 3 y 4 , y 4 y 3 respectivamente de f . f . Si quisiéramos continuar el proceso y pasar a una aproximación de incluso menor resolución podríamos promediar los cuatro elementos de f y f y obtener la señal constante A2 := (3. (3.75 75,, 3.75 75,, 3.75 75,, 3.75). 75). o

o

o

o

o

o

o

o

También deberíamos observar que este mismo resultado se s e hubiera obtenido 1 2 de proceder con A con A en vez de f de f .. La señal A señal A es el promedio de Haar a segundo nivel de f . f . La diferencia entre A entre A 1 y A 2 es la señal de detalle de Haar a segundo nivel, es decir, D2 := A := A 1 A2 = (0. (0.25 25,, 0.25 25,, 0.25 25,, 0.25). 25).

−

−

−

En definitiva podemos recuperar f a f a partir de A2 añadiendo añadiendo los detalles detalles f = A2 + D2 + D1 .

3

(a) f(t)

(b) Detalle de Haar a nivel 1

1

0.05

0.5 0

0 −0.5 −1 0

0.5

−0.05 0

1

(c) Promedio de Haar a nivel 2

0.5

1

(d) Detalle de Haar a nivel 2

1

0.05

0.5 0

0

−0.5 −1 0

0.5

−0.05 0

1

0.5

1

Figura 1.1: (a) Señal original f original f ((t). (b) Detalle de Haar a nivel 1. (c) Promedio de Haar a nivel 2. (d) Detalle de Haar a nivel 2.

En esto consiste el análisis de multirresolución . En la figura 1.1 representamos 2

10(t−0.5) f ( f (t) = 0.8e10(t sen(40t sen(40t)

muestreada a 1024 puntos en el intervalo [0 intervalo [0,, 1], 1], la señal D señal D 1 , diferencia entre f entre f y A1 (no representada) de promedio a primer nivel que consiste en promediar los elementos 1 elementos 1 y 2 , 3 y 4 , y así sucesivamente, de f de f y y sustituir los elementos de por los correspondientes promedios. La señal A señal A 2 formada por el promedio de f por los cuatro primeros elementos de f , f , situado en las cuatro primeras posiciones, luego el promedio de los elementos 5 a 8 , situado en las cuatro siguientes posiciones, y así sucesivamente. Por último la señal D2 de detalle a segundo nivel es la diferencia entre A1 y A2 . Como se puede comprobar, en el proceso anterior estamos introduciendo redundancia de información a cada paso del análisis de multirresolución que efectuamos, ya que partimos de una señal N -dimensional N -dimensional y vamos obteniendo sucesivas señales de la misma dimensión. Con la transformada wavelet que veremos con detalle en los próximos capítulos dicha redundancia se elimina y obtendremos en cada nivel una señal N -dimensional N -dimensional compuesta por varias o

o

o

o

o

o

4

Intro ducción a las Wavelets y prerrequisitos de álgebra lineal

subseñales que representan la información de las señales de promedio y detalles. A lo largo del texto damos por supuestos conocimientos básicos de algebra lineal. No obstante, en aras de la completitud, recordamos brevemente la notación y nociones que utilizaremos. Dado un sistema de vectores X = v1 , v2 , . . . , v m en RN , su envoltura lisubespacio V formado formado por las combinaciones lineales de los elementos neal es el subespacio V de X X y y lo denotamos por

{

}

V := lin := lin v1 , v2 , . . . , v m .

{

}

El producto escalar de dos vectores u = (u1 , . . . , u N ), w = (w1 , . . . , w N ) RN es u w := u := u 1 w1 + u2 w2 + + uN wN .

·

∈

···

Los vectores u vectores u y y w w son son ortogonales entre sí cuando su producto escalar es nulo. La norma euclídea de u es la raíz cuadrada de su energía, es decir



u := √ u · u = u21 + u22 + · · · + uN 2 . Un sistema de vectores X = {v1 , v2 , . . . , v m } en RN es ortonormal si cada

vector del sistema tiene norma 1 y son ortogonales entre sí. En ese caso son un sistema linealmente independiente y, al ser también sistema generador, forman una base de la envoltura lineal V lineal V de X de X ,, siendo cada elemento v elemento v V V expresable de forma única como combinación lineal de los elementos de X de de la siguiente manera v = (v v1 )v1 + (v ( v v2 )v2 + + (v (v vm )vm .

∈

·

·

···

·

Dos subespacios V y W son W son ortogonales entre sí cuando cada vector de V es ortogonal a todos los vectores de W . W . En ese caso al espacio suma de V y W generado por dichos subespacios lo denotamos como

⊕⊕⊥ W.

V

Si v Si v RN y W es W es un subespacio de RN , la proyección ortogonal de v de v sobre W sobre W es el (único) vector w W tal W tal que v w es ortogonal a W . W . También coincide con el vector w vector w W tal tal que la distancia de v de v a a w w,, dada por v w , es mínima. mínima. N N Una aplicación lineal T : R R es una transformación ortogonal si conserva el producto escalar, es decir, T decir, T v T w = v = v w para todo par v par v,, w RN . Esto equivale a que sea una isometría, o sea conserva la norma ( T v = v ,

∈

∈

∈

−

→

− 

·

·

∈    

5

para todo v RN ). En particular una transformación ortogonal preserva la energía, al ser el cuadrado de la norma. T norma. T es es transformación ortogonal si y sólo si su expresión matricial en una (cualquier) base ortonormal de RN es tal que sus filas (columnas) son una base ortonormal.

∈

Capítulo 2

Wavelets de Haar El tipo de wavelets (ondículas) que introducimos a continuación son las de Haar. Estas sirven de modelo por su sencillez y por contener la mayoría de los ingredientes que requerirán los restantes tipos de wavelets que iremos conociendo a lo largo del texto.

2.1 2.1

Scal Scalin ing g y wavelet eletss

Las wavelets de Haar de nivel 1 son:

 √ √ −   √ √ −   √ √ − 

1 , 0, . . . , 0 2 1 1 , 0, . . . , 0 w21 := 0, 0, , 2 2 .. . w11 :=

1 , 2

1 wN/2 N/2 := 0, . . . , 0,

1 , 2

1 2

.

El número de wavelets de Haar en el primer nivel es, por tanto, la mitad de la dimensión N . N . Observamos también que cada una de las wavelets w j1 , j > 1, se obtiene a partir de la wavelet w11 desplazando hacia la derecha 2 j posiciones el par de entradas no nulas 1 nulas 1// 2, 1/ 2. Por último hay que notar que conforman un sistema ortonormal y que el valor medio de cada wavelet es 0.

√ − √

7

8

Wavelets de Haar

Si tenemos una señal f RN , esta lleva asociada una subseñal primera fluctuación de Haar (coeficientes de detalle) expresada como:

∈

1 d

de

1 d1 = (d1 , . . . , d N/2 (f w11 , . . . , f wN/2 N/2 ) := (f N/2 ).

· ·

·

Es decir, es una subseñal con N/2 N/2 entradas que contiene todos los productos escalares de nuestra señal f f con cada una de las wavelets del primer nivel. Estas entradas serían: d1 =

− f 4 , . . . , d = f N N −√ − f 2 , d2 = f 3√ 1 − f N N √ . N/2 N/2 2 2 2

f 1

Por otra parte, los coeficientes de primera tendencia de f se f se definen como a1 =

f 3 + f 4 f N f 1 + f 2 N −1 + f N N , a2 = , . . . , a N/2 = . N/2 2 2 2

√

√

√

Estos también se pueden expresar como los productos escalares de f de f por por ciertas señales simples, llamadas scaling de Haar (funciones escala) de nivel 1, dadas por

√ √   √ √   √ √ 

1 1 , , 0, . . . , 0 2 2 1 1 v21 := 0, 0, , , 0, . . . , 0 2 2 .. . 1 1 1 , . vN/2 := 0 , . . . , 0 , N/2 2 2 v11 :=

Se tiene que 1 am = f = f vm , m = 1, . . . , N /2 /2,

· ·

conformando de esta manera la subseñal de primera tendencia a1 := (a1 , a2 , . . . , aN/2 N/2 ). En la figura 2.1 podemos ver la señal f ( f (t) = 10t 10t2 (2 t)6 sen(5πt sen(5πt)) en el intervalo [0, [0, 2], 2], muestreada con 1024 puntos (N ( N = = 1024) 1024) y las subseñales a1 y d1 asociadas. En las señales scaling de Haar se observan las siguientes similitudes con las respectivas wavelets: wavelets:

−

2.1 Scaling y wavelets

9

(a) f(t) 50

0

−50

0

1 00

200

3 00

400

50 0

60 0

(b) Subseñal tendencia a nivel 1

7 00

800

900

1000

(c) Subseñal detalle a nivel 1

50

2

1

0

0

−1

−50

0

1 00

20 0

300

4 00

−2

50 0

0

10 0

200

3 00

40 0

500

Figura 2.1: (a) Señal original f original f ((t). (b) Subseñal tendencia a nivel 1. (c) Subseñal fluctuación a nivel 1.

(a) Conforman Conforman un sistema ortonormal. ortonormal. 1 1 (b) Respecto los los soportes, soportes, sop( sop(vvm ) = sop(w sop(wm ), para cada m = 1, . . . , N /2 /2. 1 se obtiene a partir de la scaling v (c) Cada v Cada v m scaling v 11 desplazando hacia la derecha 2m posiciones el par de entradas no nulas 1/ 2, 1/ 2.

√ √

Pero se diferencian en que, mientras las wavelets tienen un valor medio igual a 0, en las scaling es 1/ 2. Además se cumple que

√

1 vn1 wm = 0,

·

/2. ∀ n, m = 1, . . . , N /2

Con lo cual la familia de vectores ortonormal de RN .

1 1 1 {v11, . . . , vN/2 N/2 , w1 , . . . , wN/2 N/2 } es una base

10

Wavelets de Haar

Si pasamos al nivel 2, las wavelets de Haar se definen como

     



1 1 1 1 , , , , 0, . . . , 0 w12 := 2 2 2 2 1 1 1 1 , 0, . . . , 0 w22 := 0, 0, 0, 0, , , , 2 2 2 2 .. . 1 1 1 1 2 , . wN/4 := 0 , . . . , 0 , , , N/4 2 2 2 2 Y las scaling son

v12 := v22 := .. . 2 vN/4 N/4

− −

− −

− −







1 1 1 1 , , , , 0, . . . , 0 2 2 2 2 1 1 1 1 0, 0, 0, 0, , , , , 0, . . . , 0 2 2 2 2

1 1 1 1 := 0, . . . , 0, , , , 2 2 2 2





.

Ahora definiremos las correspondientes subseñales de fluctuación y tendencia a nivel 2 como 2 d2 := (f w12 , . . . , f wN/4 N/4 )

· · · 2 a2 := (f · · v12, . . . , f · vN/4 N/4 ).

En la figura 2.2 podemos ver la señal f ( f (t) = 10t 10t2 (2 t)6 sen(5πt sen(5πt)) en el intervalo [0, [0, 2], 2], muestreada con 1024 puntos (N (N = 1024) 1024) y las subseñales a1 y d1 asociadas en la segunda fila y las subseñales a2 , d2 y d1 asociadas en la tercera tercera fila. Se puede observar, en primer lugar, que el soporte de las scaling y las wavelets a nivel 2 consta de 4 entradas consecutivas, conforman un sistema ortonormal, las wavelets a nivel 2 siguen teniendo un valor medio igual a 0 2 y en las scaling es 1/2. Las wm se obtienen desplazando hacia la derecha 4m 2 posiciones las entradas no nulas de w12 , al igual que las vm respecto v12 . Pero también hay que notar que, tanto las wavelets como las scaling a nivel 2, se

−

2.1 Scaling y wavelets

11

(a) f(t) 50

0

−50

0

1 00 200 3 00 400 (b) Subseñal tendencia a nivel 1

50 0

50

60 0

7 00 800 900 (c) Subseñal detalle a nivel 1

1 00 0

2 1

0

0 −1

−50

0

10 0 2 00 (d) Tend. nivel 2

50

0

−50 0

1 00

20 0

3 00

400 500 (e) Det. nivel 2

−2

2

2

1

1

0

0

−1

−1

−2

0

10 0

−2

20 0

0

10 0 20 0 3 00 4 00 (f) Subseñal detalle a nivel 1

500

0

10 0

500

20 0

3 00

4 00

Figura 2.2: (a) Señal original f original f ((t). (b) Subseñal tendencia a nivel 1. (c) Subseñal fluctuación a nivel 1. (c) Subseñal tendencia a nivel 2. (d) Subseñal fluctuación a nivel 2. (e) Subseñal fluctuación a nivel 1.

obtienen a partir de las scaling de nivel 1 por las formulas v12

v11 + v21 = 2 .. .

2 vN/4 N/4 =

w12

√

=

v11

− v21 √ 2

.. .

1 1 vN/2 N/2−1 + vN/2 N/2

2 wN/4 N/4 =

√ 2

(2.1) 1 vN/2 N/2−1

1 − vN/2 N/2 √ 2 .

Otra forma de ver la situación es que wavelets y scaling a nivel 2 se deducen de las wavelets y scaling a nivel 1 “dilatando” el soporte (que ahora es el doble) y “contrayendo” sus valores al multiplicarlos por el factor 1/ 1 / 2 (para preservar la energía unitaria). 2 2 2 De las igualdades (2.1) deducimos que v12 , . . . , v N/4 N/4 , w1 , . . . , w N/4 N/4 es una base ortonormal del subespacio

√

{

}

1 V 1 := lin v11 , . . . , v N/2 N/2 .

{

}

Todas estas observaciones nos acercan al análisis de multirresolución, objetivo de nuestra siguiente sección.

12

Wavelets de Haar

2.2 2.2

Análi Análisis sis de mul multir tirre reso solu luci ción ón (MR (MRA) A)

Constituye el eje central en el tratamiento de señales mediante wavelets. A partir de la señal original la descomponemos en suma de otras dos, una de promedio y otra de detalles. En el siguiente paso aplicamos el proceso anterior a la señal de promedio, y así sucesivamente. Para el proceso inverso de síntesis, empezamos por una señal de baja resolución, añadiendo los detalles sucesivamente con lo cual vamos alcanzando señales de mayor resolución, para acabar con la síntesis completa de la señal a la resolución más alta. Esto es lo que se conoce como análisis de multirresolución (MRA). 1 1 1 Hemos señalado que v11 , . . . , v N/2 N/2 , w1 , . . . , w N/2 N/2 es una base ortonormal de RN . Por tanto una señal f la f la expresaremos como

{

}

f = (f v11 )v11 +

1 1 1 1 · · · + (f · vN/2 · w11)w11 + · · · + (f · wN/2 (f · ( f · (f · N/2 )vN/2 N/2 + (f N/2 )wN/2 N/2 1 1 1 1 = a 1 v11 + · · · + aN/2 N/2 vN/2 N/2 + d1 w1 + · · · + dN/2 N/2 wN/2 N/2 dN/2 aN/2 d1 −d1 a1 a1 N/2 N/2 −dN/2 N/2 aN/2 N/2 . + √ , √ , . . . , √ , √ = √ , √ , . . . , √ , √ 2 2 2 2 2 2 2 2

· ·







Esto nos lleva al primer nivel de MRA de Haar f = A 1 + D1 , siendo A1 la señal del primer promedio A1 =

 √ √

  √ · ·

aN/2 a1 a1 N/2 aN/2 N/2 , ,..., , 2 2 2 2

√

N/2 N/2

(f v j1 )v j1 ,

=

j=1 j =1

y D 1 la señal de primer detalle D1 =

 √ −√ d1 , 2

dN/2 d1 N/2 ,..., , 2 2

dN/2 N/2 √ −√ 2

  · · N/2 N/2

(f w j1 )w j1 .

=

j=1 j =1

Por tanto la señal de promedio es una combinación lineal de las scaling, tomando los valores de la subseñal de primera tendencia como coeficientes, mientras que la señal de detalle es combinación de las wavelets, tomando las fluctuaciones de nivel 1 como coeficientes. En realidad A1 es la proyección de f sobre f sobre el espacio V 1 , y D 1 , por ser la difere diferenci nciaa entre entre f y A1 , es la proyección de f f sobre sobre el espaci espacioo W 1 := 1 lin w11 , . . . , w N/2 N/2 .

{

}

2.2 Análisis de multirresolución (MRA)

13

Si expresamos A1 y D1 en función de los valores de f tenemos f tenemos A1 = D1 =

 

 

f N f 1 + f 2 f 1 + f 2 f 3 + f 4 f 3 + f 4 N −1 + f N N f N N −1 + f N N ,..., , , , , , 2 2 2 2 2 2 f N f N f N f 1 f 2 f 2 f 1 f 3 f 4 f 4 f 3 N −1 N f N N N −1 , . , , , ,..., 2 2 2 2 2 2

−

−

−

−

−

−

Por tanto f es f es la suma de una señal de resolución más baja A1 y de la señal de detalles D 1 . Si vamos al segundo nivel del MRA de f f = A 2 + D2 + D1 , donde A2 = (f v12 )v12 +

2 2 · · · · · + (f · vN/4 (f · N/4 )vN/4 N/4 , 2 2 D2 = (f · (f · · w12)w12 + · · · + (f · wN/4 N/4 )wN/4 N/4 , 2 2 2 1 1 ya que { v12 , . . . , v N/4 N/4 , w1 , . . . , w N/4 N/4 , w1 , . . . , w N/2 N/2 } es una base ortonormal de

RN .

En este caso A2 es la señal del segundo promedio de f , f , que consiste en proyectar f (o, f (o, equivalentemente, A1 ) sobre el espacio 2 V 2 := lin v12 , . . . , vN/4 N/4

{

}

y D2 es la proyección de f (o A1 ) sobre 2 W 2 = lin w12 , . . . , w N/4 N/4 .

{

}

Expresados en función de los valores de f f tenemos que A2 =

D2 =

 

f 1 + f 2 + f 3 + f 4 f 1 + f 2 + f 3 + f 4 f 1 + f 2 + f 3 + f 4 , , , 4 4 4 f 1 + f 2 + f 3 + f 4 f N N −3 + f N N −2 + f N N −1 + f N N ,..., , 4 4



f 1 + f 2

− f 3 − f 4 , f 1 + f 2 − f 3 − f 4 , −f 1 − f 2 + f 3 + f 4 , 4

−f 1 −

4 f 2 + f 3 + f 4 ,..., 4

4

−f N N −3 − f N N −2 + f N N −1 + f N N 4



.

14

Wavelets de Haar

1

(a) f(t)

(b) D

1

0.05

0.5 0

0 −0.5 −1 0

0.5

1

−0.05 0

0.5

2

1 2

(c) A

(d) D

1

0.05

0.5 0

0

−0.5 −1 0

0.5

1

−0.05 0

0.5

1

Figura 2.3: (a) Original f Original f ((t). (b) Detalle D Detalle D 1 . (c) Promedio A Promedio A 2 . (d) Detalle D Detalle D 2 .

Así reconstruimos la señal f a f a partir de una señal de baja resolución A2 añadiendo sucesivamente los detalles D2 y D1 . En la figura 2.3 representamos la señal f ( f (t) = 20t 20t2 (1 t)4 cos(12πt cos(12πt)) muestreada a 1024 puntos en el intervalo [0, [0, 1] en 1] en (a), la señal D 1 de primer detalle en (b), la señal D señal D 2 de segundo detalle detalle en (d) y la señal A señal A2 de segundo promedio en (c).

−

Continuaríamos el proceso de forma inductiva y el MRA a nivel m de f de f es f = A m + Dm +

· · · + D 2 + D1 ,

donde m donde m k y N y N = 2k (por ejemplo, si la señal original es de N = N = 1024 datos, 1024 datos, entonces podemos llegar hasta un nivel máximo de k = 10 en el MRA). El promedio y el detalle a este nivel es

≤

Am = (f v1m )v1m +

m m · · · · · + (f · vN/2 (f · N/2 )vN/2 N/2 , m m · w1m)w1m + · · · + (f · wN/2 Dm = (f · (f · N/2 )wN/2 N/2 m

m

m

m

,

2.2 Análisis de multirresolución (MRA)

15

siendo las wavelets y scaling a nivel m definidas por 1 1 1 1 , . . . , , 0, . . . , 0) , . . . , , w1m := ( m/2 m/2 2m/2 m/2 m/2 2m/2 2m/2 2m/2

−

−

      − −            2m−1

2m−1

w2m := (0, (0, . . . , 0, 2m

.. .

1 1 1 , 0, . . . , 0) , , . . . , , . . . , m/2 m/2 m/2 m/2 2m/2 2m/2 2m/2 2m/2 1

2m−1

2m−1

m wN/2 (0, . . . , 0, N/2m := (0,

1

2

,..., m/2 m/2

1 1 ) , , . . . , m/2 m/2 m/2 2m/2 2m/2 2m/2 1

−

−

               2m−1

2m−1

v1m := (

1

1 , . . . , , 0, . . . , 0) m/2 m/2 2m/2 2m/2 2m

v2m := (0, (0, . . . , 0, 2m

.. .

1 , 0, . . . , 0) , . . . , m/2 m/2 2m/2 2m/2 1

2m

m vN/2 (0, . . . , 0, N/2m := (0,

1

1 , . . . , ). m/2 m/2 2m/2 2m/2

   2m

Como ocurría en el caso de las wavelets y scaling de nivel 2, que se expresaban en función de las scaling de nivel 1, las wavelets y scaling de nivel m son m son combinaciones lineales de las scaling de nivel m 1 bajo las mismas fórmulas de recurrencia, concretamente

−

v1m

v1m−1 + v2m−1 = 2 .. .

m vN/2 N/2m =

w1m

√

m−1 vN/2 N/2m−

1

=

v1m−1

− v2m−1 √ 2

.. .

m−1 + vN/2 −√ 1 N/2 − m

1

m wN/2 N/2m =

(2.2) m−1 vN/2 N/2m−

1

m−1 − vN/2 −√ 1 N/2 − . m

1

2 2 En definitiva, la señal de m-ésimo promedio A promedio A m , es la proyección ortogonal de f sobre f sobre el subespacio m V m := lin v1m , . . . , vN/2 N/2m ,

{

}

16

Wavelets de Haar

mientras que la señal D m de detalles a nivel m es la proyección de f sobre f sobre m W m := lin w1m , . . . , wN/2 N/2m .

{

}

Por las relaciones anteriores, los subespacios V m y W m están contenidos en V m−1 . Es más, cada uno de ellos es el complemento ortogonal del otro en V m−1 , es decir, V m ⊥ W m = V m−1 .

⊕

Esto nos da la siguiente descomposición, que es precisamente el MRA a nivel m. ⊥ W 1 . RN = V m ⊥ W m ⊥

⊕

⊕ ···⊕ Como hemos observado anteriormente, m ≤ k (N = 2k ). De manera que al

llegar al último nivel k sólo tendríamos una wavelet y una scaling w1k = (

−1 , . . . , √ −1 ), √ 1N , . . . , √ 1N , √ N N

      √ √    N/2 N/2

v1k = (

N/2 N/2

1 1 ). ,..., N N N

Con lo cual nos daría, dada una señal f , f , la resolución más baja posible Ak consiste en la constante Ak = (

f 1 + f 2 + N

· · · + f N N , . . . , f 1 + f 2 + · · · + f N N ), N

que, en definitiva, es tomar el valor promedio de toda la señal f y f y considerar la señal constante asociada. Si añadimos de forma sucesiva los detalles D detalles D k , D k−1 , . . . , D1 , a AK recuperamos la original f . f . En la descomposición del MRA para cada f lo f lo único que hay que calcular son los coeficientes que multiplican a las wavelets y scaling, que son las tendencias y fluctuaciones. Esto nos conduce de forma natural a la transformada de Haar.

2.3 2.3

Trans ransfo form rmad ada a de Haar Haar

La transformada de Haar a nivel m es sencillamente una transformación tal que a cada señal f le f le hace corresponder la señal consistente en los coeficientes

2.3 Transformada de Haar

17

de las tendencias de nivel m, junto con las fluctuaciones de nivel m e inferior ordenadas, es decir, f

−H→ (am|dm|dm−1| . . . |d1). m

Esta transformación es, en primer lugar, lineal. Además, teniendo en cuenta como se obtienen los coeficientes, es una transformación ortogonal. Escrito matricialmente, para el nivel 1 obtendríamos

  

a1

d1

  

=

  

v11 .. .

1 vN/2 N/2 w11 .. . 1 wN/2 N/2

     

f 1 .. . .. . .. . .. . f N N

   

Se observa que la matriz asociada a dicha transformación es tal que sus filas forman una base ortonormal. Es por ello que la transformación H1 es ortogonal. Esto es una gran ventaja a la hora de realizar cálculos ya que, por ejemplo, la transformada inversa 1

1

(a d )

1 H− 1

| −→ f

se calcula de forma trivial simplemente trasponiendo la matriz:

  

v11 .. . 1 vN/2 N/2 w11 .. . 1 wN/2 N/2

      T

a1

d1

  

=

   

f 1 .. . .. . .. . .. . f N N

   

18

Wavelets de Haar

De forma análoga, la transformación H transformación H m a nivel m nivel m también también es ortogonal, sería

    

am

dm .. . d1

    

=

     

v1m .. . m vN/2 N/2m w1m .. . m wN/2 N/2m .. . w11 .. . 1 wN/2 N/2

           

f 1 .. . .. . .. . .. . .. . .. . .. . .. . f N N

     

y análogamente podemos invertir el proceso de forma sencilla. A efectos gráficos, como hacer la transformada de nivel 1 produce H1 (f ) f ) = (

f N f N f 1 + f 2 f N f 2 N −1 N N −1 + f N N f 1 ), ,..., , ,..., 2 2 2 2



√



N/2 N/2

√

 

− √



N/2 N/2

√ −



entonces esta transformada viene a ser una “compresión” de la señal original en la primera mitad, multiplicada por el factor 2, mientras que la segunda mitad contiene contiene las fluctuacion fluctuaciones, es, las cuales serán “pequeñas” “pequeñas” si la señal original no tiene cambios “bruscos”. A niveles superiores, en cada paso sucesivo esencialmente comprimimos a la mitad la subseñal de tendencia que teníamos en el paso anterior y lo situamos en la primera parte, a continuación las fluctuaciones de dicha subseñal, y por último las fluctuaciones que ya había en el paso anterior. Hay dos propiedades fundamentales de esta transformada: Por una parte la conservación de la energía al ser una transformación ortogonal. Por otra parte se cumple la característica de las fluctuaciones pequeñas, es decir, si una señal es constante en el soporte de una wavelet, al multiplicarla escalarmente por ella nos saldrá 0, es decir, la correspondiente fluctuación es nula. Si la señal no “varía” mucho a lo largo del soporte de la wavelet, entonces la fluctuación asociada es “pequeña”.

√


19

La conjunción de las dos observaciones anteriores conlleva el fenómeno conocido como compactación de la energía . Concretam Concretamen ente, te, la energía energía de la señal original es la suma de las energías de las subseñales de tendencia y fluctuaciones E (f ) = E (am ) + E (dm ) +

· · · + E (d1).

Ahora bien, por la característica de las fluctuaciones pequeñas, será la tendencia quien concentrará en general la mayor parte de la energía de la transformada. Una manera de cuantificarlo es mediante el perfil de energía acumulada :



f 12 E (f ) f )

,

f 12

+ f 22

E (f )

,...,

f 12

+

···

2 + f N 1

E (f ) f )



− ,1 .

En la figura 2.4 tenemos la señal original (a), la transformada de nivel 2 (b), y los perfiles de energía acumulada de ambas ((c) y (d)). (a)

(b)

1.5

1.5

1

1

0.5

0.5

0

0

−0.5

−0.5

−1

−1

−1.5 0

1000

2000

3000

4000

−1.5 0

1000

(c) 1.5

1

1

0.5

0.5

0

0

1000

2000

3000

4000

3000

4000

(d)

1.5

−0.5 0

2000

3000

4000

−0.5 0

1000

2000

Figura 2.4: (a) Señal original. (b) Transformada de Haar de nivel 2. (c) Perfil de energía acumulada de la señal. (d) Perfil de energía acumulada de la transformada de Haar de nivel 2. Como se puede apreciar claramente, la energía se acumula antes en la transformada que en la original. De hecho en el primer cuarto del intervalo, que corresponde precisamente a la subseñal a2 , se acumula el mayor porcentaje de la energía.

20

Wavelets de Haar

Por último podemos calcular el número de operaciones necesario para realizar la transformada de Haar. La idea fundamental reside en el proceso de obtención de Hm (f ) f ), que paso a paso sigue el siguiente diagrama, conocido como algoritmo en pirámide de Mallat. f

−H→ ց

a1

1

d1

−H→ ց

a2 . . .

1

−H→ ց

am−1

1

d2

am dm

Si esto lo expresamos matricialmente tenemos

   

m

a dm .. . .. . d1

   

=

   

H 1 ( 2mN − ) 1

0

.

donde

H1 ( j) j ) =

0

   

   

I

   

√ 12

√ 12

0 ... 0

0

√ 12

√ 12

0

√ 1

− √ 1

0 0

0 0

√ 1

− √ 12

2

0 ... 0

0

0

2

2

0 0

0

√ 12

− √ 12

0

0

0

0 ... 0

√ 1

0 0

0 0

1 2

0 0

0 0

− √

√ 12

... ...

0

√

... ... ...

0

√ 12

1 2

0 0

0

√ 12

1 2

0 0

...

0

I

... ...

0

0 ... 0

0

0

0

√ 1

2

0

H 1 ( N 2 )

0

√ 1

2

...

   

√

0 0

... ... ...

√ 1

√ 1

0 0

0 0

...

√ 1

− √ 12

2

2

2

2

1 2

− √

   

               

.

f 1 f 2 .. . .. . f N N −1 f N N

j/ j /2

j/ j /2

   


21

Multiplicar por la primera matriz requiere 3N operaciones, N operaciones, por la segunda 3N/2 N/2, por la tercera 3N/4 N/4, . . . , por la m-ésima 3N/2 N/2m−1 . Esto nos da un número de operaciones no superior a 6 a 6N N ,, es decir O decir O((N ) N ) operaciones. Resultado que debe ser comparado con otras transformadas rápidas, como la FFT que requiere O(N log log N ) N ) operaciones.

Capítulo 3

Familias de wavelets ortogonales 3.1 3.1

Wavelets elets de Daubec Daubechi hies es

En esta sección describimos una amplia colección de tipos de wavelets ortogonales descubiertas por Daubechies. La diferencia respecto a las de Haar es la forma en como se definen las scaling y las wavelets. Las propiedades más importantes del MRA y la transformada seguirán siendo válidas en el nuevo contexto. Las ventajas que se obtendrán con las nuevas wavelets residen en el hecho que hay una mejora substancial en el tratamiento de señales que posean cualidades adecuadas de “suavidad”.

3.1. 3.1.1 1

Daube Daubecchies hies

db2

Conceptos básicos En primer lugar introducimos las scaling y wavelets Daubechies2 (db2), a nivel 1. v11 := (α1 , α2 , α3 , α4 , 0, 0, . . . , 0)

w11 := (β 1 , β 2 , β 3 , β 4 , 0, 0, . . . , 0)

v21 := (0, (0, 0, α1 , α2 , α3 , α4 , 0, . . . , 0) .. .

w21 := (0, (0, 0, β 1 , β 2 , β 3 , β 4 , 0, . . . , 0) .. .

1 vN/2 (0, . . . , 0, α1 , α2 , α3 , α4 ) N/2−1 := (0,

1 wN/2 (0, . . . , 0, β 1 , β 2 , β 3 , β 4 ) N/2−1 := (0,

1 vN/2 N/2 := (α3 , α4 , 0, . . . , 0, α1 , α2 )

1 wN/2 N/2 := (β 3 , β 4 , 0, . . . , 0, β 1 , β 2 ).

23

24

Familias de wavelets ortogonales

Donde las entradas en las scaling son

√ 1+ 3 α1 := √ , 4 2

√ 3+ 3 α2 := √ , 4 2

√ 3− 3 √ , α3 := 4 2

√ 1− 3 √ . α4 := 4 2

β 3 = α = α 2 ,

−α1.

y las entradas en las wavelets son β 1 = α4 ,

β 2 =

−α3,

β 4 =

Estos valores son tales que α21 + α22 + α23 + α42 = 1, α1 α3 + α2 α4 = 0. 1 1 1 Estas relaciones implican que v11 , . . . , v N/2 N/2 , w1 , . . . , w N/2 N/2 sea una base ortonormal de RN . Por otra parte también se cumple

{

}

β 1 + β 2 + β 3 + β 4 = 0, lo cual conlleva que el valor medio de las wavelets sea 0. Posteriormente resaltaremos una igualdad añadida respecto a las de Haar que cumplen los elementos de las wavelets y que les dotan de una mayor ventaja en el tratamiento de señales suaves. Dada un señal f señal f ,, definiendo las subseñales de primera tendencia y fluctuación, respectivamente, al igual que lo hicimos con las de Haar, es decir, 1 a1 := (f v11 , f v21 , . . . , f vN/2 N/2 ),

· · · · · 1 d1 := (f · · w11, f · · w21, . . . , f · wN/2 N/2 ),

se observa que los elementos de la subseñal de tendencia vuelven a ser una media ponderada, esta vez entre 4 valores consecutivos de f , f , mientras que las fluctuaciones son diferencias ponderadas. Para niveles superiores se definen las correspondientes scaling y wavelets a partir de las del nivel anterior, utilizando los coeficientes nuevos dados por las

3.1 Wavelets de Daubechies

25

fórmulas recurrentes v1m = α 1 v1m−1 + α2 v2m−1 + α3 v3m−1 + α4 v4m−1 .. . m m−1 vN/2 N/2m = α 1 vN/2 N/2m−

m−1 m−1 m−1 + α v + α v 3 4 − 1 2 −1 + α2 vN/2 N/2 wm = β 1 v m−1 + β 2 v m−1 + β 3 v m−1 + β 4 v m−1 1

m

1

1

2

1

3

4

.. . m−1 m wN/2 N/2m = β 1 vN/2 N/2m−

11 , 2

Wavelets w

1

11 y 4

w

m−1 + β 3 v1m−1 + β 4 v2m−1 . −1 + β 2 vN/2 N/2 − m

11 7

w

1

Scaling v

11 , 2

11 y 4

v

11 7

v

Figura 3.1: Para N = 214 y db2 a la izquierda las wavelet w211 w411 y w711 y a la derecha las scaling v211 , v411 y v711 . La figura 3.1 nos muestra, para N = 214 y m = 11 11,, la gráfica de varios vectores wavelet y scaling de la db2. Como en el caso de Haar, el MRA a nivel m nivel m de de una señal f señal f se se define como f = A m + Dm + donde m

· · · + D 2 + D1 ,

≤ k − 1 (N = 2k ), y las señales promedio y detalle son m m Am = (f · (f · · v1m)v1m + · · · + (f · vN/2 N/2 )vN/2 N/2 m m Dm = (f · (f · · w1m)w1m + · · · + (f · wN/2 N/2 )wN/2 N/2 . m

m

m

m

Por tanto, la señal de m-ésimo promedio Am , es también la proyección ortogonal de f sobre f sobre el subespacio m V m := lin v1m , . . . , vN/2 N/2m ,

{

}

26


mientras que la señal de detalles a nivel m, Dm , es la proyección de f sobre f sobre m W m := lin w1m , . . . , wN/2 N/2m .

{

}

En definitiva llegamos a la descomposición RN = V m

⊕⊥ W m ⊕⊥ · · · ⊕⊥ W 1. Como hemos observado anteriormente, m ≤ k − 1 (N = 2k ). De manera que al llegar al último nivel k − 1 solo tendríamos un par de wavelets y un par de

scaling. Al igual que en la transformada de Haar, la transformada db2 a nivel m es una transformación tal que a cada señal f f le hace corresponder la señal consistente en la tendencia de nivel m, junto con las fluctuaciones de nivel m e inferior ordenadas, es decir, f

db2m

−→ (am|dm|dm−1| . . . |d1).

Esta transformación que también es ortogonal y escrita matricialmente, para el nivel m, se tendría

    

am

dm .. . d1

    

=

     

v1m .. . m vN/2 N/2m w1m .. . m wN/2 N/2m .. . w11 .. . 1 wN/2 N/2

           

f 1 .. . .. . .. . .. . .. . .. . .. . .. . f N N

     

Una vez más, por ser una transformación ortogonal, se conserva la energía (E (db2 (db2m (f )) f )) = E (f ) f )) y se da la propiedad de las “fluctuaciones pequeñas”, lo que implica el fenómeno de compactación de la energía. Además la transformada se invierte de forma sencilla trasponiendo la matriz correspondiente en el nivel en el cual dispongamos de la señal transformada.


27

También aquí podemos obtener cada una de las subseñales de tendencia y fluctuación de un determinado nivel a partir de la subseñal de tendencia del nivel anterior aplicando la transformada db21 , resultando igualmente como en el caso de Haar un diagrama en forma de pirámide. f

db21

−→ ց

a1 d1

db21

−→ ց

a2 . . .

am−1

d2

db21

−→ ց

am dm

Característica de las fluctuaciones pequeñas Si una señal es lineal en el soporte de una wavelet db2, el producto escalar de la señal por dicha wavelet es 0, es decir, la correspondiente fluctuación es nula. Si la señal es “aproximadamente” lineal a lo largo del soporte de la wavelet, entonces la fluctuación asociada es “pequeña”. La razón de que, en este caso, obtengamos una mejora respecto a las wavelets de Haar es consecuencia de la siguiente igualdad 1β 1 + 2β 2 β 2 + 3β 3 β 3 + 4β 4 β 4 = 0. Indicamos la prueba para el nivel 1. El caso general es análogo. Sea i, i + 1, i + 2, i + 3 el soporte de cierta wavelet w wavelet w j1 . Si f Si f es es lineal a lo largo de dicho soporte, entonces existen a, b R tales que

}

{

∈

f l = a + b(l

1), l = i = i,, . . . , i + 3. 3. − i + 1),

Entonces f w j1 = β 1 (a + b) + β 2 (a + 2b 2b) + β 3 (a + 3b 3b) + β 4 (a + 4b 4b) = 0.

· ·

Esta observación resulta particularmente importante cuando estamos tratando una señal discreta obtenida a partir del muestreo de una señal analógica que es de clase C 2 , es decir, dos veces derivable en el intervalo y cuya derivada segunda es continua. En este caso, por el desarrollo de Taylor, sabemos que siendo f n = f ( f (nh) nh) (suponemos, por comodidad, que el muestreo es en un intervalo cuyo extremo inferior es el 0) con h el paso de muestreo, se cumple f i+n = f n + f ′ (nh) nh)ih + O(h2 ), i = 0, 1, 2, 3.

28


Y así llegamos a f w jm = O( O (h2 ). En las mismas condiciones para la señal f , f , se obtiene que f w jm = O( O (h) para las wavelets de Haar, lo cual es claramente peor.

·

· ·

Si juntamos esta última ecuación con las que ya teníamos anteriormente, obtenemos el sistema: β 12 + β 22 + β 32 + β 42 = 1, β 1 + β 2 + β 3 + β 4 = 0, β 1 β 3 + β 2 β 4 = 0, 1β 1 + 2β 2 β 2 + 3β 3 β 3 + 4β 4 β 4 = 0. A partir de estas ecuaciones se determinan los valores de los coeficientes de db2. También aquí podemos calcular el número de operaciones necesario para realizar la transformada db2. La expresión matricial de la transformada es

   

am dm .. . .. . d1

   

=

   

db21 ( 2mN − )

0

1

0

.

   

α1 0 ... α3 β 1 0 ... β 3

I

α2 α3 α4 0 α1 α2

   

...

0 α3

   

db21 ( N 2 )

0

0

0 ... α3 . . .

I

0 0

0 0

α4 0 β 2 β 3 0 β 1

0 ... β 4 0 β 2 β 3

0 0 0 ... β 4 . . .

α1 α2 0 0 0 0

β 4

0

0

β 1

0

...

0

β 2

      

   

.

f 1 f 2 .. . .. .

f N N −1 f N N

   


29

donde

db21 ( j) j ) =

   

α1 0 ... α3 β 1 0 ... β 3

α2 α3 α4 0 α1 α2

0 α3

0 ... α3 . . .

0 0

α4 0 β 2 β 3 0 β 1

0 ... β 4 0 β 2 β 3

0 0 0 ... β 4 . . .

α1 α2 0 0 0 0

β 4

0

0

β 1

0

...

0

0 0

β 2

   

   

j/ j /2

j/ j /2

Y nos da un computo total de O(N ) N ) operaciones.

3.1. 3.1.2 2


db3

El buen comportamiento que presenta la transformada db2 para señales “suaves” se puede mejorar añadiendo más condiciones, y por tanto más números en las scaling o wavelets.

Conceptos básicos Los valores que definen las scaling de db3 son α1 := 0.332670552950083 α2 := 0.806891509311092 α3 := 0.459877502118491

α4 :=

−0.135011020010255 α5 := −0.0854412738820267 α6 := 0.0352262918857095

y los correspondientes para las wavelets se toman como β 1 := α := α 6 , β 2 :=

−α5,

β 3 := α := α 4 , β 4 :=

−α3,

β 5 := α := α 2 , β 6 :=

−α1.

30


Con lo cual definimos las scaling y wavelets a nivel 1. v11 := (α1 , α2 , α3 , α4 , α5 , α6 , 0, 0, . . . , 0) (0, 0, α1 , α2 , α3 , α4 , α5 , α6 , 0, . . . , 0) v21 := (0, .. . 1 vN/2 N/2−1 := (α5 , α6 , 0, . . . , 0, α1 , α2 , α3 , α4 ) 1 vN/2 N/2 := (α3 , α4 , α5 , α6 , 0, . . . , 0, α1 , α2 )

w11 := (β 1 , β 2 , β 3 , β 4 , β 5 , β 6 , 0, 0, . . . , 0) w21 := (0, (0, 0, β 1 , β 2 , β 3 , β 4 , β 5 , β 6 , 0, . . . , 0) .. . 1 wN/2 N/2−1 := (β 5 , β 6 , 0, . . . , 0, β 1 , β 2 , β 3 , β 4 )

1 wN/2 N/2 := (β 3 , β 4 , β 5 , β 6 , 0, . . . , 0, β 1 , β 2 ).

Para, finalmente, generar las de niveles superiores mediante las correspondientes fórmulas v1m = α 1 v1m−1 + α2 v2m−1 + α3 v3m−1 + α4 v4m−1 + α5 v5m−1 + α6 v6m−1 .. . m−1 m vN/2 N/2m = α 1 vN/2 N/2m−

m−1 + α3 v1m−1 + α4 v2m−1 + α5 v3m−1 + α6 v4m−1 N/2 − −1 + α2 vN/2 wm = β 1 v m−1 + β 2 v m−1 + β 3 v m−1 + β 4 v m−1 + β 5 v m−1 + β 6 v m−1 1

1

1

1

m

2

3

4

5

6

.. . m−1 m wN/2 N/2m = β 1 vN/2 N/2m−

1

m−1 + β 3 v1m−1 + β 4 v2m−1 + β 5 v3m−1 + β 6 v4m−1 . −1 + β 2 vN/2 N/2 − m

1

A partir de aquí se definen las subseñales de tendencia, de fluctuaciones, la descomposición ortogonal del MRA y la transformada como en el caso de Haar y db2. Pero en esta ocasión se cumple una ecuación lineal más, concretamente se tienen que β 1 + β 2 + β 3 + β 4 + β 5 + β 6 = 0, 1β 1 + 2β 2 β 2 + 3β 3 β 3 + 4β 4 β 4 + 5β 5 β 5 + 6β 6 β 6 = 0, 12 β 1 + 2 2 β 2 + 3 2 β 3 + 4 2 β 4 + 5 2 β 5 + 6 2 β 6 = 0. Esta última ecuación es precisamente la que permite mejorar el tratamiento de señales con buenas condiciones de diferenciabilidad.


31

Característica de las fluctuaciones pequeñas Si una señal es cuadrática en el soporte de una wavelet db3, al multiplicarla escalarmente por ella nos saldrá 0, es decir, la correspondiente fluctuación es nula. Si la señal es “aproximadamente” cuadrática a lo largo del soporte de la wavelet, entonces la fluctuación asociada es “pequeña”. En particular, si partimos de una señal que se obtiene como muestreo de una analógica de clase C clase C 3 en un cierto intervalo, intervalo, y tomamos los puntos equidistantes a distancia h, entonces las fluctuaciones son O(h3 ). La implicación clara de este fenómeno es que, en lo que respecta a la compresión y eliminación del ruido de señales de estas características, la transformada db3 tendrá un mejor comportamien comportamiento to que la db2 o la de Haar. Sin embargo si, por ejemplo, lo que se pretendiera es localizar los puntos donde la pendiente de la señal cambia más bruscamente (donde hay “picos”), entonces db2 consigue ese objetivo de forma más eficiente que db3.

3.1. 3.1.3 3


dbJ

El resto de transformadas wavelet de Daubechies, dbJ con J = 4, 5, . . . , 10 10,, se definen de la misma manera cumpliendo la anulación de momentos de orden superior. Por lo que los números scaling α1 , . . . , α 2J son tales que α21 + α22 +

· · · + α22J = 1√ , α1 + α2 + · · · + α2J = 2 Y los números wavelet β 1 , . . . , β 2 J se toman como β 1 = α = α 2J , β 2 =

−α2J −1, . . . , β2 J −1 = α2, β 2J = −α1.

Que además han de cumplir las ecuaciones de anulación de momentos β 1 + 2 n β 2 + 3 n β 3 +

· · · + (2J (2J )n β 2J = 0, n = 0, 1, . . . , J − 1,

junto con las de ortogonalidad. La anulación de momentos es lo que permite obtener fluctuaciones pequeñas. Concretamente, si tratamos mediante dbJ una señal obtenida como muestreo de una analógica de clase C J , tomando el paso igual a h, entones las fluctuaciones son O(hJ ).

32

3.2


Coiflets

Finalmente damos una descripción de las coiflets. Siguiendo la sugerencia de Coifman, este tipo de wavelets ortogonales fueron construidas por Daubechies. El objetivo era conseguir una aproximación mejorada entre los valores de las tendencias y los de la señal original. Si definimos, en primer lugar, los números scaling scaling para coif1

−√ √ 7 , 2 √ −√ 2 7 ,

1 α1 := 16 14 α4 := 16

2

√ √ √ −√

5+ 7 , α2 := 16 2 1 7 , α5 := 16 2

√ √ √ − √

14 + 2 7 , α3 := 16 2 3+ 7 α6 := . 16 2

Entonces las scaling y wavelets a nivel 1 se construyen, a diferencia de las db3 , como v11 := (α3 , α4 , α5 , α6 , 0, . . . , 0, α1 , α2 ) v21 := (α1 , α2 , α3 , α4 , α5 , α6 , 0, 0, . . . , 0) v31 := (0, (0, 0, α1 , α2 , α3 , α4 , α5 , α6 , 0, . . . , 0) .. . 1 vN/2 N/2 := (α5 , α6 , 0, . . . , 0, α1 , α2 , α3 , α4 )

w11 := (β 3 , β 4 , β 5 , β 6 , 0, . . . , 0, β 1 , β 2 ) w21 := (β 1 , β 2 , β 3 , β 4 , β 5 , β 6 , 0, 0, . . . , 0) w31 := (0, (0, 0, β 1 , β 2 , β 3 , β 4 , β 5 , β 6 , 0, . . . , 0) .. . 1 wN/2 N/2 := (β 5 , β 6 , 0, . . . , 0, β 1 , β 2 , β 3 , β 4 ).

Las relaciones entre los numeros wavelet y los scaling son como en db3, es decir, decir, β 1 := α := α6 , β 2 :=

−α5, β 3 := α := α4 , β 4 := −α3 , β 5 := α := α 2 , β 6 := −α1 .

3.2 Coiflets

33

Y también son iguales las siguientes igualdades que se cumplían en la wavelet db3, concretamente, α21 + α22 + α32 + α42 + α25 + α26 = 1,

√

α1 + α2 + α3 + α4 + α5 + α6 = 2, β 1 + β 2 + β 3 + β 4 + β 5 + β 6 = 0, 1β 1 + 2β 2 β 2 + 3β 3 β 3 + 4β 4 β 4 + 5β 5 β 5 + 6β 6 β 6 = 0, junto con las igualdades de ortogonalidad. Sin embargo ahora tenemos que la anulación del momento para términos cuadráticos no se cumple en las wavelets, pero a cambio obtenemos las siguientes anulaciones de momentos en las scaling 0 α3 + 1α 1 α4 + 2α 2 α5 + 3α 3 α6 = 0, −2α1 − 1α2 + 0α (−2)2 α1 + ( −1)2 α2 + 0 2 α3 + 1 2 α4 + 2 2 α5 + 3 2 α6 = 0. Estas dos nuevas ecuaciones nos permiten llegar a una mejor aproximación entre entre las tenden tendencia ciass y los valores alores origin originale aless de la señal. señal. Concre Concretam tamen ente, te, si tomamos una señal f señal f que que sea el muestreo con paso h paso h de de una analógica de clase C 3 , entonces a j1 = 2f 2 j + O(h3 ), j = 1, . . . , N /2 /2.

√

Las scaling y wavelets de orden superior se definen como cabría esperar, en función de las del nivel anterior mediante las fórmulas recurrentes m−1 v1m = α1 vN/2 N/2m−

m−1 + α3 v1m−1 + N/2 − −1 + α2 vN/2 α4 v2m−1 + α5 v3m−1 + α6 v4m−1 , v2m = α1 v1m−1 + α2 v2m−1 + α3 v3m−1 + α4 v m−1 + α5 v m−1 + α6 v m−1 , 1

4

1

m

5

6

.. . m m−1 vN/2 N/2m = α 1 vN/2 N/2m−

m−1 m−1 + α3 vN/2 + −3 + α2 vN/2 N/2 − −2 N/2 − −1 m−1 m−1 m−1 α4 vN/2 + α v + α v , 5 6 − 1 2 N/2 m−1 m−1 w1m = β 1 vN/2 + β 2 vN/2 + β 3 v1m−1 + N/2 − −1 N/2 − β 4 v2m−1 + β 5 v3m−1 + β 6 v4m−1 , wm = β 1 v m−1 + β 2 v m−1 + β 3 v m−1 + m

m

2

1

m

1

1

1

m

1

2

1

3

m

1

34


β 4 v4m−1 + β 5 v5m−1 + β 6 v6m−1 , .. . m−1 m wN/2 N/2m = β 1 vN/2 N/2m−

m−1 m−1 + β 3 vN/2 + −3 + β 2 vN/2 N/2 − −2 N/2 − −1 β 4 v m−1 − + β 5 v m−1 + β 6 v m−1 . N/2 N/2m

m

1

1

1

m

1

1

2

De esta manera las tendencias de orden superior aproximan la señal original a j2

≈ 2f 4 j ,

a j3

√

≈ 2

2f 8 j , . . .

De forma análoga se pueden definir coif2, coif3, coif4, coif5. La figura 3.2 nos muestra, para N = 214 y m = 11 11 varios varios vectores wavelet de coif3 y varios vectores scaling de coif4. Wavelets w11, w11 y w11 2

4

Scaling v11, v11 y v11

7

6

8

2

Figura 3.2: Para N = 214 y m = a la izquierda las wavelets w211 , w411 y w 711 m = 11 11 a de coif3 y a la derecha las scaling v611 , v811 y v211 de coif4.

3.3

Transfo ransforma rmada da de de Four Fourier ier disc discret reta a y MRA MRA

El análisis matemático del contenido en frecuencia de las señales es lo que se conoce como el análisis de Fourier. En esta sección trataremos las nociones básicas dentro del contexto de las señales discretas y su relación con el MRA. Mediante el estudio del contenido en frecuencia se puede profundizar más en el entendimiento de las wavelets. Para ello nos centraremos fundamentalmente en el análisis de las transformadas. Por conveniencia de notación en esta sección consideramos que el conjunto de índices de las señales es 0, 1, . . . , N 1 y que sus componentes se denotan

{

−}

3.3 Transformada de Fourier discreta y MRA

35

por f ( f ( j) j ) en vez de f j . Es decir f = (f (0) f (0),, f (1) f (1),, . . . f ( N

− 1)). − 1)).

El contenido en frecuencia de una señal discreta f se f se obtiene mediante la transformada discreta de Fourier (DFT). Dicha transformada la denotaremos por f y f ˆ y sus valores se definen como N 1

f ( f ˆ(k) :=

−



f ( f ( j) j )e−i2πkj/N .

j=0 j =0

Estos valores están definidos para cualquier valor k entero, no obstante nos restringiremos a k = 0, 1, . . . , N 1 por razones de periodicidad que ahora describimos. Recordamos, en primer lugar, que la exponencial con valores puramente imaginarios se define como

− −

eiθ := cos θ + i sen θ. Esto nos lleva por tanto a la periodicidad de la DFT f ( f ˆ(k + N ) N ) = f ( f ˆ(k) para todo entero k. Si nos quedamos con los N valores N valores de la DFT tomados anteriormente entonces ya tenemos descrita la transformada y así ˆ: CN

N

→C

,

 → fˆ,

f

teniendo eso si en cuenta la numeración de subíndices que hemos considerado en la DFT por razones de conveniencia. Aparte de la propiedad de periodicidad, esta transformación es lineal e invertible, siendo su fórmula de inversión 1 f ( j) j ) = N

N 1

−



f ( f ˆ(k )ei2πkj/N ,

k =0

para j = 0, . . . , N 1. Además se cumple la igualdad de Parseval

−

E (f ) f ) =

1 ˆ E (f ) f ). N

En la figura 3.3 tenemos las señales f ( f (x) = sen sen 4πx + 3 cos 24πx 2

g(x) = e −x (1 + sen sen 24 24πx πx))

(a) (c)

36


(a)

(b)

4

1000

2

500

0

0

−2

−500 −1000

−4 −10

0

10

−10

(c)

0

10

(d)

2

60

1.5 40

1 0.5

20

0

0

−0.5

−20 −10

0

10

−10

0

10

Figura 3.3: (a) Señal f . f . (b) Valor absoluto de DFT de la señal f . f . (c) Señal g . (d) Valor absoluto de DFT de la señal g .

muestreadas en el intervalo [ intervalo [ 16 16,, 16] 16] con con N N = = 1024 puntos, 1024 puntos, y sus DFT asociadas en valor absoluto ((b) y (d) respectivamente).

−

Como se puede observar en el primer caso hay unos “picos” en la DFT de la señal (a) en los puntos 2 y 12 12,, correspondie correspondient ntes es a los periodos de las funciofunciones en los sumandos, y en el segundo caso los tenemos en 12 12,, correspondiente al seno, y también en el 0.

± ±

±

La DFT de una señal se puede calcular de forma rápida mediante un algoritmo. Este se conoce como la transformada rápida de Fourier (FFT). Si realizamos de forma directa el cálculo de la DFT de una señal f f muestreada con N valores, N valores, entonces necesitamos hacer O(N ) N ) operaciones para cada valor f ( f ˆ(n). Teniendo en cuenta que la transformada también es una señal con N valores, obtenemos un total de O(N 2 ) operaciones. Con la FFT se reduce el número de operaciones a O(N log(N log(N )) )).. Para cada señal f f se define su espectro como la señal f ˆ 2 consistente en tomar los cuadrados de los valores absolutos de la transformada. Observamos

| |


37

que, según la identidad de Parseval, como N 1

−

|

1 E (f ) f ) = f ( f ˆ(n) 2 N n=0

|

entonces la energía de f resulta f resulta ser el valor medio del espectro de f . A partir del espectro podemos analizar el efecto que produce el tratamiento de señales mediante wavelets a nivel de frecuencia. En la figura 3.4 representamos los espectros de la primera scaling (a) y primera wavelet (b) a nivel 1 de coif2. (a)

(b)

3

3

2.5

2.5

2

2

1.5

1.5

1

1

0.5

0.5

0

0

−0.5

−0.5

−1 −0.5

0

0.5

−1 −0.5

0

0.5

Figura 3.4: (a) Espectro de la primera scaling coif2 a nivel 1. (b) Espectro de la primera wavelet coif2 a nivel 1. Como se puede apreciar, el espectro de la scaling se anula para frecuencias altas y toma su máximo valor en frecuencias bajas. En el caso de la wavelet es al contrario: su espectro alcanza el máximo en frecuencias altas y se anula en la zona media correspondiente a baja frecuencia. Para una mayor ilustración, tomamos la señal 2

g (x) = e −x (1 + sen sen 24 24πx πx)) muestreada en el intervalo [ intervalo [ 16 16,, 16] 16] con con N N = = 1024 puntos, 1024 puntos, la cual ya habíamos considerado con anterioridad. En la figura 3.5 representamos en (a) el primer promedio A1 del MRA mediante coif2. En (b) tenemos la representación de los valores absolutos de la DFT de A1 . El primer detalle D1 es el gráfico (c), y la representación de los valores absolutos de la DFT la tenemos en (d). Recordemos que f se f se descomponía como f = A 1 + D 1 . Al aplicar la DFT obtenemos un filtrado de frecuencias de f de f .. En la DFT de A de A1 han sido separadas

−

38


(a)

(b)

3

60

2

40

1

20

0

0

−1

−10

0

10

−20

−10

0

(c)

10

(d)

3

30

2 20

1 0

10

−1

0

−2 −3

−10

0

10

−10

−10

0

10

Figura 3.5: (a) Promedio A1 de g de g.. (b) Valores absolutos de DFT del promedio A1 . (c) Detalle D 1 de g . (d) Valores absolutos DFT del detalle D 1 .

las bajas frecuencias y en la DFT de D 1 las altas frecuencias. Esto es lo que se conoce como filtros de paso bajo y paso alto , respectivamente. Efectivamente, esto lo podemos verificar en general a partir de las fórmulas de las señales de promedio y detalle y por la linealidad de la transformada de Fourier discreta:

 

 

A1 = (f v11 )v11 +



1 1 (f · · · · · · + (f · vN/2 N/2 )vN/2 N/2 1  1 D1 = (f · (f · · w11)w11 + · · · + (f · wN/2 N/2 )wN/2 N/2 .

Si tenemos ahora en cuenta que el efecto de la transformada de Fourier sobre una traslación cíclica de una determinada señal (como es el caso de las scaling y wavelets al mismo nivel) tan sólo varía en multiplicar por elementos de módulo unitario las entradas de la transformada de la señal sin trasladar, llegamos a que en definitiva todas las wavelets al mismo nivel tienen idéntico espectro entre sí, como ocurre con las scaling. Por tanto las igualdades anteriores nos llevan a que A que A1 es una combinación de señales cuyo espectro se anula (o es muy pequeño) en las altas frecuencias, mientras que D1 es combinación de señales con espectro nulo en las frecuencias bajas.


39

Conjuntando esto con los espectros de scaling y wavelet que habíamos obtenido anteriormente deducimos que en el primer caso se tiene un filtro paso bajo de los valores de la DFT de f , f , que solo permite los valores de frecuencia baja, mientras que en el segundo caso representa un filtro paso alto de la DFT de f que f que solo permite los valores de frecuencia alta. A niveles superiores se cumplen las correspondientes fórmulas de descomposición, pero ahora hay que tener en cuenta que, mientras que v jm sigue siendo





nula en frecuencias altas, las transformadas w jm, m 2, se anulan para frecuencias bajas, pero también en las muy altas. Es lo que se conoce por filtro paso banda . Concretamente, se establecen un par de bandas de frecuencia que son las que fundamentalmente estamos filtrando. En la figura 3.6 se representan los espectros de scaling a distintos niveles (desde 1 arriba, hasta 4 abajo) a la izquierda, y los de wavelets (también desde 1 a 4) a la derecha. A partir del segundo espectro de las wavelets wavelets podemos po demos observar observar el fenómeno fenómeno de filtrado filtrado paso banda. (a)

(b)

1

1

0

0

−1

−1

−2

−2

−3

−3

−4

−4

−5

−5

−0.5

0

≥

0.5

−0.5

0

0.5

Figura 3.6: (a) Espectros de las scaling v scaling v1k de la familia coif2 para k para k = = 1, . . . , 4. (b) Espectros de las wavelets w1k de la familia coif2 para k = 1, . . . , 4. Los gráficos se han escalado con constantes 2 constantes 2 −k y se han desplazado verticalmente.

Profundizando algo más en el análisis de frecuencias, se pueden establecer

40


las tres siguientes igualdades (ver la sección 3.1 de [2]) en el primer nivel

| | | | | | | |    

vk1 ( j) j ) 2 + vk1 ( j + j + N/2) N/2) 2 = 2,

wk1 ( j) j ) 2 + wk1 ( j + j + N/2) N/2) 2 = 2,

vk1 ( j) j )wk1 ( j) j ) + vk1 ( j + j + N/2) N/2)wk1 ( j + j + N/2) N/2) = 0, 0,

k = 1, . . . , N /2 /2, j = 0, . . . , N /2 /2 1, donde z¯ denota el complejo conjugado de z (a + bi = bi = a a bi). bi). De las dos primeras igualdades deducimos que las bajas frecuencias son nulas (o casi nulas) si y solo si las altas tienen una gran magnitud, y viceversa. De la tercera tercera igualdad concluimos concluimos que una scaling scaling tiene nulas nulas o casi nulas altas frecuencias si y solo si las wavelets tienen bajas frecuencias nulas o de pequeña magnitud. Aquí hay que resaltar que las igualdades anteriores se verifican para cual1 1 1 quier base ortonormal v11 , . . . , v N/2 cada vk1+1 (respectiN/2 , w1 , . . . , w N/2 N/2 tal que cada v vamente, w vamente, w k1+1 ) se obtiene a partir del v k1 (resp., w (resp., w k1 ) desplazando cíclicamente sus entradas dos posiciones. Esto nos permite una forma sencilla de construcción de familias de wavelets ortogonales, ya que se trataría de fijar un sistema de vectores cumpliendo esas tres igualdades y aplicar a continuación la fórmula de inversión de la transformada de Fourier.

−

−

{

}

Capítulo 4

Compresión de señales y reducción de ruido En los capítulos anteriores hemos observado el fenómeno de la compactación de la energía al aplicar la transformada wavelet ya que esta se concentra fundamentalme mentalmente nte en la parte parte de la tendencia, tendencia, mientras mientras que las fluctuaciones fluctuaciones tienen tienen magnitud menor. Esto constituye, en gran parte, el punto clave en la compresión de señales y reducción del ruido. El proceso, esencialmente común en ambos casos, lo podríamos esquematizar de la siguiente manera:

Paso 1: Transformada. Realizar la transformada wavelet de la señal (contaminada o no) a un determinado nivel.

Paso 2: Cuantización. Los valores obtenidos en la señal transformada están en un determinado intervalo y deben ser almacenados con precisión finita, lo cual requiere de una “cuantización” de dichos valores pasando a estar en un conjunto finito fijo. Este paso puede combinarse con una previa umbralización, precisa en el caso de reducción del ruido y opcional para la compresión, que consiste en anular aquellos valores que se encuentren en un cierto intervalo alrededor del 0. Estos dos procesos conllevan pérdida de información.

Paso 3: Codificación. Se trata de codificar en modo binario la señal obtenida en el paso anterior de manera que ocupe lo mínimo posible. Aquí 41

42

Compresión de señales y reducción de ruido

no hay pérdida de información. Este paso no es necesario si tan sólo pretendemos la reducción del ruido. En el primer paso de la transformada, según las condiciones de suavidad de la señal, convendrá elegir las wavelets en consonancia con el suficiente número de momentos nulos. A continuación describiremos los pasos de cuantización y codificación con detalle.

4.1 4.1

Cuan Cu anti tiza zaci ción ón y umbr umbral aliz izad ado o

En el proceso de cuantización suponemos los valores de la señal transformada en un intervalo de la forma [ forma [ α, α]. Incluso podemos restarle a la señal su valor medio para que esté centrada en el 0. A continuación, si pretendemos cuantizar a q niveles q niveles (q (q entero entero positivo par), tomamos δ := 2α/q y y se define la función de cuantización escalar uniforme

−

Q(x) :=



0 sig( sig(x)(k )(k + 1/ 1 /2)δ 2)δ

si x δ si kδ < x

| | ≤

1) δ | | ≤ (k + 1)δ

Es decir, los valores que caen en cada subintervalo de la partición son redondeados al punto medio. Al aplicarle a la señal la función de cuantización Q(x), los valores los tomará ahora en el conjunto finito

− − q

q

1

q

− 3 α , . . . , 0, . . . , q − 3 α, q − 1 α α, − q q q



.

Si nuestro objetivo final es la compresión de la señal, la cuantización puede combinarse con una previa umbralización o corte de aquellos valores que no superen un cierto umbral. Para elegir el valor umbral hay que buscar un criterio. Una forma es basarnos en la energía. Concretamente, si ordenamos en forma decreciente los valores absolutos de la transformada de una señal f y1

≥ y2 ≥ · · · ≥ yN

entonces vamos calculando sucesivamente los cocientes y12 y 12 + y22 , ,... E (f ) E (f ) f )

4.1 Cuantización y umbralizado

43

hasta que lleguemos a un valor lo suficientemente cercano a 1 como para asegurar que captamos un porcentaje de energía adecuado. Por ejemplo podemos establecer el valor 0.9999 9999.. De esta manera al aplicar la transformada inversa posteriormente llegaremos a una aproximación de la señal original que contenga el 99 99..99 99 por por cien de la energía. Sea j el primer índice tal que y12 + y22 + + y j2 E (f ) f )

···

≥ 0.9999 9999..

Tomamos, por tanto, como valor umbral T := y := y j . En definitiva, una vez fijado el umbral T , T , se trata de anular los coeficientes que no superan el umbral. Esto se hace mediante una función de umbralización. El caso más simple consiste en anular los coeficientes por debajo del umbral y dejar como están los restantes, lo que se concreta en la llamada función de

umbralización fuerte Df (x) :=



0 x

si x T . si x > T

| | ≤ | |

Otro procedimiento alternativo, que tan sólo difiere en como afecta a los coeficientes que superan el umbral, es mediante la umbralización suave Ds (x) :=

 

0 x T x + T

−

si x T si x > T . si x < T

| | ≤ −

La umbralización es un paso preciso si nuestro objetivo es la reducción del ruido. Por ruido nos referiremos a la alteración no deseada de los valores de la señal original. Supondremos que dicho ruido es aditivo, es decir que g = f = f + x donde f donde f es es la señal original, x original, x es el ruido añadido, y g es la señal contaminada o alterada. Supondremos también que la señal de ruido x es x es aleatoria, también llamado ruido blanco, y de media 0. Al realizar la transformada wavelet de la señal g los coeficientes provenientes del ruido x, de altas frecuencias por su componente aleatoria, se concentran en la porción de las fluctuaciones. Un método simple de reducción del ruido consistiría en examinar visualmente la porción de fluctuaciones en la transformada para establecer el umbral de corte.

44


Sin embargo se puede llegar a un proceso más sofisticado basándose en métodos estadísticos. La compresión y la reducción del ruido blanco están estrechamente relacionados. Veremos como seleccionar el umbral de ruido en el caso (bastante usual) de que el ruido sea gaussiano. Por otra parte, para poder efectuar una reducción del ruido de forma adecuada hay que partir de la situación en que el umbral T necesario T necesario para dicha reducción es menor que el que se precisa para capturar la energía de la señal original (sin contaminar) de forma eficiente, en otro caso podríamos estar eliminando coeficientes significativos de la señal. Esto va a depender del tipo de wavelet que estemos utilizando. En el ejemplo de la figura 4.1 utilizamos la transformada db10. Se ha utilizado la señal generada con la función f ( f (t) =



t2 (2 t)2 sen8πt sen8πt 2(2 t)(3 t)2 sen20πt sen20πt

− −

−

si 0 si 2

≤t≤2 ≤t≤3

,

muestreada con N = N = 4096 puntos. 4096 puntos. A esta señal se le ha añadido ruido blanco gaussiano. (a)

(b)

2

2

1

1

0

0

−1

−1

−2

0

1000

2 00 0

30 0 0

4 0 00

−2

0

1000

(c) 2

1

1

0

0

−1

−1

0

1 0 00

2 00 0

30 0 0

4000

(d)

2

−2

2 0 00

3000

4 0 00

−2

0

10 0 0

2 0 00

3 0 00

4000

Figura 4.1: (a) Señal con ruido blanco. (b) Transformada db10 de 10 niveles con umbral de 0.227 227.. (c) Valores de la transformada que han sobrevivido al umbral. (d) Reconstrucción de la señal.

±

4.1 Cuantización y umbralizado

45

Si el ruido es gaussiano, es posible escoger el valor umbral T T de forma automática sin un conocimiento previo de la señal. En la figura 4.2 presentamos el ruido gaussiano (a), su histogramas de frecuencias de valores (b), la transformada db10 del ruido (c), y el histograma de dicha transformada (d). (b)

(a) 5

0.03 0.02

0

0.01 0

−5

0.5

1

1.5

2

−0.01

−100

−50

0

50

100

50

100

4

x 10

(c)

(d)

5

0.03 0.02 0.01

0

0 −5

0.5

1

1.5

2

−0.01

−100

−50

0

4

x 10

Figura Figura 4.2: 4.2: (a) Ruido Ruido blanco blanco gaussi gaussiano ano.. (b) Histog Histogram ramaa de frecue frecuencia nciass del ruido. (c) Transformada db10 del ruido. (d) Histograma de frecuencias de la transformada. Como podemos observar los histogramas presentan la típica forma de campana gaussiana. gaussiana. Una de las consecuenci consecuencias as de estar trabajando con transform transformaaciones ortogonales es que se preserva la naturaleza gaussiana del ruido cuando aplicamos la transformación. A partir de los histogramas deducimos que la media es 0 y la desviación típica es σ = 0.056 056.. La idea es tomar T tomar T de de forma que haya una alta probabilidad de estar justo por encima del máximo nivel de los coeficientes que provienen del ruido. Un resultado de Donoho y Johnstone (ver la sección 10.2.2 de [3] para más detalles) cuantifica el valor óptimo del umbral mediante la fórmula T = σ



2log N .

En el ejemplo de la figura 4.1 la fórmula nos da un valor de T = 0.227 227..

46


Nos queda la cuestión de como estimar la desviación típica σ típica σ del ruido que sea añadido a una cierta señal. La estrategia es la siguiente: Tomamos una porción de la transformada de la señal contaminada que consista fundamentalmente de valores provenientes del ruido (ese es el caso en general si tomamos una porción de la primera fluctuación, ya que los valores transformados de la señal sin contaminar que están en esa primera fluctuación deben tener una magnitud pequeña). Estimamos σ Estimamos σ a a partir de esa porción de señal seleccionada.

4.2 4.2

Codi Codific fica ación ción

En el paso de la cuantización, junto con la posible umbralización, gran parte de los coeficientes pasan a ser nulos. Esta información se almacena en lo que llamaremos mapa de significancia . Los valores serán 1 (si el correspondiente índice sobrevivió al corte) o 0 (si el índice correspondía a un valor menor que el umbral). Este mapa de significancia contiene muchos ceros, en general, con lo cual se comprime de forma eficiente con los métodos usuales de compresión sin pérdida de información. Un buen ratio de compresión lo alcanzaremos si tenemos en cuenta la entropía . Hay que hacer notar que no todos los valores de la señal transformada y cuantizada se repiten con la misma frecuencia. Esto lo podemos verificar mediante un histograma de frecuencias. En consecuencia, asignar el mismo número de bits a todos los valores es un gasto innecesario. Ciertamente podremos alcanzar una mayor compresión si asignamos menos bits a los valores que se repiten con mayor frecuencia, reservando más bits a los menos frecuentes. Este proceso se puede precisar de forma matemática utilizando resultados fundamentales en el campo de la teoría de la información . A continuación expondremos expondremos algunas algunas nociones nociones básicas. básicas. Supongamos que queremos codificar la siguiente secuencia de 16 números, del 1 al 4: 2412212332213232 (4.1) Esto lo podríamos hacer, por ejemplo, con 2 bits por número con las asignaciones 1 00,, 2 00 01,, 3 01 10,, 4 10 11,, resultando la secuencia de 32 11 bits: 01110001010001101001010010011001 (4.2)

→ 

→ 

→ 

→ 

Ahora bien, podemos observar que en la secuencia anterior de números hay algunos mucho más frecuentes (como el 2) que otros. Una forma de codificación

4.2 Codificación

47

mucho más inteligente consistiría en asignar menos bits a los valores más frecuentes y más bits a los menos frecuentes, pudiendo así obtener una secuencia de bits codificada más corta. Por ejemplo, en la secuencia de números anterior vemos que el 2 se repite 8 veces, el 3 lo hace 4 veces, el 1 tres veces y el 4 una vez. Si realizamos la siguiente asignación de bits 1 110,, 2 110 0, 3 10,, 10 4 111,, nos resulta la siguiente secuencia codificada de 28 bits: 111

→

→

→

0111110001100101000110100100

→

(4.3)

Esto produce un ahorro del 12 del 12..5 %. Aquí hay que hacer una observación fundamental mental en lo que respecta a sistemas sistemas de codificación codificación con asignación de longitud variable de bits: Deben ser decodificables es decir, no tiene que existir ambigüedad en la decodificación. ¿Como conseguimos este objetivo? Basta con que se cumpla la condición del prefijo: En una codificación de longitud variable ningún símbolo puede tener asignada una secuencia de bits que sea prefijo de la asignación de otro símbolo. Si nos fijamos en la codificación anterior, el 2 tiene asignado el 0, que no es prefijo de las secuencias del 1, 3 o 4, y el 3 tiene asignado 10, que no es prefijo de las secuencias del 1 o el 4, por tanto se cumple la condición del prefijo. Esta codificación del ejemplo se conoce como código Huffman, descubierta por David A. Huffman en 1951 y que hoy es utilizada en la transmisión comprimida de datos por internet, módem o fax! Esta codificación se puede sistematizar para secuencias arbitrarias de símbolos mediante árboles binarios. Si partimos de una determinada señal que, después del proceso de cuantización, toma sus valores en el conjunto finito vk , k = 1, . . . , n , y que pk kn=1 son las frecuencias relativas de repetición de la sucesión vk kn=1 . En este caso cada pk [0, [0, 1] y p1 + p + p2 + + pn = 1.

{

∈

{ }

}

{ }

···

Sea L Sea L k la longitud de la secuencia de bits que se le asigna al valor vk mediante un determinado sistema de codificación. La longitud media de bits asociada es L := p := p1 L1 + p + p2 L2 +

· · · + pnLn.

Por ejemplo, si tomamos la señal formada por la secuencia (4.1), mediante la codificación binaria de longitud fija (4.2), obtenemos L = L = 2. Pero si utilizamos la codificación de Huffman (4.3), llegamos a L = 1.75 75.. Un resultado de Shannon nos indica que L no puede ser inferior a la entropía . Este valor de entropía lo representaremos por y entonces tenemos

H

48


que L

≥ H := − p1 log2 p1 − p2 log2 p2 − · · · − pn log2 pn.

La cuestión es, en definitiva, como asignar las longitudes Lk para optimizar el valor L, teniendo en cuenta la restricción dada por el teorema de Shannon. Técnicas de codificación, como la codificación Huffman descrita, consiguen situar L tuar L entre entre y + 1. De hecho la aproximación es mejorable (con tendencia asintótica a ) si agrupamos los símbolos en bloques adecuadamente y posteriormente se aplica la codificación Huffman. Cabe mencionar que existen códigos correctores de errores que introducen algo de redundancia, pero que minimizan la transmisión de errores producidos por alteración de bits. Veamos un ejemplo en la figura 4.3. La señal representada en (a) consiste en la grabación de la palabra greasy , cuantizada a 8 a 8 bpp, tomando 16384 tomando 16384 puntos. puntos. La entropía es = 5.43 43.. Si codificamos directamente la señal con un buen método de codificación podemos esperar 5 esperar 5..93 93 bpp, bpp, que representan un total de unos 97000 unos 97000 bits bits al multiplicar por el número de puntos. Esto no nos conduce a una compresión muy efectiva teniendo en cuenta los 8 bpp 8 bpp originales. Mediante la codificación utilizando la transformada wavelet llegamos a una compresión importante. Tomamos la transformada coif5 a nivel 14 14 (gráfico (gráfico (c)). A la derecha se representan los histogramas de la señal original (b), y de la cuantización a 8 a 8 bits bits con umbralizado previo de los valores transformados (d). La entropía, en el segundo caso, es de 4 de 4..34 34.. Según el estimador establecido, un buen método de codificación de la señal transformada y cuantizada nos dará una longitud media de bits no superior a 4.84 84.. Ahora tenemos en cuenta que hay 3922 hay 3922 valores valores que superan el umbral en la transformada, lo que nos llevaría a un total de 3922 4.84 18922 bits. 18922 bits.

H H H

H

·

≈

Cuando se aplica la transformada inversa a los valores de la transformada cuantizados, la señal resultante es prácticamente indistinguible de la original. De hecho, todavía se puede mejorar la compresión. Recordemos que habíamos cuantizado con umbralizado previo toda la transformada de la señal anterior a 8 bpp, para después codificar en función de la entropía del histograma correspondiente. Si, como método alternativo, hubiéramos cuantizado la parte de la tendencia a 8 a 8 bpp, bpp, la de las fluctuaciones a 6 a 6 bpp bpp (lógico, pues el intervalo donde éstos valores se encuentran es claramente más pequeño que el de la tendencia), y además hubiéramos analizado por separado los histogramas de la tendencia y de las fluctuaciones a distintos niveles, observaríamos que la

4.3 Medidas de la aproximación

49

(a)

(b) 0.3

100 0.2

50

0.1

0 −50

0

−100 5000

10000

15000

−0.1

−100

−50

(c)

0

50

100

50

100

(d) 0.4

500

0.3 0.2 0 0.1 0 −500

5000

10000

15000

−0.1

−100

−50

0

Figura 4.3: (a) Señal original “greasy”. (b) Histograma de la cuantización de los niveles de intensidad de “greasy” a 8 bits. (c) Transformada coif5 de “greasy” de 14 niveles. (d) Histograma de la cuantización a 8 bits con umbralizado previo de los valores de la transformada

entropía en las fluctuaciones es menor, con lo cual podemos codificar por separado las subseñales subseñales de fluctuación fluctuación cuantizad cuantizadas as llegando llegando a una cantidad cantidad de bpp menor. Por ejemplo, con esta modificación aplicada a la transformada coif5 de nivel 4 de la señal anterior, llegamos a una compresión de 11305 11305 bits. bits.

4.3 4.3

Medi Me dida dass de la apro aproxi xima maci ción ón

Una vez completados los tres pasos del proceso anterior, mediante la transformada inversa obtenemos una aproximación de la señal original. En un buen número de aplicaciones es conveniente cuantificar el grado de aproximación entre señales. Existen distintas formas de medir este grado de aproximación, dos de las medidas más utilizadas son el error cuadrático medio (error RMS) y la relación señal-ruido (SNR). Dadas dos señales f y g , el error cuadrático medio entre f y g está dado

50


por RMS( RMS(f, g) :=

 − (f 1

g1

)2

+ (f (f 2

− g2

)2

+

N

· · · + (f (f N N −

gN )2

=

 √ − −

E (f g) , N

y la relación señal-ruido, medida en decibelios, por la expresión SNR( SNR(f, g) := 10 log10

 

N 2 i=1 f i

N i=1 (f i

− gi)2 .

Ambas medidas son muy utilizadas en el campo del tratamiento de la señal. El error RMS, debido a su sencillez, será la fórmula más habitual que utilizaremos para estimar el grado de aproximación entre señales.

Capítulo 5

Wavelet packets En esta sección estudiaremos transformadas ortogonales basadas en wavelets pero con una estructura en árbol diferente, de hecho lo que hacemos es completar el árbol binario. Este tipo de transformadas son las wavelet packets . En ellas seguiremos calculando las tendencias al igual que lo hacíamos con las wavelets, pero difiere en el hecho que no guardamos las fluctuaciones obtenidas en pasos anteriores sino que también las transformamos. A modo introductorio veremos como es el proceso con la transformada wavelet packet asociada a las transformada de Haar, la cual se conoce como transformada de Walsh . Sea H Sea H 1 la transformada de Haar de primer nivel. Al actuar sobre una señal f producimos producimos una subseñal de tendencia a tendencia a 1 y otra de fluctuaciones d fluctuaciones d 1 , es decir, H1 (f ) f ) = (a1 d1 ).

|

La transformada de Walsh de nivel 1, llamémosle Wh1 , actúa de la misma manera Wh1 (f ) f ) = (a1 d1 ) = H1 (f ) f ).

|

Para el siguiente nivel aplicaremos una transformada de Haar a la subseñal de tendencia y también a la de fluctuaciones: ahí radica la diferencia. Wh2 (f ) f ) = (H1 (a1 ) H1 (d1 )) = (a (a2 d2 H1 (d1 )). )).

|

| |

Si denotamos H1 (d1 ) = (d1,a d1,d ), donde d1,a es la subseñal de tendencia a primer nivel de la señal d señal d 1 y d 1,d es la subseñal de fluctuaciones a nivel 1 de la señal d1 , entonces Wh2 (f ) f ) = (a2 d2 d1,a d1,d ).

|

| | |

51

52

Wavelet packets

Cada una de las 4 subseñales anteriores tiene la misma longitud N/4 N/4. El proceso se repite para niveles superiores de manera que vamos calculando en cada paso las tendencias y fluctuaciones de cada una de las subseñales que tenemos en el paso anterior. Para el nivel 3 tendríamos Wh3 (f ) f ) = (H1 (a2 ) H1 (d2 ) H1 (d1,a ) H1 (d1,d ))

|

| | = (a3 |d3 |d2,a |d2,d |d1,a,a |d1,a,d |d1,d,a |d1,d,d ).

Comparando la transformada de Haar con la de Walsh, con esta última se suele conseguir una mayor compresión de la energía. La razón de ello reside en que al computar fluctuaciones de fluctuaciones podemos obtener más términos nulos o cercanos a cero. Matricialmente podemos expresar la transformada de Walsh de primer nivel al igual que la de Haar, y para niveles superiores se descompone en un producto de matrices ortogonales en bloques, los cuales son matrices de Haar. ¯ k) a la matriz ortogonal k Concretamente, si denotamos por H(k H( ortogonal k k asociada a la transformada de Haar H1 , entonces tenemos

×

 

.. . Wh j (f ) f ) .. .

 

=

 

¯ jN 0 ... H( ) 2 − ¯ jN 0 H( ) ... 2 − .. . 0 ...

0 0

1

1

0

0

... ...

0 ¯ jN H( ) 2 − 1



 

...

¯ N ) H( 0 2 ¯ N ) 0 H( 2

   ¯ N )) H(N H(

.. . f .. .

 

Como estamos realizando un producto de matrices ortogonales, se seguirá cumpliendo la ley de conservación de la energía. Por otra parte invertir el proceso consistirá consistirá simplemen simplemente te en trasponer trasponer matrices. matrices. El nivel máximo de transformada que podemos alcanzar, si N = 2n , es igual a n = log 2 N . N . Por tanto a lo sumo tenemos un producto de log2 N matrices. Cada una de ellas genera un número de operaciones no superior a 3N . N . En definitiva esto nos da un cómputo total de O(N log log N ) N ) operaciones para realizar la transformada de Walsh de un determinado nivel, al igual que con la transformada rápida de Fourier, pero superior al número de operaciones necesario con la transformada wavelet que era O(N ) N ).

53

El proceso que hemos descrito con la transformada de Walsh es análogo si tomamos una transformada wavelet distinta (por ejemplo dbI, coifJ). Posibles ventajas al aplicar la transformada wavelet packet, comparada con la transformada wavelet, se concretan en un ratio de compresión generalmente mejor. Cabría ponderar en cada caso las posibles ventajas de estas mejoras con el hecho de tener que realizar un número superior de operaciones. Concluimos diciendo que el método adoptado por el FBI para la compresión de huellas digitales, llamado método WSQ, es un híbrido entre la transformada wavelet y la wavelet packet.

Capítulo 6

La transformada wavelet continua En este capítulo analizaremos la transformada wavelet continua. En primer lugar el espacio de señales donde trabajaremos serán las señales definidas en toda la recta real de energía energía finita. finita. Concretam Concretament ente, e, señales señales f f L2 (R), es decir tales que E (f ) f ) :=

 |

∈ ∈

f (t) 2 dt <

|

R

∞.

Los valores que toman las señales pueden ser complejos (f ( f : R producto escalar definido en L2 (R) está dado por f g :=

· ·



f ( f (t)g ∗ (t)dt, f,g

R

→

∈ L2(R),

donde g ∗ (t) denota el conjugado de g (t). La transformada de Fourier para señales de módulo integrable es f ( f ˆ(w) :=



f ( f (t)e−iwt dt.

R

Esta transformada se extiende por completitud a L2 (R). El producto de convolución de f y g es f g(u) :=

∗ ∗



f ( f (t)g(u

R

55

− t)dt, u ∈ R.

C). El

56

La transformada wavelet continua

En lo que sigue supondremos que la señal scaling cumple que tiene media 1, es decir,

 

φ(t)dt = dt = 1,

R

y que la señal wavelet tiene media cero, es decir, ψ(t)dt = dt = 0.

R

Además también supondremos que ambas están normalizadas de forma que E (φ) = E (ψ) = 1. Por ejemplo, la wavelet de Haar se define como

 − 

1,

ψH (t) :=

0 t 1/2 1, 1/2 < t 1 . 0, t [0, [0, 1]

≤ ≤ ≤ ∈

La correspondiente función scaling es φH (t) :=

1, t 0, t

[0, 1] ∈ [0, [0, 1] ∈ [0,

.

Otro ejemplo interesante viene dado por la wavelet sombrero mejicano ψ (t) :=

  √ −

2 π 1/4

t2 σ2

1

3σ

e

−t 2 2σ 2

.

La figura 6.1 muestra su representación gráfica. A partir de la wavelet madre ψ obtenemos obtenemos las wavel wavelets ets “hijas” “hijas” mediant mediantee cambios de escala y traslaciones, es decir,

  − √

1 ψu,s (t) := ψ s

t

u

s

, u

∈ R, s > 0. 0 .

Con esta definición las hijas también están normalizadas y tienen media cero. La transformada wavelet continua de una señal f se f se define como W f ( f (u, s) := f := f ψu,s

· ·



1 ∗ = f ( f (t) ψ s R

√

− t

u

s

dt = dt = f f ψ¯s (u), u

∗ ∗

∈ R, s > 0, 0 ,

57

1

0.8

0.6

0.4

0.2

0

−0.2

−0.4 −8

−6 −6

−4

−2

0

2

4

6

8

Figura 6.1: Wavelet madre “sombrero mejicano”.

donde

1 ∗ ψ ψ¯s (t) := s

√

− 

t . s

En realidad la cantidad W f ( f (u, s) mide la variación de f en f en un entorno de u de magnitud proporcional a s. En la figura 6.2 tenemos una señal f ( f (t) y debajo la transformada wavelet mediante el sombrero mejicano. Ha de tenerse en cuenta que la notación en el entorno matemático MATLAB (usado para generar las figuras) es C a,b f (b, a) = W f ( f (posición, posición, escala) escala). a,b = W f ( Evidentemente se observa que hay una redundancia de información al computar la transformada wavelet continua puesto que pasamos de una señal unidimensional a una bidimensional. Bajo ciertas condiciones generales se puede establecer una fórmula de inversión. Para ello supondremos a partir de ahora que la wavelet tiende a cero en el infinito con suficiente rapidez, lo cual lo cuantificaremos mediante la expresión



(1 + t ) ψ (t) dt <

R

|||

|

∞.

En estas condiciones se cumple la siguiente fórmula de inversión para cualquier señal f L2 (R):

∈ ∈

1 f ( f (t) = C ψ

  [0, [0,

∞]

  − √

1 W f ( f (u, s) ψ s R

t

u

s

du

ds , t s2

∈ R.

58


1.5 1 0.5 0 −0.5 −1 −1.5

500

1000

1500

2000

2500

3000

3500

4000

Absolute Values of Ca,b Ca,b Coefficients for a = 512 445.7219 388.0234 337.794 337.794 294.0668 ..

a s e l a c s

2.297397 3.482202 5.278032 8 12.12573 18.37917 27.85762 42.22425 64 97.00586 147.0334 222.8609 337.794 512

500

1000

1500 2000 2500 time (or space) b

3000

3500

4000

Figura 6.2: Señal original y transformada wavelet continua con mexh.

La constante C ψ es la condición de admisibilidad dada por la expresión C ψ :=



|ψˆ(t)|2 dt, t

[0, [0,

∞]

la cual es finita en nuestro contexto de trabajo. También se cumple la siguiente ley de conservación de la energía:

 | R

1 f ( f (t) 2 dt = dt = C ψ

|

  | [0, [0,

∞]

W f ( f (u, s) 2 du

|

R

ds . s2

Por otra parte, a partir de la scaling padre podemos definir las scaling “hijos” correspondientes mediante

  − √

1 t u φu,s (t) := φ s s

, u

0 . ∈ R, s > 0.

La señal de aproximación de f a f a escala s es entonces Lf ( Lf (u, s) := f := f φu,s , u

· ·

∈ R.

59

A partir de la fórmula de inversión anterior se puede probar que, para una escala fija s0 , obtenemos que 1 f ( f (t) = C ψ 1 s0 C ψ

   [0,s [0,s0 ]

  − √ − ∈

1 W f ( f (u, s) ψ s R

Lf (u, s0 )

R

√ 1s0 φ

t

t

u

s

u

s0

du, t

du

ds + s2 R.

En la expresión anterior la segunda integral correspondería a la aproximación y la primera a los detalles. Con la transformada wavelet continua estamos también realizando un filtrado de la señal. Para ello tenemos en cuenta la fórmula (  f g)(w )(w) = f ˆ(w)ˆg(w).

∗ ∗

Fijando una escala s > 0, 0 , entonces  W f ( f (w) = (f ψ¯s )(w )(w) = f ( f ˆ(w)ψ¯s (w).





∗ ∗

Si consideramos, por ejemplo, la wavelet sombrero mejicano, la cual es simétrica, entonces ψˆ¯s es real, y la igualdad anterior se convierte en un producto de transformadas de Fourier. Multiplicar por la transformada de la wavelet sombrero mejicano a distintas escalas significa realizar un filtrado de frecuencias sobre la señal. En la figura 6.3 damos los gráficos de espectros del sombrero mejicano mejicano para diversas diversas escalas. escalas. 6

5

4

3

2

1

0 −8

−6

−4

−2

0

2

4

6

8

Figura 6.3: Espectros de la wavelet sombrero mejicano a distintas escalas.

60


Al multiplicar los anteriores espectros por f estamos f ˆ estamos dando una descomposición de esta última en distintas bandas de frecuencia. Las bandas correspondientes a frecuencias altas son las asociadas a escalas pequeñas. Es decir, hay una relación inversa entre los valores de escala y los de frecuencia. Veamos un ejemplo más que analizaremos con la wavelet de Morlet. La wavelet de Morlet madre se define por := e ψ(t) := e

−t 2 2

cos(5t cos(5t).

La figura 6.4 muestra su representación gráfica. 1 0.8 0.6 0.4 0.2 0 −0.2 −0.4 −0.6 −0.8 −1 −8

−6

−4

−2

0

2

4

6

8

Figura 6.4: Wavelet madre Morlet. Tomemos la señal discreta f generada f generada al muestrear la señal (2 sen( sen(a a1 x) + cos(a cos(a2 x))e ))e−b

1

(x c1 )2

−

+

sen(a sen(a3 x)e−b

3

(x c 3 ) 2

−

+

(sen(a (sen(a1 x) + 2 co cos( s(a a2 x))e ))e−b

2

(x c 2 ) 2

−

en el intervalo [0, [0, 1] con N N = 2048 analógica, donde a1 = 40 40π π , a2 = 80 80π π, a3 = 60 60π π, b1 = 80 80π π , b2 = 40 40π π , b3 = 50 50π π, c1 = 0.25 25,, c2 = 0.75 y c3 = 0.5. La figura 6.5 muestra la gráfica de la señal y debajo el escalograma (transformada wavelet continua) de dicha señal tomando la wavelet de Morlet. Podemos observar como se detectan las bandas de frecuencia que corresponden a a a a 1 y a 2 en la primera y tercera subseñales. Las bandas de frecuencia que corresponden a a3 se detectan en la segunda subseñal. Es decir, en definitiva el escalograma provee un retrato tiempo-frecuencia de la señal.

61

4 2 0 −2 −4 0

0.2

0.4

0.6

0.8

1

Absolute Values of Ca,b Ca,b Coefficients Coefficients for a = 1 1.5 2 2.5 3 ...

a s e l a c s

128 118 108 98 88 78 68 61 56 51 46 41 36 31 26 21 16 11 6 1

200

400

600

800 1000 1200 time (or space) b

1400

1600

1800

Figura 6.5: Señal original y su escalograma con morl.

2000

Apéndice A

MATLAB y la “Wavelet Toolbox” La “Wav “Wavelet elet Toolbox” oolbox” es una herram herramien ienta ta para para el tratam tratamien iento to de señales señales e imágenes mediante wavelets. wavelets. Se ejecuta dentro del entorno MATLAB MATLAB por lo que es necesario tener instalado previamente este programa (así como la toolbox) en nuestro ordenador. Aquellos lectores que deseen aspectos básicos sobre el uso de MATLAB pueden consultar su manual [4]. Para el uso de la herramienta específica de wavelets recomendamos también el manual oficial [5].

A.1

MATLAB

No pretendemos hacer aquí una descripción de los comandos y las posibilidades que tiene el programa MATLAB, tan solo vamos a comentar aquellas órdenes que serán utilizadas para poder analizar señales que pueden estar en formato MATLAB, o puede ser necesario manipularlas para que MATLAB trabaje con ellas .

A.1.1 A.1.1

Brev Breve (muy (muy brev breve) e) int introdu roducci cción ón

El programa MATLAB utiliza aritmética de doble precisión (“double”) para todas sus operaciones, es decir, todos los números son representados internamente usando 8 bytes. Esto nos obliga a preparar nuestros datos antes de poder analizarlos en el entorno gráfico. 63

64

MATLAB y la “Wavelet Toolbox”

Cuando se comienza el programa, éste crea un espacio de trabajo llamado “work space”. Para saber qué variables tenemos en este espacio teclee who

Si además queremos información del tamaño de estas variables y del tipo de datos que contienen podemos teclear whos

Suponiendo que tengamos una variable a que contiene una matriz de n de n filas filas por m columnas, si escribimos a(i:j,k:l)

el programa extrae la submatriz con filas desde la i a la j la j y columnas desde la k a la l. Para los vectores fila tan sólo hay que especificar la columna y para los vectores columna sólo la fila. Si deseamos guardar una sesión para después continuar con nuestro traba jo exactamente en el mismo punto debemos usar los comandos save y load seguido del nombre del archivo. Este archivo estará en formato MATLAB y su extensión será *.mat Si olvidamos cómo es la sintaxis de una orden siempre se puede recurrir a la ayuda “on line” del programa escribiendo help help comand comando o

Alternativamente se puede usar la ayuda picando con el ratón en el botón “Help” de la barra de menús del programa.

A.1.2 A.1.2

Repre Represen sentac tación ión de señale señaless

Todas las señales que analicemos estarán representadas como vectores fila o columna. Estas señales pueden proceder de aparatos de medida, representación de funciones, etc. Naturalmente antes de analizar la señal hay que preparar los datos y cambiarlos a un formato que sea apropiado para MATLAB. A continuación veremos algunos casos concretos que suelen ser habituales.

A.1 MATLAB

65

Datos en un fichero ASCII Supongamos que nos llegan los valores numéricos en ASCII de una señal que deseamos tratar en un fichero llamado sig1.txt, es decir, el fichero sig1.txt tiene 1024 datos ordenados en una columna y cada fila corresponde a un dato. En la dirección URL http://www.upv.es/frechet/wavelets pueden encontrarse los materiales del curso, entre otros, el fichero sig1.txt Debe guardarse en la carpeta C:\Mis C:\Mis Documentos Documentos. Para cargarlo en MATLAB usaremos la orden load ’C:\Mis ’C:\Mis Documentos\s Documentos\sig1.tx ig1.txt’ t’ -ascii -ascii

Si inspeccionamos ahora el espacio de trabajo de MATLAB veremos que ha creado una variable llamada sig1 con los valores del fichero. Además es un vector columna y el formato de cada valor es de doble precisión. Si lo deseamos ahora podemos duplicar esa variable. Por ejemplo para duplicarla con nombre nsig1 escribiremos nsig1 = sig1;

El punto y coma al final de la orden es importante para que el programa suprima la respuesta de la evaluación de la orden y no escriba todos los datos que representa la variable. A continuación salvaremos la variable en un fichero con formato MATLAB escribiendo save ’C:\Mis ’C:\Mis Documentos\n Documentos\nsig1’ sig1’ nsig1

Observación A.1 Para no tener que especificar siempre la ruta para llegar al fichero (tanto para leer como para grabar) podemos añadir esa ruta al path del programa MATLAB. Esto se puede hacer seleccionando la ruta deseada en el campo Current Directory. Si no aparece en la lista puede buscarse picando en el cuadrado con tres puntos que hay a la derecha. Seleccionar la carpeta C:\Mis C:\Mis Documentos Documentos. Ahora tenemos un fichero nsig1.mat que puede ser cargado en el entorno gráfico de la “wavelets toolbox” para el análisis de la señal. También podemos guardar los datos en otro fichero en formato ASCII. save save ’nsig1 ’nsig1.tx .txt’ t’ nsig1 nsig1 -ascii -ascii

Si los datos de la señal vienen en un fichero pero escritos en una fila, el procedimiento es el mismo aunque hay que trasponer la matriz para convertirla en un vector columna. La orden para trasponer una matriz es

66


matriz’;

Una vez los datos son un vector fila o columna, representarlos gráficamente es sencillo con la orden plot. plot(nsig1); (véase la figura A.1)

1 0.8 0.6 0.4 0.2 0 −0.2 −0.4 −0.6 −0.8 −1 0

200

400

600

800

1000

1200

Figura A.1: Un ejemplo de plot().

Datos creados a partir de una función Algunas Algunas veces puede ser conven convenien iente te “fabricar” “fabricar” la señal a partir partir de una función dada (para docencia, pruebas, ejemplos, simulaciones, etc). Para ello evaluaremos la función en ciertos puntos y con el resultado construiremos la señal. Veamos como se hace con un ejemplo práctico.

Ejemplo A.2 Generar una señal con 1024 datos evaluados con la función g(x) = 20x 20x2 (1

− x)4 cos(12πx cos(12πx))

en el intervalo [0, [0, 1]. 1]. La forma trivial de hacerlo es con una estructura bucle for. Cuando el número de evaluaciones es grande se recomienda construir el vector usando un método llamado “vectorization”, el cuál es mucho más efectivo y más rápido. Ilustraremos el segundo método por ser el más estándar. También hay que notar que este método crea señales que son vectores fila. Pueden transponerse si se desea para que sean vectores columna.

A.1 MATLAB

67

x = 0:1/10 0:1/1023: 23:1; 1; signal = (20*x.^2).* (20*x.^2).*((1-x) ((1-x).^4).* .^4).*cos(1 cos(12*pi*x 2*pi*x); );

La diferencia entre el operador multiplicación * y el “punto multiplicación” .* es importante. Cuando MATLAB lee * hace una multiplicación usual de matrices (incluido el caso de escalar por matriz); cuando utilizamos .* el programa realiza un producto de matrices o vectores componente a componente. Esta misma regla se aplica a las potencias ^ y .^. La representación gráfica de los vectores se puede hacer con el comando Y que queremos representar se fija con la orden plot. La porción del plano X Y que axis(). plot(signal); axis([ axis([0 0 1024 1024 -0.5 -0.5 0.5]); 0.5]); (véase la figura A.2)

0.5 0.4 0.3 0.2 0.1 0 −0.1 −0.2 −0.3 −0.4 −0.5

0

200

400

600

800

1000

Figura A.2: Señal creada a partir de la función g (x). En MATLAB es usual escribir una orden por línea aunque si deseamos escribir varios comandos en la misma línea se deben separar con ;. Por ejemplo, los comandos para dibujar la figura A.2 también se pueden escribir plot(s plot(sign ignal) al); ; axis([ axis([0 0 1024 1024 -0.5 -0.5 0.5]); 0.5]);

Finalmente hay que mencionar que para eliminar las variables definidas se debe usar el comando clear clear var, donde var es el nombre de la variable que se desea eliminar. Cuando se quiere eliminar todas las variables definidas en el entorno de trabajo se usa clear sin especificar ninguna variable. variable. Es conveniente conveniente usarlo antes de empezar un ejercicio que es independiente de los anteriores. De esta manera nos evitaremos interferencias con variable usadas previamente.

68

A.2 A.2


La “Wa “Wavelet elet Toolbo oolbox” x”

La “Wav “Wavelet elet Toolbox” oolbox” es una una herram herramien ienta ta para para el tratam tratamien iento to de señales señales e imágenes mediante wavelets. wavelets. Se ejecuta dentro del entorno MATLAB MATLAB por lo que es necesario tener instalado previamente este programa (así como la toolbox) en nuestro ordenador.

A.2.1 A.2.1

Wavele avelets ts 1-D desde desde la línea línea de comandos comandos (prim (primer eros os pasos)

Supondremos que la señal está almacenada en un fichero sig.mat en un vector columna. Obtener el fichero sig.mat de la página web del curso y guardarlo en la carpeta Mis Documentos Documentos o crear la variable a partir de la función g (x) = 20x 20x2 (1

− x)4 cos(12πx cos(12πx))

muestreada con 2048 puntos en el intervalo [0, [0, 1]. 1]. clear; load load ’sig’; ’sig’; s = sig; ls = length length(s) (s); ;

Para calcular la primera tendencia y la primera fluctuación deberemos efectuar los productos escalares con las wavelets y las scaling de nivel 1. [C,L] [C,L] = wavede wavedec(s c(s,1, ,1,’ha ’haar’ ar’)[a )[a b];

La función wavedec calcula los productos escalares de la señal con las wavelets y las scaling de nivel 1 y almacena los valores de la tendencia y la fluctuación en un solo vector llamado C; el vector L contiene las longitudes. La orden necesita tres argumentos, el primero es la señal, el segundo los niveles de descomposición y el tercero es la wavelet a utilizar. Ahora debemos extraer los coeficientes de la tendencia y la fluctuación para formar las subseñales. cA1 = appcoef(C,L, appcoef(C,L,’haar’ ’haar’,1); ,1); cD1 = detcoe detcoef(C f(C,L, ,L,1); 1);

Para representar gráficamente los valores se usa la orden plot que ya mencionamos en A.1.2. Por ejemplo, para representar la tendencia y la fluctuación escribiremos

A.2 La “Wavelet Toolbox”

69

plot(cA1); axis([ axis([1 1 length length(cA (cA1) 1) -1 1]); 1]); plot(cD1); axis([ axis([1 1 length length(cD (cD1) 1) -0.25 -0.25 0.25]) 0.25]); ; (véase la figura A.3) 1

0.25

0.8

0.2

0.6

0.15

0.4

0.1

0.2

0.05

0

0

−0.2

−0.05

−0.4

−0.1

−0.6

−0.15

−0.8

−0.2

−1

2 00

400

600

800

1000

−0.25

200

400

600

8 00

1 000

Figura A.3: Tendencia y fluctuación de la señal sig.mat. Aunque en la figura A.3 aparecen las dos gráficas una al lado de la otra, en realidad MATLAB las muestra de una en una en una misma ventana que MATLAB abre para mostrar el resultado de la orden plot. Para que sean representadas en una misma ventana se usa subplot como se muestra a continuación. subplot(1,2,1); plot(cA1); title(’Tendencia’); axis([ axis([1 1 length length(cA (cA1) 1) -1 1]); 1]); subplot(1,2,2); plot(cD1); title(’Fluctuación’); axis([ axis([1 1 length length(cD (cD1) 1) -0.25 -0.25 0.25]) 0.25]); ;

Para construir la primera señal promedio y la primera señal de detalle del análisis de multirresolución debemos usar estos valores como coeficientes en las bases de wavelets y de scaling. Para esta tarea se dispone del comando wrcoef. Acepta cinco argumentos: el primero indica si queremos aproximaciones o detalles, el segundo y tercero son los coeficientes y las longitudes de las subseñales que produjo la orden wavedec, el cuarto es la wavelet a utilizar y por

70


último el nivel de aproximación o detalle que deseamos calcular. Finalmente representaremos las dos señales en una única ventana. A1 = wrcoef(’a’, wrcoef(’a’,C,L,’h C,L,’haar’,1 aar’,1); ); D1 = wrcoef(’d’, wrcoef(’d’,C,L,’h C,L,’haar’,1 aar’,1); ); subplot(1,2,1); plot(A1); title(’Promedio’); axis([1 ls -1 1]); subplot(1,2,2); plot(D1); title(’Detalle’); axis([ axis([1 1 ls -0.1 -0.1 0.1]); 0.1]); (véase la figura A.4)

Promedio

Detalle

1

0.1

0.8

0.08

0.6

0.06

0.4

0.04

0.2

0.02

0

0

−0.2

−0.02

−0.4

−0.04

−0.6

−0.06

−0.8

−0.08

−1

500

1000

1500

2000

−0.1

500

1000

1 500

20 00

Figura A.4: Reconstrucción de las señales promedio y detalle de primer nivel de la señal sig.mat.

Observación A.3 Cuando el número de comandos es grande y se va a utilizar repetidamente es conveniente hacer uso de la herramienta que permite guardar una serie de comandos en un archivo M-file. Estos archivos tienen extensión *.m y contienen las órdenes que escribiríamos en la línea de comandos de MATLAB. Para ejecutar un archivo de comandos de MATLAB debemos escribir el nombre del archivo (no es necesario escribir la extensión) en la línea de comandos del programa. En este tipo de archivos es conveniente incluir comentarios que ayuden a recordar las acciones de algunos comandos, de esta manera será mucho más sencillo modificarlos. Los comentarios siempre deben comenzar por


71

% comentario comentario... .... El programa no evaluará el texto que se escriba desde el %

hasta el final de la línea. Si la orden que pretendemos escribir es demasiado larga y deseamos partirla en varias líneas, esto se puede hacer escribiendo ... en el punto de la línea donde deseemos cortar la orden. También se puede hacer uso de herramientas de depuración cuando nuestras órdenes no realizan la tarea deseada aunque no entraremos en la explicación de estas herramientas. El uso fundamental que nosotros daremos a estos ficheros es la posibilidad de escribir un conjunto de órdenes de MATLAB en un proceso automatizado que puede ser usado repetidas veces con posiblemente pequeñas modificaciones. Finalmente una observación respecto de los ficheros M-file que vamos a incluir en este texto. Como ya hemos mencionado es habitual escribir una orden por línea e incluir comentarios para hacer el código más comprensible. Nosotros escribiremos en algunos casos varias órdenes por línea con la idea de que los códigos no ocupen tantas líneas y en la impresión quede excesivo espacio “en blanco”. Las subseñales tendencia (señales promedio) y subseñales fluctuación (señales detalle) a niveles superiores se obtienen de manera similar. Veamos el ejemplo anterior llegando a nivel 3 y haremos uso del archivo de comandos ana_nivel_3.m.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% File File ana_ni ana_nivel vel_3. _3.m m %% Análisis wavelet de una señal hasta nivel 3 %% La variable s contiene la señal a analizar w = ’haar’; % Fijamos la wavelet n = 3; % Niveles [C,L] [C,L] = wavede wavedec(s c(s,n, ,n,w); w); % Descom Descompos posici ición ón de la señal señal % Extraemos los valores de la tendencia y de % las fluctuaciones del vector C. cA1 = appcoef(C,L appcoef(C,L,w,1); ,w,1); cA2 = appcoef(C,L appcoef(C,L,w,2); ,w,2); cA3 = appcoef(C,L appcoef(C,L,w,3); ,w,3); cD1 = detcoe detcoef(C f(C,L, ,L,1); 1); cD2 = detcoe detcoef(C f(C,L, ,L,2); 2); cD3 = detcoe detcoef(C f(C,L, ,L,3); 3);

72


% Representam Representamos os estos coeficientes coeficientes. . figure; subplot(3,2, subplot(3,2,1); 1); plot(cA1); plot(cA1); title(’Tend title(’Tendencia encia a axis([ axis([1 1 length length(s) (s)/2 /2 -1.5 -1.5 1.5]); 1.5]); subplot(3,2, subplot(3,2,2); 2); plot(cD1); plot(cD1); title(’Fluc title(’Fluctuació tuación n axis([ axis([1 1 length length(s) (s)/2 /2 -0.03 -0.03 0.03]) 0.03]); ; subplot(3,2, subplot(3,2,3); 3); plot(cA2); plot(cA2); title(’Tend title(’Tendencia encia a axis([ axis([1 1 length length(s) (s)/2 /2 -1.5 -1.5 1.5]); 1.5]); subplot(3,2, subplot(3,2,4); 4); plot(cD2); plot(cD2); title(’Fluc title(’Fluctuació tuación n axis([ axis([1 1 length length(s) (s)/2 /2 -0.03 -0.03 0.03]) 0.03]); ; subplot(3,2, subplot(3,2,5); 5); plot(cA3); plot(cA3); title(’Tend title(’Tendencia encia a axis([ axis([1 1 length length(s) (s)/2 /2 -1.5 -1.5 1.5]); 1.5]); subplot(3,2, subplot(3,2,6); 6); plot(cD3); plot(cD3); title(’Fluc title(’Fluctuació tuación n axis([ axis([1 1 length length(s) (s)/2 /2 -0.03 -0.03 0.03]) 0.03]); ;

nivel 1’); a nivel 1’); nivel 2’); a nivel 2’); nivel 3’); a nivel 3’);

% Recons Reconstru truimo imos s a contin continuac uación ión las señale señales s promed promedio io y detall detalles es A1 = wrcoef(’a’, wrcoef(’a’,C,L,w, C,L,w,1); 1); A2 = wrcoef(’a’, wrcoef(’a’,C,L,w, C,L,w,2); 2); A3 = wrcoef(’a’, wrcoef(’a’,C,L,w, C,L,w,3); 3); D1 = wrcoef(’d’, wrcoef(’d’,C,L,w, C,L,w,1); 1); D2 = wrcoef(’d’, wrcoef(’d’,C,L,w, C,L,w,2); 2); D3 = wrcoef(’d’, wrcoef(’d’,C,L,w, C,L,w,3); 3); % Repres Represent entamo amos s estas estas señale señales s en una nueva nueva ventan ventana. a. subplo subplot(3 t(3,2, ,2,1); 1); plot(A plot(A1); 1); title( title(’Pr ’Prome omedio dio a nivel nivel 1’); 1’); axis([ axis([1 1 length length(s) (s) -0.5 -0.5 0.5]); 0.5]); subplo subplot(3 t(3,2, ,2,2); 2); plot(D plot(D1); 1); title( title(’De ’Detal talle le a nivel nivel 1’); 1’); axis([ axis([1 1 length length(s) (s) -0.02 -0.02 0.02]) 0.02]); ; subplo subplot(3 t(3,2, ,2,3); 3); plot(A plot(A2); 2); title( title(’Pr ’Prome omedio dio a nivel nivel 2’); 2’); axis([ axis([1 1 length length(s) (s) -0.5 -0.5 0.5]); 0.5]); subplo subplot(3 t(3,2, ,2,4); 4); plot(D plot(D2); 2); title( title(’De ’Detal talle le a nivel nivel 2’); 2’); axis([ axis([1 1 length length(s) (s) -0.02 -0.02 0.02]) 0.02]); ; subplo subplot(3 t(3,2, ,2,5); 5); plot(A plot(A3); 3); title( title(’Pr ’Prome omedio dio a nivel nivel 3’); 3’); axis([ axis([1 1 length length(s) (s) -0.5 -0.5 0.5]); 0.5]); subplo subplot(3 t(3,2, ,2,6); 6); plot(D plot(D3); 3); title( title(’De ’Detal talle le a nivel nivel 3’); 3’); axis([ axis([1 1 length length(s) (s) -0.02 -0.02 0.02]) 0.02]); ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Cuando ejecutamos este fichero obtendremos dos ventanas con los gráficos de la figura A.5.


Tendencia a nivel 1

73

Fluctuacion a nivel 1

Promedio a nivel 1 0.5

1

Detalle a nivel 1 0.02

0.02 0.01

0

0

−1

−0.02

0

0 −0.01

−0.5 −0.02 200400600 8001000 500 10001500 100015002000 2000 Fluctuacion a nivel 2 Promedio a nivel 2

200 400600 8001000 Tendencia a nivel 2

0.5 1

500 10001500 100015002000 2000 Detalle a nivel 2

0.02

0.02 0.01

0

0

−1

−0.02

0

0 −0.01

−0.5 −0.02 200 400600 8001000 500 10001500 100015002000 2000 Fluctuacion a nivel 3 Promedio a nivel 3

200 400600 8001000 Tendencia a nivel 3

0.5 1

500 10001500 100015002000 2000 Detalle a nivel 3

0.02

0.02 0.01

0

0

−1

−0.02

0

0 −0.01

200 400600 8001000

−0.5 200 400600 8001000

−0.02 500 100015002000

500 100015002000

Figura A.5: Análisis de tres niveles de la señal sig.mat. Las dos columnas de la izquierda representan las subseñales tendencia y fluctuación y las dos de la derecha las señales promedio y detalles.

Ejercicio A.4 Dadas las funciones 1. f ( f (x) = 1 2. f ( f (x) = x 2 (1

− x) 3. f ( f (x) = x 4 (1 − x)6 cos(64πx cos(64πx)) 4. f ( f (x) = 1(0. 1(0.2 < x < 0. 0 .3) − 3(0. 3(0.4 < x < 0. 0 .5) + 2(0. 2(0.5 < x < 0. 0 .8) 5. f ( f (x) = sgn(sin(12πx sgn(sin(12πx)) ))

74


realizar con cada una de ellas las siguientes tareas:

• Crear señales de longitud 2048 que representen cada función en el intervalo [0, [0, 1]. 1].

• Obtener las subseñales tendencia y fluctuación utilizando la wavelet de Haar hasta nivel cuatro. Imprimir una hoja con las ocho gráficas.

• En la señal fluctuación de nivel 1 (representado por cD1 en la notación de los ejemplos anteriores), obtener un vector que tenga un 1 en la posición iésima si la correspondiente i-ésima coordenada de cD1 es no nula (soporte o mapa de “significancia” de la subseñal cD1).

• Para cada función contar cuantos unos hemos obtenido en el paso an-

terior. Para hacer esto basta con sumar las coordenadas del soporte o mapa de significancia.

• Ordenar las funciones por número de elementos no nulos en la subseñal fluctuación de nivel 1. Relacionar el orden obtenido con la naturaleza de las señales analizadas.

Para realizar este ejercicio se puede tomar el fichero ana_nivel_3.m y modificarlo adecuadamente.

A.2.2 A.2.2

Wavele avelets ts 1-D desde desde el entorn entorno o gráfico gráfico (primero (primeross pasos) pasos)

La “Wavelet Toolbox” también tiene un entorno gráfico que resulta potente y sencillo de usar. Debe hacerse notar que el método explicado anteriormente ofrece más control sobre el análisis de una señal mediante wavelets que hacer el análisis en el entorno gráfico. En el otro lado de la balanza tenemos que el entorno gráfico es más sencillo de utilizar. Por ejemplo, desde la herramienta gráfica no es posible visualizar las subseñales tendencia y fluctuación. Vamos ahora a repetir los ejercicios anteriores desde el entorno gráfico. Es importante notar que las señales a analizar deben estar disponibles en archivos *.mat antes de poder analizarlas desde la herramienta gráfica. Tanto si estamos “creando” una señal a partir de una función o si nos llega en fichero *.txt u otro formato, los datos deben ser tratados adecuadamente hasta obtener una señal fila o columna que será guardada en un archivo *.mat. Este archivo podrá ser leído posteriormente desde el entorno gráfico. En primer lugar debemos cargar la herramienta gráfica.

A.3 Programando algunas funciones en MATLAB

75

wavemenu

A continuación se abrirá una ventana para seleccionar la acción que deseamos realizar. A partir de este momento el entorno gráfico se usa de la manera habitual con el ratón. Por este motivo no nos excederemos de todas y cada una de las opciones en esta breve introducción. 1. Seleccionam Seleccionamos os Wavelet 1D. 2. Cargar Cargar la señal grabada grabada en el fichero fichero sig.mat con Load—signal. 3. Seleccionam Seleccionamos os la wave wavelet let a utilizar utilizar (haar) y los niveles niveles (3). 4. Seleccionam Seleccionamos os Analyze. 5. Usar los distintas distintas maneras de visualizar visualizar las señales promedio y los detalles (en la figura A.6 tenemos el análisis con método Full decomposition de visualización). 6. Seleccionar Seleccionar Histograms y ver los histogramas de frecuencias de las distintas señales o subseñales. Notar que en este caso se pueden obtener histogramas tanto de la señal aproximación (señal detalle) como de la subseñal tendencia (fluctuación).

Ejercicio A.5 Con la herramienta gráfica de la “Wavelet Toolbox” tomar las señales creadas en el ejercicio A.4 y hacer un análisis de las mismas, experimentando las distintas opciones del entorno gráfico.

A.3

Progra Programan mando do alguna algunass func funcion iones es en MA MATLA TLAB B

Cuando se pretende hacer una tarea repetidas veces es conveniente definir una función en MATLAB que realice esta tarea y que pueda ser invocada de una manera sencilla y rápida. Es conveniente seguir algunas recomendaciones básicas.

• Asegurarse de que MATLAB no posee ya programada esa tarea. • Escoger un nombre para la función que sea nemotécnico.

76


Decomposition at level 3 : s = a3 + d3 + d2 + d1 . 0.4 0.2

s

0 −0.2 −0.4 0.4 0.2

a

3

0 −0.2 −0.4 0.02

d

3

0

−0.02

0.01

d

2

0

−0.01

0.01

d

1

0

−0.01 200

400

60 0

80 0

1000

1 200

1400

16 0 0

1 800

2000

Figura A.6: Análisis wavelet de la señal sig.mat con la haar y 3 niveles (Full decomposition display mode).

• Estudiar detenidamente los argumentos de la función (“input”) y el resultado (“output”) de la misma.

• Escribir las correspondientes explicaciones de ayuda. • Utilizar el máximo número de funciones de MATLAB que el programa original tiene.

Con estas recomendaciones y la ayuda que el programa ofrece para el comando function hacer los siguientes ejercicios.

Ejercicio A.6 Programar una función llamada energy que acepte como argumento una señal y calcule, como resultado, la energía de la señal.

A.3 Programando algunas funciones en MATLAB

77

Ejercicio A.7 Programar una función llamada cumenergy que acepte como argumento una señal y calcule, como resultado, su perfil de energía acumulada. Esto es, el resultado será otra señal del mismo tamaño que represente el perfil de energía acumulada (normalizada) de la señal original.

Ejercicio A.8 Programar una función llamada rms que tenga como argumento dos señales y que calcule, como resultado, el error cuadrático medio (RMS) entre las dos señales. Como estas funciones serán de gran utilidad en el resto de los ejercicios incluimos aquí unos archivos con la programación de las funciones. En este tipo de archivo la primera línea será el nombre de la función que debe coincidir con el nombre del archivo y la sintaxis de la función. Archivo: Archivo: energy.m energy.m functi function on [y] = energy energy(x) (x) %ENE %ENERG RGY Y Comp Comput utes es the the ener energy gy of a sign signal al. . % y = ENER ENERGY GY(x (x) ) ret returns urns the the ener energy gy of x. % The The ener energy gy comp comput uted ed as the the sum sum of the the % squa square re of the the ele element ments s of x. %

See See also also CUME CUMENE NER RGY. GY.

y = norm(x)^2;

Archivo: Archivo: cumenergy.m cumenergy.m functi function on [y] = cumene cumenergy rgy(x) (x) %CUMEN %CUMENERG ERGY Y Comput Computes es the normal normalize ized d cumula cumulativ tive e energy energy profil profile. e. % y = CUME CUMENE NERG RGY( Y(x) x) crea create tes s a vect vector or y with with the the norm normal aliz ized ed % cumu cumula lati tive ve energ energy y prof profil ile e of vecto vector r x. %

See also ENERGY.

y = cumsum(x.^2) cumsum(x.^2)/(norm /(norm(x)^2 (x)^2); );

78


Archivo: Archivo: rms.m functi function on [e] = rms(x, rms(x,y) y) %RMS Computes the RMS error of two signals. % e = RMS( RMS(x, x,y) y) retu returns rns the the RMS RMS erro error r bet between ween x and and y. % Sign Signal als s x and and y must ust have have same same size size. . e = sqrt(norm(xsqrt(norm(x-y)^2/ y)^2/length length(x)); (x));

Apéndice B

Ejercicios Los ejercicios incluidos en este apéndice complementan los capítulos dedicados a la exposición teórica. Se han clasificado en secciones que de alguna manera siguen el desarrollo teórico del curso. Antes de comenzar los ejercicios resultará conveniente trabajar el apéndice A en el que se dan los conceptos básicos sobre el uso de MATLAB y su “Wavelet Toolbox”.

B.1 B.1

Trans ransfo form rmad ada a de Haar Haar de de una una seña señall

Dada una señal s, podemos calcular su transformada de Haar desde la línea de comandos. La señal transformada queda almacenada en el vector (C) si utilizamos la notación de los ejercicios anteriores. Esta señal puede ser dibujada con plot. Debe notarse que no es posible dibujar la transformada desde la herramienta herramienta gráfica. Veamos un ejemplo: supongamos que la variable s contiene la señal del fichero sig.mat que también puede ser creada como se explicó en la sección A.2.1. [C,L] = wavedec(s,1 wavedec(s,1,’haar ,’haar’); ’); plot(C); axis([ axis([0 0 length length(s) (s) -1 1]); 1]); (véase la figura B.1)

Si queremos hacer dos transformadas consecutivas y observar la acumulación de energía en la parte correspondiente a la tendencia haremos lo siguiente. [C,L] = wavedec(s,2 wavedec(s,2,’haar ,’haar’); ’);

79

80

Ejercicios

1 0.8 0.6 0.4 0.2 0 −0.2 −0.4 −0.6 −0.8 −1 0

500

1000

15 1500

2000

Figura B.1: Transformada de Haar de 1 nivel de la señal sig.mat.

subplot(2,2,1); plot(s); axis([ axis([0 0 length length(s) (s) -1 1]); 1]); title(’Señal’); subplot(2,2,2); plot(cumenergy(s)); axis([ axis([0 0 length length(s) (s) -0.5 -0.5 1.5]); 1.5]); title(’Perfi title(’Perfil l energía energía acum.’); acum.’); subplot(2,2,3); plot(C); axis([ axis([0 0 length length(s) (s) -1 1]); 1]); title(’Trans title(’Transforma formada da nivel 2’); subplot(2,2,4); plot(cumenergy(C)); axis([ axis([0 0 length length(s) (s) -0.5 -0.5 1.5]); 1.5]); title(’Perfi title(’Perfil l energía energía acum. transf.’); transf.’); (véase la figura B.2)

Ejercicio B.1 Con las señales generadas en el ejercicio A.4 calcular transformadas de Haar de diferentes niveles y construir gráficos similares al de la figura B.2.

B.2 B.2

Comp Compre resi sión ón de de seña señale less (pri (prime mero ross paso pasos) s)

Veamos un ejemplo de compresión de señales siguiendo el esquema dado en el capítulo 4. Supondremos que la variable s contiene la señal del fichero sig.mat

B.2 Compresión de señales (primeros pasos)

81

Señal

Perfil energía acum.

1

1.5

0.5

1

0

0.5

−0.5

0

−1 0

500

1000

1500

2000

−0.5 0

Transformada nivel 2 1.5

0.5

1

0

0.5

−0.5

0

500

1000

1500

1000

1500

2000

Perfil energía acum. transf.

1

−1 0

500

2000

−0.5 0

500

1000

1500

2000

Figura B.2: Señal sig.mat y su transformada de Haar de nivel 2. A la derecha tenemos sus perfiles de energía acumulada.

que puede bajarse de la página web del curso o crearse muestreando la función dada al comienzo de la sección A.2.1 (ya se usó anteriormente). Para la compresión usaremos el fichero compress.m que incluimos a continuación. Este fichero puede ser fácilmente modificado para posteriores ejercicios. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Ficher Fichero o compre compress. ss.m m %% Arch Archiv ivo o de coma comand ndos os para para comp compri rimi mir r una una seña señal l %% La variable externa s contiene la señal %%%%%%%% %%%%%%%% Parámetros Parámetros a fijar %%%%%%%%%%% %%%%%%%%%%%%%%%%% %%%%%%%%%%% %%%%% ls = length(s); % Longitud de la señal w = ’haar’; % Fijamos la wavelet n = 10; % Niveles de descomposición thr_met = ’h’; % Método threslholding en_con en_cons s = 0.9999 0.9999; ; % Porcen Porcentaj taje e energí energía a a conser conservar var [C,L] [C,L] = wavede wavedec(s c(s,n, ,n,w); w); % Descom Descompos posici ición ón % Reordenamos Reordenamos valores valores absolutos absolutos transformad transformada a (decreciente (decrecientes) s)

82

Ejercicios

C_dec = abs(sort(-a abs(sort(-abs(C)) bs(C))); ); % Índice Índice hasta hasta acumul acumular ar un porcen porcentaj taje e prefij prefijado ado de energí energía a % Los 1 de esta variable son los índices de los valores que sobran ind_sobran ind_sobran = find(cumener find(cumenergy(C_d gy(C_dec)>= ec)>=en_con en_cons); s); % Valor correspondiente al primer índice que sobra es el umbral umbral umbral = C_dec( C_dec(ind ind_so _sobra bran(1 n(1)); )); % Valor Valor umbral umbral (thres (threshol hold) d) % Anulamos valores de la transformada que no pasan el umbral C_sob = wthresh(C,t wthresh(C,thr_met hr_met,umbra ,umbral); l); s_rec s_rec = wavere waverec(C c(C_so _sob,L b,L,w) ,w); ; % Recomp Recompone onemos mos la señal señal % Gráfico de las dos señales figure; subplot(2,1, subplot(2,1,1); 1); plot(s); plot(s); axis([1 axis([1 ls min(s) min(s) max(s)]); max(s)]); title(’Orig title(’Original’) inal’); ; subplot(2,1, subplot(2,1,2); 2); plot(s_rec); plot(s_rec); axis([1 axis([1 ls min(s_rec) min(s_rec) max(s_rec)]) max(s_rec)]); ; title(’Recon title(’Reconstrucc strucción’) ión’); ; % Calcul Calculamo amos s alguno algunos s datos datos sobre sobre la compre compresió sión n % Error RMS entre la señal original y la reconstruida error error = rms(s, rms(s,s_r s_rec) ec); ; map = C_sob~=0; % Mapa de significancia val_si val_sig g = sum(ma sum(map); p); % Valore Valores s signif significa icante ntes s % Factor de compresión (real y en forma aproximada) [comp,err] [comp,err] = sprintf(’%d: sprintf(’%d:%d’,ls %d’,ls,val_ ,val_sig); sig); [comp_ap,err [comp_ap,err] ] = sprintf(’%d: sprintf(’%d:%d’,ro %d’,round(ls und(ls/val_s /val_sig),1) ig),1); ; %% Resume Resumen n de result resultado ados s [l_long_orig [l_long_orig,err] ,err] = sprintf(’Lon sprintf(’Longitud gitud original: original: %d \n’,ls); \n’,ls); [l_val_sig,e [l_val_sig,err] rr] = sprintf(’Val sprintf(’Valores ores significati significativos: vos: %d \n’,val_sig \n’,val_sig); ); [l_com [l_comp,e p,err] rr] = sprint sprintf(’ f(’Fac Factor tor compre compresió sión n de %s \n’,co \n’,comp) mp); ; [l_comp_ap,e [l_comp_ap,err] rr] = sprintf(’\t sprintf(’\t (aproximadam (aproximadamente ente de %s) \n’,comp_ap \n’,comp_ap); ); [l_rms [l_rms,er ,err] r] = sprint sprintf(’ f(’Err Error or RMS: RMS: %d \n’,er \n’,error ror); ); sprintf(’%s%s%s%s%s’,... l_long_orig,l_val_sig,l_comp,l_comp_ap,l_rms) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Cuando ejecutamos este fichero obtendremos los gráficos de la figura B.3 y además nos muestra el siguiente resumen de resultados. Longitud Longitud original: original: 2048 Valores Valores significati significativos: vos: 573

B.2 Compresión de señales (primeros pasos)

83

Factor Factor compre compresió sión n de 2048:5 2048:573 73 (aproximada (aproximadamente mente de 4:1) Error RMS: 1.764582e-0 1.764582e-003 03 original 0.4 0.2 0 −0.2 −0.4

200

400

600

800

1000

1200

1400

1600

1800

2000

1400

1600

1800

2000

reconstrucción 0.4 0.2 0 −0.2 −0.4

200

400

600

800

1000

1200

Figura B.3: En el gráfico de la izquierda tenemos la señal original sig.mat. A la derecha se representan la reconstrucción de la señal después de comprimirla.

Ejercicio B.2 Comprimir cada una de las funciones del ejercicio A.4 conservando el 99.99 % de la energía. Utilizar Utilizar la wavelet wavelet de Haar y el nivel máximo de transformadas que sea posible. Para cada señal comprimida hacer el siguiente análisis. 1. Dibujar Dibujar la señal original y la señal reconstru reconstruida ida a partir de la compresión. compresión. 2. Calcular Calcular el error cometido cometido (RMS). 3. Contar Contar el número úmero de elemen elementos tos de la transf transform ormada ada que son nulos, nulos, esto esto es, el complementario de los que han sobrevivido al corte. 4. Calcular Calcular el factor de compresión compresión de la señal (tomaremo (tomaremoss como factor de compresión número de datos que tiene la señal original y el número de datos no nulos de la transformada que sobrevivieron al umbral, expresados en formato N:M; además se suele normalizar M a 1).

84

Ejercicios

Ejercicio B.3 Repetir el ejercicio anterior utilizando la herramienta gráfica. Con el método de “Global threshold” los resultados deben ser similares. Intentar mejorar la compresión con el método de “By level thresholding”. Mejorar en el sentido de misma energía retenida pero mayor número de ceros en la transformada después de eliminar aquellos valores que no pasan el umbral o umbrales fijados.

Ejercicio B.4 Construir señales de 1024 datos en el intervalo [0 intervalo [0,, 10] 10] para para cada una de las siguientes funciones. (a) f ( f (x) = x (b) f ( f (x) = 2[2 < 2[2 < x < 4]

− 2[5 2[5 < < x < 7] + 2[8 < 2[8 < x < 9] (c) f ( f (x) = 2[2 < 2[2 < x < 4] − x[5 < [5 < x < 7] + 2[8 < 2[8 < x < 9] (d) f ( f (x) = x(10 x (10 − x) 1. Calcular Calcular el factor factor de compresión compresión para cada señal utilizand utilizandoo la wavelet wavelet de Haar y reteniendo reteniendo el 99.99 % de la energía original. original. Calcular Calcular el error entre entre la señal original y la reconstruida a partir de la transformada después de eliminar los valores que no pasaron el umbral. 2. Calcul Calcular ar el factor factor de compre compresió sión n para para cada cada señal señal utiliza utilizand ndoo la wave wavelet let (db2) y reteniendo el 99.99 % de la energía original. Calcular el error entre entre la señal original y la reconstruida a partir de la transformada después de eliminar los valores que no pasaron el umbral. 3. Calcul Calcular ar el factor factor de compre compresió sión n para para cada cada señal señal utiliza utilizand ndoo la wave wavelet let (coif1) y reteni reteniend endoo el 99.99 99.99 % de la energí energíaa origin original. al. Calcular Calcular el error error entre la señal original y la reconstruida a partir de la transformada después de eliminar los valores que no pasaron el umbral. 4. Basado Basado en los cálculos cálculos anteriores anteriores argumentar argumentar qué wavel wavelet et es más apropiada para comprimir cada señal.

Ejercicio B.5 Para las funciones (b) y (c) del ejercicio B.4 tomar los 50 coeficientes mayores en valor absoluto de la transformada de Haar y calcular su transformada inversa. ¿Qué señal reconstruida es la que mejor se aproxima a la original y por qué? Calcular el error cometido en el intervalo [2. [2.5, 3.5] y en [5. [5.5, 6.5]. 5]. ¿Por qué son los dos errores similares en el intervalo [2. [2.5, 3.5] y 5] y no en [5. [5.5, 6.5]? 5]?

B.3 Otras wavelets ortogonales: Daubechies y Coiflets

B.3 B.3

85

Otra Otrass wavelets elets ortog ortogon onal ales: es: Daubec Daubechi hies es y CoiCoiflets

En esta sección pretendemos “visualizar” la mejoría que supone el uso de la transformada wavelet con Coiflets respecto de las Daubechies. Esta mejoría hay que entenderla en el sentido de que la tendencia “reproduzca” mejor la señal original. Antes de ello necesitamos explicar algunas consideraciones propias del programa MATLAB. Tomemos una señal y calculemos transformadas de nivel 1 con las wavelets db2 y coif1. A continuación inspeccionamos las longitudes de las subseñales tendencia. tendencia. Observ Observaremos aremos que estas longitude longitudess no son la mitad de la señal original. Esto ocurre porque MATLAB añade algunos coeficientes “extra” para asegurar la perfecta reconstrucción de la señal original. Estos problemas tienen relación con los valores de la señal en los bordes. Hay varios métodos para minimizar esta perturbación aunque nosotros deseamos precisamente que el programa no añada ningún coeficiente. En MATLAB existe una variable global que le dice al programa qué método debe utilizar para minimizar la distorsión producida en los bordes de la señal. Se trata del “Discr “Discrete ete Wavele Wavelet t Transf Transform orm Extens Extension ion Mode” Mode”. La orden dwtmode, sin argumentos nos dirá qué modo de extensión está activo. Para que MATLAB calcule tendencias y fluctuaciones de exactamente la mitad de la señal original debemos asegurarnos en primer lugar que estamos trabajando con señales de longitud par (en nuestro caso una potencia de 2). También debe usarse el DWT Extens p or tanto Extension ion Mode Mode llamado per. Debemos por asegurarnos de que está activo o activarlo con el siguiente comando. dwtmode(’per’)

Ejercicio B.6 Consideremos la función g (x) = 20x 20x2 (1

cos(12πx)) − x)4 cos(12πx

definida en el intervalo [0, [0, 1]. 1]. 1. Construir Construir una una señal de 214 datos muestreando g (x) en el intervalo [0, [0, 1]. 1]. 2. Extraer Extraer la subseñal subseñal tendencia tendencia de nivel nivel 2 para la wavelet wavelet db2. 3. Construir Construir la señal que represen representa ta la función 2 función 2gg(4x (4x) en el intervalo [0 intervalo [0,, 1]/ 1]/4 con 214 /22 datos.

86

Ejercicios

4. Comparar Comparar la subseñal subseñal tendencia tendencia con esta nuev nuevaa señal, señal, tanto tanto visualmente visualmente como calculando el error entre ellas (error máximo y RMS). 5. Repetir los últimos últimos tres pasos con la wavel wavelet et coif1. 6. Comparar Comparar los resultados resultados obtenidos. obtenidos. Para la realización de los ejercicios B.6, B.7 y B.8 se puede usar el archivo de comandos daub_coif.m. Este archivo se muestra a continuación para el ejercicio B.6. Cuando es ejecutado aparece la figura B.4 y este resumen de resultados. Error Error Error Error Error Error Error Error

máximo máximo para para máximo máximo para para rms para para la rms para para la

la db2 2.1720 2.172051e 51e-00 -003 3 la coif1 coif1 2.3026 2.302663e 63e-00 -008 8 db2 8.9669 8.966900e 00e-00 -004 4 coif1 coif1 1.9181 1.918157e 57e-00 -009 9

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Fichero: Fichero: daub_coif.m daub_coif.m %% Órdenes para ilustrar la mejoría de las coiflets %% resp respec ecto tos s de las las daub daubec echi hies es en algu alguno nos s caso casos s %%%% %% El DWT Extension Mode debe ser ’per’ %%%% clear; %% Señal( Señal(t) t) t1 = 0:(1-0)/(2^ 0:(1-0)/(2^14-1): 14-1):1; 1; s1 = 20*(t1.^2). 20*(t1.^2).*((1-t *((1-t1).^4) 1).^4).*cos( .*cos(12*pi* 12*pi*t1); t1); n = 2; % Niveles w1 = ’db2’; % Wavelet [C1,L1 [C1,L1] ] = wavede wavedec(s c(s1,n 1,n,w1 ,w1); ); % Descom Descompos posici ición ón de la señal señal tend_w tend_w1 1 = appcoe appcoef(C f(C1,L 1,L1,w 1,w1,n 1,n); ); % Tenden Tendencia cia % (sqrt(2)^n) (sqrt(2)^n)*Señal *Señal((2^n) ((2^n)*t) *t) en 1/(2^n)*dom 1/(2^n)*dominio inio t2 = t1(1:(2^n): t1(1:(2^n):length length(t1))/ (t1))/(2^n); (2^n); s2 = (sqrt(2)^n) (sqrt(2)^n)*(20*( *(20*((2^n*t (2^n*t2).^2) 2).^2)... ... .*((1-(2^n*t2)).^4).*cos(12*pi*(2^n*t2))); figure;

B.3 Otras wavelets ortogonales: Daubechies y Coiflets

87

subplot(3,1, subplot(3,1,1); 1); plot(s2); plot(s2); axis([1 axis([1 length(s2) length(s2) min(s2) min(s2) max(s2)]); max(s2)]); [tit,e [tit,err] rr] = sprint sprintf(’ f(’Grá Gráfic fica a de sqrt(2 sqrt(2)^% )^%d d g(2^%d g(2^%d x)’,n, x)’,n,n); n); title(tit); subplot(3,1, subplot(3,1,2); 2); plot(tend_w1 plot(tend_w1); ); axis([1 axis([1 length(tend length(tend_w1) _w1) min(tend_w1) min(tend_w1) max(tend_w1) max(tend_w1)]); ]); [tit_w [tit_w1,e 1,err] rr] = sprint sprintf(’ f(’Ten Tenden dencia cia de g(x) g(x) a nivel nivel %d con %s’,n, %s’,n,w1) w1); ; title(tit_w1); max_error_w1 = max(abs(tend_w1-s2)); rms_error_w1 rms_error_w1 = rms(tend_w1, rms(tend_w1,s2); s2); w2 = ’coif1’; % Wavelet [C2,L2 [C2,L2] ] = wavede wavedec(s c(s1,n 1,n,w2 ,w2); ); % Descom Descompos posici ición ón de la señal señal tend_w tend_w2 2 = appcoe appcoef(C f(C2,L 2,L2,w 2,w2,n 2,n); ); % Tenden Tendencia cia subplot(3,1, subplot(3,1,3); 3); plot(tend_w2 plot(tend_w2); ); axis([1 axis([1 length(tend length(tend_w2) _w2) min(tend_w2) min(tend_w2) max(tend_w2) max(tend_w2)]); ]); [tit_w [tit_w2,e 2,err] rr] = sprint sprintf(’ f(’Ten Tenden dencia cia de g(x) g(x) a nivel nivel %d con %s’,n, %s’,n,w2) w2); ; title(tit_w2); max_error_w2 = max(abs(tend_w2-s2)); rms_error_w2 rms_error_w2 = rms(tend_w2, rms(tend_w2,s2); s2); %% Resultados Resultados [l_max_err_w [l_max_err_w1,err 1,err] ] = sprintf(... sprintf(... ’Error ’Error máximo máximo para para la %s %d \n’,w1 \n’,w1,ma ,max_e x_erro rror_w r_w1); 1); [l_max_err_w [l_max_err_w2,err 2,err] ] = sprintf(... sprintf(... ’Error ’Error máximo máximo para para la %s %d \n’,w2 \n’,w2,ma ,max_e x_erro rror_w r_w2); 2); [l_rms_err_w [l_rms_err_w1,err 1,err] ] = sprintf(... sprintf(... ’Error ’Error rms para para la %s %d \n’,w1 \n’,w1,rm ,rms_e s_erro rror_w r_w1); 1); [l_rms_err_w [l_rms_err_w2,err 2,err] ] = sprintf(... sprintf(... ’Error ’Error rms para para la %s %d \n’,w2 \n’,w2,rm ,rms_e s_erro rror_w r_w2); 2); sprintf(’%s%s%s%s’,... l_max_err_w1,l_max_err_w2,l_rms_err_w1,l_rms_err_w2) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Ejercicio B.7 Consideremos la función g(x) = 20x 20x2 (1

cos(64πx)) + 30x 30x2 (1 − x)4 sen(30πx sen(30πx)) − x)2 cos(64πx

definida en el intervalo [0, [0, 1]. 1].

88

Ejercicios

2

2

Gráfica de sqrt(2) g(2 x) 0.5 0 −0.5 500

1000 1500 2000 2500 3000 Tendencia de g(x) a nivel 2 con db2

3500

4000

500

1000 1500 2000 2500 3000 Tendencia de g(x) a nivel 2 con coif1

3500

4000

500

1000

3500

4000

0.5 0 −0.5

0.5 0 −0.5 1500

2000

2500

3000

Figura B.4: De arriba a abajo tenemos la función 2g (4x (4x), la tendencia a nivel 2 con la db2 y finalmente la tendencia al mismo nivel con la coif1

1. Constr Construir uir una señal señal de 214 datos muestreando g (x) en el intervalo [0, [0, 1]. 1]. 2. Extraer Extraer la subseñal subseñal tendencia tendencia de nivel nivel 3 para la wavelet wavelet db2.

√

3. Constr Construir uir la señal señal que repres represen enta ta la funció función n 2 2g (8x (8x) en el intervalo [0, [0, 1]/ 1]/8 con 214 /23 datos. 4. Comparar Comparar la subseñal subseñal tendencia tendencia con esta nueva nueva señal, señal, tanto tanto visualmente visualmente como calculando el error entre ellas. 5. Repetir los últimos últimos tres pasos con la wavele wavelett coif1. 6. Comparar Comparar los resultados resultados obtenidos. obtenidos. db3, db10 db10. ¿Puedes obEjercicio B.8 Repetir el ejercicio con las wavelets db3,

servar el mismo fenómeno?

B.4

Dibuja Dibujand ndo o wa wavelets elets and scalin scaling g

Nuestro propósito es dibujar de una manera sencilla las wavelets y las scaling de los distintos niveles. Para ello se proponen una serie de ejercicios de experi-

B.4 Dibujando wavelets and scaling

89

mentación, con dimensión “pequeña”. La idea es obtener algunas conclusiones por exploración exploración directa sobre los resultado resultadoss y, finalment finalmente, e, generalizar generalizar estos resultados a cualquier dimensión. Antes de comenzar necesitamos fijar le “DWT Extension Mode” a ’per’. dwtmode(’per’)

Ejercicio B.9 Supongamos que estamos en el espacio RN con N = 23 , es decir, todas nuestras señales tendrán 8 coordenadas. Utilizar el fichero llamado wave_scal.m como base para realizar las siguientes acciones. 1. Fijar Fijar la wavele avelett a utilizar utilizar como la de ’haar’. 2. Fijar el número número de niveles de descomposición descomposición a 1. 3. Construir Construir una señal señal de zeros de longitud longitud N . N . 4. Calcular Calcular su transformada transformada wavel wavelet et de Haar de nivel 1. El vector vector C contendrá los coeficientes de la transformada y L es el vector de longitudes (en este caso con valores 22 , 22 , 23 ). 5. Cambia Cambiarr el vector ector C de la transformada de la señal al vector de la base canónica e canónica e1 , es decir, contiene un 1 en la posición 1 y ceros en las restantes 7 posiciones. 6. Calcular Calcular la transform transformada ada wavel wavelet et inversa inversa de Haar de nivel 1 de la señal C. Dibujar la señal reconstruida con la orden stem y anotar los resultados observados. 7. Repetir los dos últimos pasos con los restantes restantes vectores vectores de la base canónica ei , con i = 2, . . . , 8. Cuando se ejecuta el fichero wave_sca.m se obtiene la figura B.5. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Ficher Fichero o wav_sc wav_sca.m a.m %% Para Para constr construir uir scalin scaling g y wavele wavelets ts %%%%%%%% %% El DWT Extension Mode debe ser ’per’ %%%%%%%% clear;

90

Ejercicios −1

−1

H1 (e1)

H1 (e5)

1

1

0.5

0

0 0

2

4 −1 H1 (e2)

6

8

−1

1

1

0.5

0

0 0

2

4 −1 H1 (e3)

6

8

−1

1

1

0.5

0

0 0

2

4

−1 H1 (e4)

6

8

1

2

4 −1 H1 (e6)

6

8

0

2

4 −1 H1 (e6)

6

8

0

2

4 −1 H1 (e8)

6

8

0

2

4

6

8

1

0.5 0 0

−1

0

0 2

4

6

8

−1

Figura B.5: Valores de las señales reconstruidas utilizando la transformada inversa de Haar con N = 23 .

k = 3; N = 2^k; m = 1; % Niveles w = ’haar’; % Wavelet % Construimos L a partir de una señal de ceros s = zeros( zeros(1,N 1,N); ); [C,L] = wavedec(s,m wavedec(s,m,w); ,w); C = zeros( zeros(siz size(C e(C)); )); % Anulam Anulamos os los valore valores s i = 1; C(i) = 1; s_rec = waverec(C,L,w); figure; subplot(4,2, subplot(4,2,1); 1); stem(s_rec); stem(s_rec); title(’H_1^-1(e_1)’); C=zero C=zeros(s s(size ize(C) (C)); ); % Anulam Anulamos os los valore valores s i = 2; C(i) = 1; s_rec = waverec(C,L,w); subplot(4,2, subplot(4,2,3); 3); stem(s_rec); stem(s_rec);

B.4 Dibujando wavelets and scaling

91

title(’H_1^-1(e_2)’); C = zeros( zeros(siz size(C e(C)); )); % Anulam Anulamos os los valore valores s i = 3; C(i) = 1; s_rec = waverec(C,L,w); subplot(4,2, subplot(4,2,5); 5); stem(s_rec); stem(s_rec); title(’H_1^-1(e_3)’); C = zeros( zeros(siz size(C e(C)); )); % Anulam Anulamos os los valore valores s i = 4; C(i) = 1; s_rec = waverec(C,L,w); subplot(4,2,7); stem(s_rec); title(’H_1^-1(e_4)’); C = zeros( zeros(siz size(C e(C)); )); % Anulam Anulamos os los valore valores s i = 5; C(i) = 1; s_rec = waverec(C,L,w); subplot(4,2, subplot(4,2,2); 2); stem(s_rec); stem(s_rec); title(’H_1^-1(e_5)’); C = zeros( zeros(siz size(C e(C)); )); % Anulam Anulamos os los valore valores s i = 6; C(i) = 1; s_rec = waverec(C,L,w); subplot(4,2, subplot(4,2,4); 4); stem(s_rec); stem(s_rec); title(’H_1^-1(e_6)’); C = zeros( zeros(siz size(C e(C)); )); % Anulam Anulamos os los valore valores s i = 7; C(i) = 1; s_rec = waverec(C,L,w); subplot(4,2, subplot(4,2,6); 6); stem(s_rec); stem(s_rec); title(’H_1^-1(e_6)’); C = zeros( zeros(siz size(C e(C)); )); % Anulam Anulamos os los valore valores s i = 8; C(i) = 1; s_rec = waverec(C,L,w); subplot(4,2, subplot(4,2,8); 8); stem(s_rec); stem(s_rec); title(’H_1^-1(e_8)’); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Ejercicio B.10 Repetir el ejercicio ejercicio anterior anterior con 2 niveles de transformada, transformada, es decir, el vector de longitudes L tiene valores 2 valores 2,, 2, 22 , 23 . Observar los resultados

92

Ejercicios

y buscar diferencias o similitudes con el ejercicio anterior.

Ejercicio B.11 Repetir el ejercicio anterior con 3 niveles de transformada. Observar los resultados y buscar diferencias o similitudes con el ejercicio anterior.

Ejercicio B.12 Generali Generalizar zar los result resultado adoss anter anterior iores. es. Esto Esto es, suponga supongamos mos que tenemos el espacio RN , con N con N = 2n (tendremos n (tendremos n niveles niveles de transformada wavelet como máximo). Si tomamos un vector cualquiera de la base canónica de RN , por ejemplo ei (con i = 1, . . . , N ), ) , y calculamos su transformada wavelet inversa de nivel m (m n) ¿Qué señal obtendremos? Para contestar el ejercicio completar el siguiente esquema

≤

ei , i = 1, . . . , N

−1 DW T m

−→

  

. . .. . .. . . . . .. . .. . . . . .. . .. . . .. .

i = . = . . . . . . . . . i = . = . . . . . . . . . i = . = . . . . . . . . . .. .

. . . . . . . . . i = . = . . . . . . . . .

Ejercicio B.13 Tomar ahora N = 210 = 1024. 1024. 1 , w1 , v 2 , (a) Para Para la wavel wavelet et db2 calcular las señales scaling y wavelets v25 20 4 3 5 6 w4 , w2 , v4 .

(b) Repetir el apartado anterio anteriorr con la wavelet wavelet coif1.

B.5 B.5

Comp Compre resi sión ón de de seña señale less y cuan cuanti tiza zaci ción ón

En esta sección ampliaremos las herramientas de compresión de señales incluyendo el tema de la cuantización uniforme y codificación. El proceso de cuantizar los datos puede ser programado sin dificultad en MATLAB y así lo haremos en el ejercicio B.14. Para el proceso de codificación (al codificar también cuantizaremos) vamos a utilizar funciones que no pertenecen a la “Wavelet Toolbox” sino a la “Communications Toolbox” que también deberá estar disponible en nuestro entorno MATLAB. Nos aseguraremos que el “DWT Extension Mode” se ha vuelto a fijar a su valor por defecto, es decir, a ’sym’. dwtmode(’sym’)

B.5 Compresión de señales y cuantización

93

Ejercicio B.14 Vamos a comprimir comprimir una señal tal y como se hizo en la sección B.2 pero experimentaremos el proceso de cuantización con distinto número de intervalos, comparando el resultado. Supondremos que la variable s contiene la señal del fichero sig.mat que puede bajarse de la página web del curso o crearse muestreando la función dada al comienzo de la sección A.2.1. Para la compresión y cuantización usaremos el fichero compress2.m que incluimos a continuación. Este fichero puede ser fácilmente modificado para posteriores pruebas cambiando el número de interv intervalos alos en la cuantizaci cuantización. ón. Comprimiremos Comprimiremos la señal conservando conservando el 99.99 % de la energía y usaremos para la cuantización 32, 64, 128 y 256 intervalos. Los resultados pueden verse en la figura B.6. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Fichero Fichero compress2.m compress2.m %% Compresión de una señal %% Incluye el tema de la cuantización %% pero pero no codifi codificac cación ión % Arch Archiv ivo o de coma comand ndos os para para comp compri rimi mir r una una seña señal l % La variable s contiene la señal %%%%%%%% %%%%%%%% Parametros Parametros a fijar %%%%%%%%%%% %%%%%%%%%%%%%%%%% %%%%%%%%% %%% ls = length(s); % Longitud de la señal med = mean(mean(s)); % Media de la señal w = ’haar’; % Fijamos la wavelet n = 10; % niveles de descomposición thr_met = ’h’; % método threslholding Econ Econs s = 0.99 0.9999 99; ; % Porc Porcen enta taje je ener energí gía a a cons conser erva var r q = 32; 32; % Númer Número o inter interva valo los s cuant cuantiz izac ació ión n (par) (par) s = s - med; % Centramos valores [C,L] [C,L] = wavede wavedec(s c(s,n, ,n,w); w); % Descom Descompos posici ición ón % Reordenamos Reordenamos valores valores absolutos absolutos transformad transformada a (decreciente (decrecientes) s) C_dec = abs(sort(-a abs(sort(-abs(C)) bs(C))); ); % Busc Buscam amos os el índi índice ce hast hasta a acum acumul ular ar un porc porcen enta taje je alto alto de ener energí gía a % Los 1 de esta variable son los índices de los valores que sobran ind_sobran ind_sobran = find(cumener find(cumenergy(C_d gy(C_dec)>= ec)>=Econs) Econs); ; % El valor correspondiente al primer índice que sobra es el umbral umbral umbral = C_dec( C_dec(ind ind_so _sobra bran(1 n(1)); )); % Valor Valor umbral umbral (thres (threshol hold). d).

94

Ejercicios

% Anulamos valores de la transformada que no pasan el umbral C_sob = wthresh(C,t wthresh(C,thr_met hr_met,umbra ,umbral); l); %% Cuantizació Cuantización n alpha alpha = max(ab max(abs(C s(C_so _sob)) b)); ; delta delta = 2*alph 2*alpha/q a/q; ; %% Cent Centro ros s inte interv rval alos os en vect vector or z a = -alpha -alpha+de +delta lta/2: /2:del delta: ta:-3* -3*del delta/ ta/2; 2; b = 3*delt 3*delta/2 a/2:de :delta lta:al :alpha pha-de -delta lta/2; /2; z = [a 0 b]; %% La mitad de la longitud de cada intervalo en w lon = [(delta/2)* [(delta/2)*ones(1 ones(1,lengt ,length(a)) h(a)) delta (delta/2)*ones(1,length(b))]; %% Cuanti Cuantizam zamos os la transf transform ormada ada C_cuan C_cuant t = C_sob; C_sob; for i=1:le i=1:lengt ngth(C h(C_cu _cuant ant), ), for j=1:length( j=1:length(z), z), if abs(C_cuant( abs(C_cuant(i)-z(j i)-z(j))<=l ))<=lon(j) on(j) C_cuant(i)=z C_cuant(i)=z(j); (j); end; end; end; s_rec s_rec = wavere waverec(C c(C_cu _cuant ant,L, ,L,w); w); % Recomp Recompone onemos mos la señal señal s_rec = s_rec + med; % Añadimos media % Gráfico de las dos señales figure figure; ; subplo subplot(2 t(2,1, ,1,1); 1); plot(s plot(s); ); axis([ axis([1 1 ls min(s) min(s) max(s) max(s)]); ]); title(’origi title(’original’) nal’); ; subplot(2,1, subplot(2,1,2); 2); plot(s_rec) plot(s_rec); ; axis([1 axis([1 ls min(s_rec) max(s_rec)]); title(’reconstrucción’); % Calcul Calculamo amos s alguno algunos s datos datos sobre sobre la compre compresió sión n % Error RMS entre la señal original y la reconstruida error error = rms(s, rms(s,s_r s_rec) ec); ; map = C_sob~=0; % Mapa de significancia val_si val_sig g = sum(ma sum(map); p); % Valore Valores s signif significa icante ntes s % Factor Factor de compre compresió sión n y energí energía a conser conservad vada a [comp,err] [comp,err] = sprintf(’%d: sprintf(’%d:%d’,ls %d’,ls,val_ ,val_sig); sig); [comp_aprox, [comp_aprox,err] err] = sprintf(’%d: sprintf(’%d:%d’,r %d’,round(l ound(ls/val_ s/val_sig),1 sig),1); ); %% Resume Resumen n de result resultado ados s [l_long_orig [l_long_orig,err] ,err] = sprintf(’Lon sprintf(’Longitud gitud original: original: %d \n’,ls); \n’,ls); [l_val_sig,e [l_val_sig,err] rr] = sprintf(’Val sprintf(’Valores ores significati significativos: vos: %d \n’,val_sig \n’,val_sig); ); [l_com [l_comp,e p,err] rr] = sprint sprintf(’ f(’Fac Factor tor compre compresió sión n de %s \n’,co \n’,comp) mp); ; [l_com [l_comp_a p_apro prox,e x,err] rr] = sprint sprintf(’ f(’\t \t (aprox (aprox. . de %s) \n’,co \n’,comp_ mp_apr aprox) ox); ;


95

[l_rms [l_rms,er ,err] r] = sprint sprintf(’ f(’Err Error or RMS: RMS: %d \n’,er \n’,error ror); ); sprintf(’%s%s%s%s%s’,... l_long_orig,l_val_sig,l_comp,l_comp_aprox,l_rms) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

original 0.4 0.2 0 −0.2 −0.4 0.4 0.2 0 −0.2 −0.4 0.4 0.2 0 −0.2 −0.4 0.4 0.2 0 −0.2 −0.4

600

1400 160 1600 1800 1800 2000 2000 800 1000 1000 1200 1200 1400 reconstrucción reconstrucción con q=32

20 0

400 6 00 00

1200 0 1400 1400 1600 1600 1800 800 2000 800 1000 120 reconstrucción reconstrucción con q=64

20 0

400 6 00 00

1200 0 1400 1400 1600 1600 1800 800 2000 800 1000 120 reconstrucción reconstrucción con q=128

4 00 200 40

200 40 4 00

600

20 0

600

1200 0 1400 1400 1600 1600 1800 800 200 2000 800 1000 120 reconstrucción reconstrucción con q=256

0.4 0.2 0 −0.2 400

80 0

1000 120 1200 0 1400 1400 1600 1600 1800 800 2000

Figura B.6: De arriba abajo tenemos la señal original y la reconstruida a partir de la compre compresió sión, n, reteni reteniend endoo el 99.99 99.99 % de la energí energíaa y usand usandoo cuan cuantizació tización n con q intervalos intervalos para q = 32 32,, 64 64,, 128 128,, 256 256..

Los ejercicios que vamos a realizar ahora se basarán en señales de audio (ficheros *.wav) que pueden obtenerse de la página web del curso (otros ficheros *.wav de los que disponga el lector también servirían). El método que vamos a explicar también se puede aplicar a otras señales que no sean de audio.

96

Ejercicios

Ejercicio B.15 Pretendemos comprimir la señal de audio greasy.wav. para ello seguiremos los siguientes pasos. El fichero comp_audio.m contiene las órdenes y comentarios necesarios. 1. Cargar Cargar la señal de audio en MA MATLAB TLAB.. Puede Puede reproducirse reproducirse si se desea. 2. Comprimir Comprimir la señal. Usar la wavel wavelet et coif5 y 14 niveles de transformada. Retene Retenerr el 99.99 % de la energí energíaa origin original al de la señal. señal. 3. Cuantizar Cuantizar la señal al mismo número número de bits que la original. Usar cuanticuantización uniforme. 4. Manipular Manipular los coeficientes coeficientes de la transformada transformada,, después después de anular anular los valores que no pasaron el umbral fijado, para obtener los bits necesarios para transmitir el mapa de significancia y el número de coeficientes no nulos a transmitir. 5. Calcul Calcular ar los bits por dato dato (bpp) (bpp) que son necesari necesarios os para para transm transmitir itir la señal. 6. Decodificar Decodificar la señal transforma transformada da y cuantizad cuantizada. a. 7. Recomponer Recomponer la señal original original con la transform transformada ada wavelet wavelet inversa inversa sobre la decodificación de los coeficientes transmitidos. Calcular el error RMS entre la señal reconstruida y la original. Reproducir si se desea la señal. %%%%%%%%%%%%%%%%%%%%% %% Fichero Fichero comp_audio. comp_audio.m m %% Compresión de una señal de audio %% Incluye el tema de la cuantización %% y tamb tambié ién n la codi codifi fica caci ción ón %% Utiliz Utiliza a para para ello ello alguno algunos s comand comandos os %% de la tool toolbo box x de comu comuni nica caci cion ones es %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % La variable file debe ser el nombre del fichero *.wav % % % %

Signif Significa icado do de las variab variables les s = señal de audio a comprimir fs = frecuencia en herzios nb = bits por dato


97

clear; clear; [s,fs, [s,fs,nb] nb] = wavrea wavread(f d(file ile); ); ls = length(s); % Longitud señal w = ’coif5’; % Wavelet a utilizar n = 14; % Niveles de transformada en_r en_ret et = 0.99 0.9999 99; ; % Ener Energí gía a rete reteni nida da thr_me thr_met t = ’h’; ’h’; % Método Método thresh threshold olding ing [C,L] [C,L] = wavede wavedec(s c(s,n, ,n,w); w); % Transf Transform ormada ada wavele wavelet t % Proces Proceso o de compre compresió sión n dec = abs(so abs(sort( rt(-ab -abs(C s(C))) ))); ; % Reorde Reordenam namos os left left = find(c find(cume umener nergy( gy(dec dec)>= )>=en_ en_ret ret); ); % Índice Índices s que sobran sobran thr = dec(le dec(left( ft(1)) 1)); ; % Valor Valor del umbral umbral (thres (threshol hold). d). Cthr Cthr = wthres wthresh(C h(C,th ,thr_m r_met, et,thr thr); ); % Anulam Anulamos os valore valores s de la transf transform ormada ada % Proces Proceso o de cuanti cuantizac zación ión (codif (codifica icació ción) n) peak peak = max(ab max(abs(C s(Cthr thr)); )); % Máximo Máximo (peak (peak value) value) Cquant = uencode(Cth uencode(Cthr,nb,p r,nb,peak); eak); % Cuantizamos Cuantizamos valores transformad transformada a % Result Resultado ados s de la compre compresió sión n map = Cthr~=0; % Mapa de significancia last_s last_sig_ ig_val val = max(fi max(find( nd(map map~=0 ~=0)); )); % Último Último índice índice signif significa icativ tivo o sig_val = sum(map); % Número de valores a transmitir de la transf. bpp = (last_ (last_sig sig_va _val+s l+sig_ ig_val val*nb *nb)/l )/ls; s; % Bits Bits por dato dato [bpp_orig,er [bpp_orig,err]=sp r]=sprintf( rintf(’Origi ’Original: nal: %3.3f bpp \n’,nb); \n’,nb); [bpp_comp,er [bpp_comp,err]=sp r]=sprintf( rintf(’Compr ’Compresión esión: : %3.3f bpp \n’,bpp); \n’,bpp); Cdecod = udecode(Cqu udecode(Cquant,nb ant,nb,peak) ,peak); ; % Decodificaci Decodificación ón s_rec = waverec(Cde waverec(Cdecod,L, cod,L,w); w); % Reconstrucc Reconstrucción ión %% Compar Comparaci ación ón y result resultado ados s error error = rms(s, rms(s,s_r s_rec) ec); ; % Entre Entre origin original al y recons reconstru truida ida [error_men,e [error_men,err]=s rr]=sprintf printf(’Erro (’Error r cometido: cometido: %d \n’,error); \n’,error); wavpla wavplay(s y(s_re _rec,f c,fs); s); % Reprod Reproducc ucción ión de audio audio si se desea desea sprintf(’%s%s%s’,bpp_orig,bpp_comp,error_men) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Ejercicio B.16 Repetir el ejercicio B.15 pero utilizar un umbral de corte para

98

Ejercicios

retener retener el 96 % de la energía energía de la señal. señal. Comparar Comparar los resultados. resultados.

Ejercicio B.17 Repetir el ejercicio B.15 utilizando una wavelet db15. Comparar parar los resultados. resultados.

Ejercicio B.18 Repetir el ejercicio ejercicio B.15 con la señal de audio call_back.wav u otras de las que disponga el lector.

B.6 B.6

Redu Reducc cció ión n del del ruid ruido o en una una seña señall

Esta sección está dedicada al proceso de reducción del ruido en una señal unidimensional. En primer lugar haremos unos ejercicios para visualizar el comportamiento del ruido blanco gaussiano cuando le aplicamos una transformada wavelet.

Observación B.19 En este caso los ejercicios se pueden realizar igual con el “DWT Extension Mode” activo a un valor o a otro ya que no utilizaremos las longitudes de las subseñales tendencia y fluctuación. La cuestión de que la señal original y su transformada tengan distinta longitud no es problema para dibujarlas conjuntamente. Sí puede ser un problema a la hora de realizar operaciones entre ellas. Por ejemplo, para operar (sumar, restar, ...) con una señal original y su transformada, ambas deben tener la misma longitud y necesariamente “DWT Extension Mode” debería ser fijado a ‘per’ y la longitud de la señal original debe ser par.

B.6.1 B.6.1

Comport Comportami amien ento to del del ruido ruido blanco blanco a trav través de una una transtransformada wavelet

En primer lugar vamos a ver cómo se comporta el ruido blanco cuando se calcula su transformada wavelet.

Ejercicio B.20 Las órdenes básicas necesarias para este ejercicio pueden encontrarse en el fichero white_noise.m; con pequeñas modificaciones puede ser utilizado para resolver todos los apartados. 1. Construir Construir una señal de ruido blanco blanco de longitud 1024. 1024. 2. Calcul Calcular ar su transf transform ormada ada haar de un nivel. Dibujar conjuntamente la señal original y la señal transformada (véase la figura B.7).

B.6 Reducción del ruido en una señal

99

3. Calcular Calcular su transform transformada ada haar de tres niveles. Dibujar conjuntamente la señal original y la señal transformada. 4. Calcular Calcular su transf transforma ormada da db2 de un nivel. Dibujar conjuntamente conjuntamente la señal original y la señal transformada. 5. Calcular Calcular su transformada transformada db2 de tres niveles. Dibujar conjuntamente la señal original y la señal transformada. 6. Calcular Calcular su transformada transformada coif1 de un nivel. Dibujar conjuntamente la señal original y la señal transformada. 7. Calcular Calcular su transformada transformada coif1 de tres niveles. Dibujar conjuntamente la señal original y la señal transformada.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Fichero Fichero white_noise white_noise.m .m %% Genera Generació ción n ruido ruido blanco blanco gaussi gaussiano ano y su transf transform ormada ada clear; N = 1024; s = wgn(N,1,1); w = ’haar’; % Wavelet n = 1; % Niveles [C,L] = wavedec(s,n wavedec(s,n,w); ,w); figure; subplot(2,1,1); plot(s,’k.’,’markersize’,3); title(’Ruido title(’Ruido original’); original’); axis([ axis([1 1 length length(s) (s) -5 5]); 5]); subplot(2,1,2); plot(C,’k.’,’markersize’,3); title(’Trans title(’Transformad formada a Haar’); Haar’); axis([ axis([1 1 length length(C) (C) -5 5]); 5]); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

100

Ejercicios

Ruido original 5

0

−5

100

20 0

3 00

4 00

5 00

600

700

800

9 00

1000

700

800

900

100 0

Transformada Haar 5

0

−5

100

200

300

400

50 0

6 00

Figura B.7: En la parte superior tenemos una señal formada por ruido blanco gaussiano y debajo su transformada Haar de 1 nivel.

B.6. B.6.2 2

Redu Reducc cció ión n del del ruid ruido o en una una señ señal al

Veamos como reducir la presencia de ruido en una señal. Supondremos que el ruido es gaussiano de media cero (ruido blanco). Si se dispone de señales de estas características se puede proceder a su reducción directamente; si no se dispone de señales con ruido lo añadiremos nosotros utilizando funciones que pertenecen a la “Communications Toolbox” de MATLAB.

Ejercicio B.21 Para este ejercicio se puede utilizar el fichero denoise.m para la parte que concierne a la reducción de ruido en una señal contaminada. Para crear la señal y añadirle ruido utilizamos el fichero make_noisy_signal.m. 1. Construir Construir una una señal de 211 puntos a partir de la función f ( f (x) = 4x2 (1

− x)3(2 − x)2 cos(18x cos(18x(1 + x))

en el intervalo [0,2]. 2. Añadir Añadir ruido blanco blanco gaussiano gaussiano de 25 dBW de potenc p otencia ia a esta señal.


101

3. Calcular Calcular la transformada transformada de nivel 9 de la señal contamin contaminada. ada. Utilizar Utilizar la wavelet coif5. 4. Obtener Obtener la desviación desviación típica típica σ σ de de una porción de la fluctuación fluctuación de primer nivel.

√

5. Fijar Fijar el umbra umbrall a σ 2log N . N . 6. Eliminar Eliminar los valores valores de la transformada transformada que no pasan el umbral. umbral. Utilizar Utilizar “hard threshold”. 7. Reconstrui Reconstruirr la señal a partir de los valores valores de la transform transformada ada que sobrevivieron al umbral. 8. Comparar Comparar la señal reconstruida reconstruida y la señal original. original. Calcular Calcular para ello el error RMS, el porcentaje de energía retenida y finalmente dibujarlas. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Fichero Fichero make_noisy_ make_noisy_signal signal.m .m %% Fich Ficher ero o para para cons constr trui uir r una una seña señal l y añad añadir irle le ruid ruido o blan blanco co gaus gaussi sian ano o clear; a = 0; % Initial value b = 2; % End value N = 2^11; % Total number of values x = a:(b-a a:(b-a)/( )/(N-1 N-1):b ):b; ; s = (4*(x.^2)).*((1-x).^3).*((2-x).^2).*cos((18*x).*(1+x)); s = awgn(s,25); % Añadimos ruido blanco de 25 dB %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Ficher Fichero o denois denoise.m e.m %% Eliminación del ruido de una señal %% La vari variab able le exte extern rna a s cont contie iene ne la seña señal l cont contam amin inad ada a w = ’coif5’; % Wavelet a utilizar n = 9; % Niveles de transformada thr_ thr_me met t = ’h’; ’h’; % Méto Método do thre thresh shol oldi ding ng: : Hard Hard

102

Ejercicios

% Numero de veces la desviación típica en umbral coef = sqrt(2*log(l sqrt(2*log(length( ength(s))); s))); [C,L] [C,L] = wavede wavedec(s c(s,n, ,n,w); w); % Transf Transform ormada ada wavele wavelet t % Cálcul Cálculo o del umbral umbral (thres (threshol hold) d) cD1 = detcoe detcoef(C f(C,L, ,L,1); 1); % Extrae Extraemos mos la primer primera a fluctu fluctuaci ación ón des_ti des_tip p = std(cD std(cD1); 1); % Desvia Desviació ción n típica típica thr = coef*d coef*des_ es_tip tip; ; % Umbral Umbral (thres (threshol hold) d) Cthr Cthr = wthres wthresh(C h(C,th ,thr_m r_met, et,thr thr); ); % Umbral Umbraliza izamos mos la transf transform ormada ada s_den s_den = wavere waverec(C c(Cthr thr,L, ,L,w); w); % Recons Reconstru truimo imos s %% Resulta Resultados dos y represe representa ntació ción n gráfic gráfica a figure figure; ; hold hold on; axis([ axis([1 1 length length(s) (s) min(s) min(s) max(s) max(s)]); ]); plot(s plot(s’,’ ’,’y’) y’); ; plot(s_den’, plot(s_den’,’k’); ’k’); hold off; error error = rms(s, rms(s,s_d s_den) en); ; % Error Error [error_men,e [error_men,err] rr] = sprintf(’Err sprintf(’Error: or: %d \n’,error); \n’,error); ene_ret ene_ret = 100*energy( 100*energy(s_den) s_den)/energ /energy(s); y(s); % Porcentaje Porcentaje energía energía retenida retenida [ene_ret_men [ene_ret_men,err] ,err] = sprintf(’Ene sprintf(’Energía rgía retenida: retenida: %2.3f %% \n’,ene_ret \n’,ene_ret); ); sprintf(’%s%s’,error_men,ene_ret_men) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Cuando ejecutamos los códigos de los ficheros make_noisy_signal.m y denoise.m se obtiene el gráfico de la figura B.8. El error RMS entre la señal con ruido y la filtrada es 0.0567; el porcentaje de energía retenida es de 87.842 87.842 %.

Ejercicio B.22 Repetir el ejercicio B.21 utilizando “soft threshold”. Comparar los resultados.

Ejercicio B.23 Repetir el ejercicio B.21 utilizando “hard threshold” y la wavelet db10. Niveles de transformada a vuestra decisión. Comparar los resultados.

Ejercicio B.24 Repetir el ejercicio B.21 utilizando “soft threshold” y la wavelet db10. Niveles de transformada a vuestra decisión. Comparar los resultados. La eliminación del ruido blanco en una señal ha sido ampliamente investigada. Hay múltiples maneras de seleccionar el umbral, algunas basadas en


103

0.4 0.3 0.2 0.1 0 −0.1 −0.2 −0.3 −0.4 200

40 400

600

800

1000

1200

1400

1600

1800

2000

Figura B.8: Señal original contaminada (en gris), superpuesta tenemos la señal obtenida después de reducir el ruido (en negro).

algoritmos realmente complicados. Algunas de estas técnicas de seleccionar el umbral han sido implementadas en la “Wavelet Toolbox”. No entraremos en detalles, aunque sí usaremos las órdenes programadas para saber al menos que existen y comparar resultados. Los siguientes ejercicios van encaminados en este sentido.

Ejercicio B.25 Repetir el ejercicio B.21 utilizando el comando wden. Es útil ver la ayuda a la selección del umbral (comando thselect y experimentar con ella). El fichero denoise_auto.m puede ser utilizado como base. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Fichero Fichero denoise_aut denoise_auto.m o.m %% Reducción del ruido de una señal de forma automática %% La vari variab able le exte extern rna a s cont contie iene ne la seña señal l cont contam amin inad ada a ls = length(s); % Longitud w = ’coif5’; % Wavelet a utilizar n = 9; % Niveles de transformadas %% Eliminamos el ruido de una con la orden de la toolbox. %% Fijamo Fijamos s alguno algunos s datos datos antes antes tptr = ’rigrsure’; % Método de selección del umbral sorh = ’h’; % Soft or hard threshold scal = ’one’;

104

Ejercicios

[s_den,C,L] [s_den,C,L] = wden(s,tptr, wden(s,tptr,sorh,s sorh,scal,n cal,n,w); ,w); % Reducción Reducción ruido %% Resulta Resultados dos y represe representa ntació ción n gráfic gráfica a figure; hold hold on; on; axis([1 axis([1 length(s) length(s) min(s) max(s)]); plot(s’,’y’); plot(s_den’,’k’); hold hold off; off; error error = rms(s, rms(s,s_d s_den) en); ; % Error Error [error_men,e [error_men,err] rr] = sprintf(’Err sprintf(’Error: or: %d \n’,error); \n’,error); ene_ret ene_ret = 100*energy( 100*energy(s_den) s_den)/energ /energy(s); y(s); % Porcentaje Porcentaje energía energía retenida retenida [ene_ret_men [ene_ret_men,err] ,err] = sprintf(’Ene sprintf(’Energía rgía retenida: retenida: %2.3f %% \n’,ene_ret \n’,ene_ret); ); sprintf(’%s%s’,error_men,ene_ret_men) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Cuando ejecutamos el códigos del denoise_auto.m con la misma señal del ejercicio B.21 se obtiene el gráfico de la figura B.9. El error RMS entre la señal con ruido y la filtrada es 0.0753; el porcentaje de energía retenida es de 76.416 76.416 %. 0.4

0.3

0.2

0.1

0

−0.1

−0.2

−0.3

−0.4 200

400

600

800

1000

1200

1400

1600

1800

2000

Figura B.9: Señal original contaminada (en gris), superpuesta tenemos la señal obtenida después de reducir el ruido con la orden wden (en negro).

Ejercicio B.26 Repetir Repetir el ejerci ejercicio cio B.21 B.21 utiliz utilizand andoo la herram herramien ienta ta gráfica gráfica.. Probar con cada uno de los métodos de selección del umbral. Si se desea se puede ir guardando los resultados para posteriores comparaciones.

B.7 Transformada de Fourier Discreta

B.7

105

Transfo ransforma rmada da de Fourier ourier Discreta Discreta

B.7.1

“Fast “Fast Fourier ourier Transform” ransform”

Ejercicio B.27 Calcularemos la transformada de Fourier de una señal y a continuación reconstruiremos la señal transformada inversa de Fourier. El fichero fourier_1.m contiene las órdenes necesarias y el resultado se puede ver en la figura B.10. 1. Crear Crear una señal señal a partir partir de la función función h(t) = t sen4t sen4t en [0, [0, 2π] con 2048 puntos de muestreo. 2. Calcular Calcular su transform transformada ada de Fourier ourier y represen representarla tarla.. 3. Calcular Calcular su transform transformada ada inversa inversa y represen representarla tarla.. Debido a la aritmética finita aparecen números complejos con parte imaginaria casi nula, para evitar este efecto tomaremos la parte real.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Fichero Fichero fourier_1.m fourier_1.m %% Tran Transf sfor orma mada da de Four Fourie ier r disc discre reta ta y su inve invers rsa a clear; a = 0; % Initial value b = 2*pi; % End value N = 2^11; % Total number of values t = a:(b-a a:(b-a)/( )/(N-1 N-1):b ):b; ; s = t.*sin t.*sin(4* (4*t); t); figure; subplot(3,1, subplot(3,1,1); 1); plot(t,s); plot(t,s); axis([a axis([a b min(s) min(s) max(s)]); max(s)]); title(’Orig title(’Original’) inal’); ; s_dft s_dft = fft(s) fft(s); ; % Transf Transform ormada ada Discre Discreta ta de Fourie Fourier r subplot(3,1,2); plot(t,abs(fftshift(s_dft))); title(’Trans title(’Transforma formada da discreta discreta de Fourier’); Fourier’);

106

Ejercicios

Original 5 0 −5 0

1

2 3 4 Transformada discreta de Fourier

5

6

1

2 3 4 Reconstrucción Reconstrucción a partir de la DFT

5

6

5

6

3000 2000 1000 0 0 5 0 −5 0

1

2

3

4

Figura B.10: Señal original en el gráfico superior y debajo la reconstrucción después de haber calculado la transformada discreta de Fourier y su inversa.

axis([a axis([a b -0.1*max(ab -0.1*max(abs(s_df s(s_dft)) t)) 1.1*max(abs 1.1*max(abs(s_dft (s_dft))]); ))]); s_rec s_rec = ifft(s ifft(s_df _dft); t); % Transf Transform ormada ada Invers Inversa a Discre Discreta ta de Fourie Fourier r s_re s_rec c = real real(s (s_r _rec ec); ); % Toma Tomamo mos s part parte e real real, , imag imagin inar ario ios s son son casi casi cero cero subplot(3,1, subplot(3,1,3); 3); plot(t,s_rec plot(t,s_rec); ); axis([a axis([a b min(s_rec) min(s_rec) max(s_rec)]) max(s_rec)]); ; title( title(’Re ’Recon constr strucc ucción ión a partir partir de la DFT’); DFT’); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Ejercicio B.28 Análisis de frecuencias de una suma de senos. El fichero de nombre fourier_2.m contiene las órdenes necesarias y el resultado se puede ver en la figura B.11. 1. Crear Crear una señal a partir partir la función función f ( f (t) = 2 co cos( s(44πt) πt ) + 0. 0.5 sen( sen(24 24πt πt)) en el intervalo [ 16 muestreada con 1024 puntos. 16,, 16] 16] muestreada

−


107

2. Calcular Calcular su transform transformada ada de Fourier ourier y realizar realizar su análisis análisis de frecuencias frecuencias relacionándolo con la función f ( f (t) utilizada para generar la función. Original 8 6 4 2 0 −2 −4 −6 −8

−15

−1 0

−5

0

5

10

15

10

15

Transformada discreta de Fourier 1000 800 600 400 200 0 −1 5

−1 0

−5

0

5

Figura B.11: Señal original en el gráfico superior y debajo los valores (en módulo) de la transformada discreta de Fourier.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Ficher Fichero o fourie fourier_2 r_2.m .m %% Anális Análisis is de frecue frecuenci ncias as con la transf transform ormada ada discre discreta ta de Fourie Fourier r clear; a = -16; % Initial value b = 16; % End value N = 2^10; % Total number of values t = a:(b-a a:(b-a)/( )/(N-1 N-1):b ):b; ; s = 2*cos(4*pi* 2*cos(4*pi*t)+0.5 t)+0.5*sin(2 *sin(24*pi*t 4*pi*t); ); figure; subplot(2,1, subplot(2,1,1); 1); plot(t,s); plot(t,s); axis([ axis([-16 -16 16 -8 8]); 8]); title( title(’Or ’Origi iginal nal’); ’); s_dft s_dft = fft(s) fft(s); ; % Transf Transform ormada ada Discre Discreta ta de Fourie Fourier r

108

Ejercicios

subplot(2,1,2); plot(t,abs(fftshift(s_dft))); title(’Trans title(’Transforma formada da discreta discreta de Fourier’); Fourier’); axis([-16 axis([-16 16 -0.1*max(abs -0.1*max(abs(s_dft (s_dft)) )) 1.1*max(abs 1.1*max(abs(s_dft (s_dft))]); ))]); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Ejercicio B.29 Análisis de frecuencias de una señal. 1. Crear Crear una señal a partir partir la función función g(t) =

1 + cos(24πt cos(24πt)) 1 + 4t 4t2

en el intervalo [ 16 16,, 16] 16] muestreada muestreada con 1024 puntos.

−

2. Calcular Calcular su transform transformada ada de Fourier ourier y realizar realizar su análisis análisis de frecuencias frecuencias relacionándolo con la función g (t) utilizada para generar la función. 3. ¿Puede ¿Puede obtenerse obtenerse tanta tanta información información como en el caso anterior? anterior?

B.7.2 B.7.2

El espect espectro ro de seña señales les wav wavele elets ts y scali scaling ng

Ejercicio B.30 Calcular y representar gráficamente el espectro de la primera scaling y la primera wavelet wavelet de nivel 1 de la coif2 (desplazar la transformada de Fourier si es necesario). El fichero spectrum.m contiene las órdenes necesarias y la figura B.12 el resultado. Espectro de la scaling

Espectro de la wavelet

3

3

2.5

2.5

2

2

1.5

1.5

1

1

0.5

0.5

0

0

−0.5

−0.5

−1 −0.5

0

0. 5

−1 −0.5

0

0.5

Figura B.12: Para la coif2 a la izquierda el espectro de la scaling v11 y a la derecha de la wavelet w11 .


109

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Ficher Fichero o spectr spectrum. um.m m %% Para construir el gráfico de los espectros de scaling y wavelets %%%%%%%% %% El DWT Extension Mode debe ser ’per’ %%%%%%%% clear; k = 10; N = 2^k; % Longitud n = 1; % Niveles w = ’coif2’; % Wavelet % Construimos L a partir de una señal de ceros s = zeros( zeros(1,N 1,N); ); [C,L] = wavedec(s,n wavedec(s,n,w); ,w); t = -.5: -.5:1/ 1/(N (N-1 -1): ):.5 .5; ; % Para Para gráf gráfic icos os C = zeros(1,N); % Anulamos i = 1; C(i) = 1; v = waverec(C,L,w); v_dft v_dft = fft(v) fft(v); ; v_spec v_spec = abs(v_ abs(v_dft dft).^ ).^2; 2; figure; subplot(1,2,1); plot(t,fftshift(v_spec)); title( title(’Es ’Espec pectro tro de la scalin scaling’) g’); ; axis axis([ ([-. -.5 5 .5 -1 3]); 3]); C = zeros(1,N); % Anulamos i = N/2+1; C(i) = 1; w = waverec(C,L,w); w_dft w_dft = fft(w) fft(w); ; w_spec w_spec = abs(w_ abs(w_dft dft).^ ).^2; 2; subplot(1,2,2); plot(t,fftshift(w_spec)); title( title(’Es ’Espec pectro tro de la wavele wavelet’) t’); ; axis axis([ ([-. -.5 5 .5 -1 3]); 3]); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Ejercicio B.31 Obtención y representación de espectros de wavelets y scaling de distintos niveles. 1. Calcular Calcular y represen representar tar gráficamen gráficamente te los espectros espectros de las señales scaling primeras desde nivel 1 a 4 de la coif2.

110

Ejercicios

2. Calcular Calcular y represen representar tar gráficamen gráficamente te los espectros espectros de las señales señales wavel wavelet et primeras desde nivel 1 a 4 de la coif2.

B.7.3 B.7.3

Anális Análisis is de frecue frecuenci ncias as y transfo transforma rmada da wav wavelet elet

Veamos a continuación que el uso conjunto de la transformada wavelet y el análisis de frecuencias de la transformada de Fourier discreta puede mejorar nuestro análisis de frecuencias del ejercicio B.29.

Ejercicio B.32 Vamos a volver a realizar el análisis de frecuencia de la señal del ejercicio B.29 pero antes manipularemos la señal con la transformada wavelet. Las órdenes se pueden encontrar en el fichero fourier_wav.m y el resultado en la figura B.13. 1. Crear una señal señal que represente represente la función g(t) =

1 + cos(24πt cos(24πt)) 1 + 4t 4t2

en el intervalo [ 16 16,, 16] 16] con con 1024 puntos.

−

2. Obtener Obtener su primera señal promedio promedio A1 con la wavelet coif2. 3. Obtener Obtener su primera primera señal de detalle D 1 con la wavelet coif2. 4. Realizar Realizar el análisis análisis de frecuencia frecuenciass con la transform transformada ada discreta discreta de Fourier de la señal original y de las señales promedio y detalle. Obtener conclusiones (desplazar la transformada si es necesario).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Fichero Fichero fourier_wav fourier_wav.m .m %% Anális Análisis is de frecue frecuenci ncias as con la transf transform ormada ada discre discreta ta de Fourie Fourier r %% combin combinada ada con wavele wavelets ts clear; a = -16; % Initial value b = 16; % End value N = 2^10; % Total number of values t = a:(b-a a:(b-a)/( )/(N-1 N-1):b ):b; ; s = (1+2*cos(24 (1+2*cos(24*pi*t) *pi*t))./(1+ )./(1+4*t.^2 4*t.^2); );


111

Original s

DFT de s

5

40

0

20

−5

0 −10

0 1 Aprox. A

10

5

− 10

0 1 DFT de A

10

− 10

0 DFT de D 1

10

− 10

0

10

40

0

20

−5

0 − 10

0 Detalle D 1

10

5

40

0

20

−5

0 − 10

0

10

Figura B.13: En la columna de la izquierda tenemos la señal original, la aproximación A1 y el detalle D 1 . En la columna de la derecha se puede ver sus correspondientes transformadas discretas de Fourier (valores absolutos).

w = ’coif2’; % Wavelet n = 1; % Niveles figure; subplot(3,2, subplot(3,2,1); 1); plot(t,s); plot(t,s); axis([ axis([a a b -8 8]); 8]); title( title(’Or ’Origi iginal nal s’); s’); s_dft s_dft = fft(s) fft(s); ; % Transf Transform ormada ada Discre Discreta ta de Fourie Fourier r subplot(3,2,2); plot(t,abs(fftshift(s_dft))); title(’DFT de s’); axis([a axis([a b -0.1*max(ab -0.1*max(abs(s_df s(s_dft)) t)) 1.1*max(abs 1.1*max(abs(s_dft (s_dft))]); ))]); % Descomposic Descomposición ión wavelet wavelet [C,L] = wavedec(s,n wavedec(s,n,w); ,w); A1 = wrcoef(’a’, wrcoef(’a’,C,L,w, C,L,w,1); 1); D1 = wrcoef(’d’, wrcoef(’d’,C,L,w, C,L,w,1); 1); subplot(3,2, subplot(3,2,3); 3); plot(t,A1); plot(t,A1);

112

Ejercicios

axis([ axis([a a b -8 8]); 8]); title( title(’Ap ’Aprox rox. . A^1’); A^1’); A1_d A1_dft ft = fft( fft(A1 A1); ); % Tran Transf sfor orma mada da Disc Discre reta ta de Four Fourie ier r de A^1 A^1 subplot(3,2,4); plot(t,abs(fftshift(A1_dft))); title( title(’DF ’DFT T de A^1’); A^1’); axis([a axis([a b -0.1*max(ab -0.1*max(abs(A1_d s(A1_dft)) ft)) 1.1*max(abs 1.1*max(abs(A1_df (A1_dft))]); t))]); subplot(3,2, subplot(3,2,5); 5); plot(t,D1); plot(t,D1); axis([ axis([a a b -8 8]); 8]); title( title(’De ’Detal talle le D^1’); D^1’); D1_d D1_dft ft = fft( fft(D1 D1); ); % Tran Transf sfor orma mada da Disc Discre reta ta de Four Fourie ier r de D^1 D^1 subplot(3,2,6); plot(t,abs(fftshift(D1_dft))); title( title(’DF ’DFT T de D^1’); D^1’); axis([a axis([a b -0.1*max(ab -0.1*max(abs(D1_d s(D1_dft)) ft)) 1.1*max(abs 1.1*max(abs(D1_df (D1_dft))]); t))]); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

B.8 B.8

Dete Detecc cció ión n de sub subseñ señal ales es de de cort corta a durac duració ión n

Veamos que a través de la correlación de señales se pueden detectar señales de “corta duración” en señales “grandes”. Esta técnica puede ser utilizada para detectar irregularidades o patrones en señales (por ejemplo electrocardiogramas). Finalmente veremos que el tratamiento de las señales con wavelets puede mejorar la detección. Dadas dos señales f y g , se define la señal de correlación (f : g ) cuyo k-ésimo valor es (f : g) g )k := f := f 1 gk + f 2 gk+1 +

· · · + f N N gk+N −1.

Aquí hay que precisar que tomamos g j+ cada j (es decir, hacemos j +N = g j para cada j g periódica). La correlación normalizada por f es 1 (f : g) g ). E (f ) Supongamos que f es f es una señal formada en los primeros puntos por una porción (subseñal) de g, siendo 0 en el resto de puntos. Partimos también del supuesto que la porción que hemos tomado de g tiene una energía superior

B.8 Detección

113

a cualquier otra porción de g de la misma longitud. A partir de la desigualdad de Cauchy para productos escalares se deduce que el valor máximo de la correlación coincide justo cuando estamos al principio del soporte de la porción tomada de g de g.. Al normalizar dicho valor máximo es 1 es 1.. Esta situación puede darse, por ejemplo, cuando pretendemos localizar una anomalía cardíaca en un electrocardiograma. Para una detección más clara de la posición de la subseñal en la señal g , podemos realizar previamente un filtrado mediante el MRA asociado a cierta wavelet. A un determinado nivel m, como g = A = A m + Dm +

· · · + D1 ,

tomamos g˜ = g = g

− Am = D m + · · · + D1.

Hacemos lo propio con f con f ,, produciendo f . f ˜. De esta manera hemos filtrado permitiendo el paso alto y paso banda, que preserva las características de la subseñal y reduce posibles interferencias. Finalmente evaluamos la correlación normalizada 1 ˜ (f : g˜). ˜ E (f ) f ) El modo en que contribuye la DFT (o, más bien, la FFT) a este proceso es de la siguiente manera. Se cumple la relación entre correlación y DFT dada por (  f : g) g ) = f ˆ f ˆg. gˆ. Para reducir el número de operaciones procedemos de este modo: 1. Comput Computamo amoss la DFT de f y g (mediante FFT). 2. Multiplica Multiplicamos mos f ˆ y gˆ. 3. Computamo Computamoss la DFT inversa inversa de f ˆ f ˆgˆ y obtenemos (f : g) g ). Esto requiere requiere O(N log(N log(N )) )) operac operacio ione nes, s, que que hay hay que que comp compar arar ar con con las las O(N 2 ) operaciones si hiciéramos la correlación directamente.

Ejercicio B.33 Program Programar ar una funció función n para para MA MATLA TLAB B llamad llamadaa corr que acepte como argumentos dos señales de la misma longitud como vectores fila y como resultado calcule la correlación entre las dos señales. Un ejemplo de cómo puede realizarse esto se muestra en el fichero corr.m.

114

Ejercicios

Archivo: Archivo: corr.m functi function on [z] = corr(x corr(x,y) ,y) %CORR Computes the correlation of x and y % z = CORR CORR(x (x,y ,y) ) crea create tes s a vect vector or z with with the corre correla lati tion on % of vec vectors x and y. % Vect Vector ors s x and and y MUST UST BE vect vector or file files. s. z = zeros( zeros(siz size(x e(x)); )); for i=1:length(z i=1:length(z), ), z(i) = dot(x,circs dot(x,circshift(y hift(y,[1 ,[1 -i+1])); -i+1])); end;

Ejercicio B.34 Supongamos que se tienen dos señales, una de un electrocardiograma real y la otra es el aislamiento de una arritmia cardíaca. Intentaremos detectar y localizar si el electrocardiograma tiene la arritmia. El fichero detect.m contiene las órdenes necesarias para realizar este ejercicio. El resultado gráfico se puede ver en la figura B.14. 1. Obtener Obtener la señal que simula simula un electrocardio electrocardiograma grama de ecg.mat. Representarlo. 2. Obtene Obtenerr la señal señal que simula simula una arritm arritmia ia de arritmia.mat. Representarla. 3. Calcul Calcular ar la correl correlaci ación ón entre entre la arritm arritmia ia y el electroc electrocard ardiog iogram ramaa para para detectar detectar el latido latido “defectuos “defectuoso”. o”. Represen Representar tar la correlación correlación.. 4. Para Para detectar mejor mejor esta anomalía anomalía se propone normalizar normalizar la correlación, correlación, eliminar los valores negativos y tomar el cuadrado de los valores restantes. Representar la señal resultante.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Ficher Fichero o detect detect.m .m %% Detecc Detección ión de patron patrones es cortos cortos en señale señales s grande grandes s clear;

B.8 Detección

115

Electrocardiograma 1 0 −1 200

400

600

800

1000 1200 Arritmia

1400

1600

1800

2000

200

400 600 800 1000 1200 1400 1600 Correlación arritmia y electrocardiogra electrocardiograma ma

1800

2000

200

400

600 800 1000 1200 1400 Valores positivos al cuadrado

1600

1800

2000

200

400

600

1600

1800

2000

1 0 −1

50 0 −50

40 20 0

800

1000

1200

1400

Figura B.14: De arriba a abajo tenemos una simulación de una electrocardiograma, una simulación de una arritmia, la correlación entre ambas y para resultar más el efecto de detección hemos representado el cuadrado de tan sólo los valores positivos de la correlación.

load load ecg; ecg; load arritmia; arritmia; figure; subplot(4,1, subplot(4,1,1); 1); plot(ecg); plot(ecg); axis([1 axis([1 length(ecg) length(ecg) -2*max(ecg) -2*max(ecg) 2*max(ecg)] 2*max(ecg)]); ); title(’Electrocardiograma’); subplot(4,1,2); plot(arritmia); axis([1 axis([1 length(arri length(arritmia) tmia) -2*max(arri -2*max(arritmia) tmia) 2*max(arritm 2*max(arritmia)]) ia)]); ; title(’Arritmia’); a_corr_e a_corr_e = corr(arritm corr(arritmia,ecg ia,ecg); ); % Correlación Correlación subplot(4,1,3); plot(a_corr_e); axis([1 axis([1 length(a_co length(a_corr_e) rr_e) -2*max(a_co -2*max(a_corr_e) rr_e) 2*max(a_corr 2*max(a_corr_e)]) _e)]); ; title(’Corre title(’Correlació lación n arritmia arritmia y electrocardi electrocardiograma ograma’); ’);

116

Ejercicios

final = (max(a_corr (max(a_corr_e,0). _e,0).^2)/en ^2)/energy( ergy(arritm arritmia); ia); % Positivos Positivos y cuadrado cuadrado subplot(4,1, subplot(4,1,4); 4); plot(final); plot(final); axis([1 axis([1 length(fina length(final) l) -1 1.1*max(fin 1.1*max(final)]); al)]); title(’Valor title(’Valores es positivos positivos al cuadrado’); cuadrado’); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Ejercicio B.35 Repetir el ejercicio anterior con las señales ecg_big.mat y arritmia_big.mat. Observar que el precio computacional al calcular la corre-

lación entre señales con un número grande de datos es muy alto. Si todavía se es escéptico en este punto se propone repetir el proceso tomando las señales ecg_big_big.mat y arritmia_big_big.mat. A continuación vamos a dar un método alternativo para calcular la correlación entre dos señales que no sea tan computacionalmente caro. Comenzaremos programando una función alternativa para calcular la correlación que haga uso de transformada discreta de Fourier.

Ejercicio B.36 Repetir el ejercicio B.34 utilizando el siguiente “atajo”, en el sentido de rapidez computacional. 1. Obtene Obtenerr las señale señaless ecg_big.mat y arritmia_big.mat. 2. Calcular Calcular las transform transformadas adas discretas discretas de Fourier ourier de las dos señales señales (usar (usar la FFT). 3. Multiplica Multiplicarr término a término término la señal conjugada conjugada de la transform transformada ada de la arritmia_big con la transformada del ecg_big. 4. Calcular Calcular la transform transformada ada inversa inversa de Fourier ourier (usar (usar la IFFT) IFFT) del resultado (es conveniente eliminar las partes imaginarias que aparecen por los problemas de trabajar con aritmética finita). 5. Compro Comprobar bar que se obtien obtienee lo mismo mismo que si usamos usamos la correl correlaci ación ón entre entre la arritmia_big y el ecg_big directamente. directamente.

Ejercicio B.37 Programar una función para MATLAB llamada corrfft que acepte como argumento dos señales de la misma longitud como vectores fila y como resultado calcule la correlación entre las dos señales. Para calcular la correlación la función hará uso de la relación existente entre la transformada

B.8 Detección

117

discreta de Fourier y la correlación, es decir, primero deberá calcular las transformadas discretas de Fourier de las señales y después multiplicar término a término la conjugada de la primera transformada con la segunda transformada. Finalmente se volverá a calcular la transformada discreta inversa de Fourier del resultado. Un ejemplo de cómo puede realizarse esto se muestra en el fichero corrfft.m. Archivo: Archivo: corrfft.m corrfft.m functi function on [z] = corrff corrfft(x t(x,y) ,y) %CORRFFT Computes the correlation of x and y % z = CORR CORRFF FFT( T(x, x,y) y) crea create tes s a vect vector or z with with % the the corr correl elat atio ion n of vect vector ors s x and and y. % Vect Vector ors s x and and y MUST UST BE vect vector or file files. s. % This This funct functio ion n uses uses FFT and IFFT IFFT to compu compute te % the the corr correl elat atio ion. n. z = real(ifft(co real(ifft(conj(ff nj(fft(x)). t(x)).*fft(y *fft(y))); )));

Ejercicio B.38 Es importante comprobar la mejora computacional que supone calcular la correlación haciendo uso de la transformada discreta de Fourier. En los ejercicios B.33 y B.37 hemos automatizado el cálculo de la correlación programando funciones para MATLAB, ahora podemos calcular correlaciones por ambos métodos observando el tiempo que necesitamos para ello. Hacer pruebas con las señales arritmia y ecg, con las señales arritmia_big y ecg_big y finalmente con las señales arritmia_big_big y ecg_big_big. Sin problemas se podrá observar la diferencia en el tiempo de cálculo de la correlación a partir de la definición y haciendo uso de la transformada discreta de Fourier.

Ejercicio B.39 En este ejercicio vamos a ver que los resultados de la correlación entre señales también se pueden mejorar si añadimos componentes de la transformada wavelet de señales. 1. Obtene Obtenerr las señale señaless ecg.mat y arritmia.mat. 2. Con la wave wavelet let coif5 obtener el promedio A4 de ambas señales. Restar este promedio a cada seña, es decir, tomaremos D 1 + D2 + D3 + D4 .

118

Ejercicios

3. Calcular Calcular la correlación correlación de las señales detalle obtenidas. obtenidas. 4. Normalizar Normalizar,, eliminar eliminar negativo negativoss y elevar elevar al cuadrado. cuadrado. Representar Representar gráfigráficamente lo obtenido. El fichero detect_wav contiene las órdenes necesarias para realizar el ejercicio. El resultado que se muestra en la figura B.15, se debe comparar con lo obtenido en el ejercicio B.34. 1

2

3

4

Detalles D +D +D +D del electrocardiogram electrocardiograma a 1 0 −1 200

400

1600 600 800 1000 1200 1400 1 2 3 4 Detalles D +D +D +D de la arritmia

1800

2000

1800

2000

2 0 −2

200

400 600 800 1000 1200 1400 1600 Correlación detalles arritmia y electrocardiogram electrocardiograma a

40 20 0 −20 −40 200

400

600 800 1000 1200 1400 Valores positivos al cuadrado

1600

1800

2000

200

400

600

1600

1800

2000

20 10 0 800

1000

1200

1400

Figura B.15: De arriba a abajo tenemos: los detalles D1 + + D 4 de la simulación de un electrocardiograma; los mismos detalles de la simulación de una arritmia; la correlación entre ambas y, finalmente, para resaltar más el efecto de detección hemos representado el cuadrado de tan sólo los valores positivos de la correlación.

·· · · ·

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Fichero Fichero detect_wav. detect_wav.m m %% Detecc Detección ión de patron patrones es cortos cortos en señale señales s grande grandes. s. %% Combin Combina a correl correlaci ación ón con técnic técnicas as wavele wavelets ts clear; load arritmia; arritmia;

B.8 Detección

119

load load ecg; ecg; w = ’coif5’; % Wavelet n = 4; % Niveles % Extraemos detalles 1 a 4 de arritmia [C_arrit,L_a [C_arrit,L_arrit] rrit] = wavedec(arri wavedec(arritmia, tmia,n,w); n,w); A4_arrit A4_arrit = wrcoef(’a’, wrcoef(’a’,C_arri C_arrit,L_ar t,L_arrit,w, rit,w,n); n); D1_4_arrit D1_4_arrit = arritmia-A4_ arritmia-A4_arrit; arrit; % Extraemos detalles 1 a 4 del electro. [C_ecg,L_ecg [C_ecg,L_ecg] ] = wavedec(ecg, wavedec(ecg,n,w); n,w); A4_ecg = wrcoef(’a’, wrcoef(’a’,C_ecg, C_ecg,L_ecg, L_ecg,w,n); w,n); D1_4_ecg = ecg - A4_ecg; figure; subplot(4,1,1); plot(D1_4_ecg); axis([1 axis([1 length(D1_4 length(D1_4_ecg) _ecg) -2*max(D1_4 -2*max(D1_4_ecg) _ecg) 2*max(D1_4_e 2*max(D1_4_ecg)]) cg)]); ; title(’Detal title(’Detalles les D^1+D^2+D^3+ D^1+D^2+D^3+D^4 D^4 del electrocardi electrocardiograma ograma’); ’); subplot(4,1,2); plot(D1_4_arrit); axis([1 axis([1 length(D1_4 length(D1_4_arrit _arrit) ) -2*max(D1_4 -2*max(D1_4_arrit _arrit) ) 2*max(D1_4_ 2*max(D1_4_arrit) arrit)]); ]); title(’Detal title(’Detalles les D^1+D^2+D^3+ D^1+D^2+D^3+D^4 D^4 de la arritmia’); arritmia’); a_corr_e a_corr_e = corrfft(D1_ corrfft(D1_4_arri 4_arrit,D1_4 t,D1_4_ecg); _ecg); % Correlación Correlación (usa fft) subplot(4,1,3); plot(a_corr_e); axis([1 axis([1 length(a_co length(a_corr_e) rr_e) -2*max(a_co -2*max(a_corr_e) rr_e) 2*max(a_corr 2*max(a_corr_e)]) _e)]); ; title(’Corre title(’Correlació lación n detalles detalles arritmia arritmia y electrocardi electrocardiograma ograma’); ’); % Positi Positivos vos y cuadra cuadrado do final = (max(a_corr (max(a_corr_e,0). _e,0).^2)/en ^2)/energy( ergy(D1_4_a D1_4_arrit); rrit); subplot(4,1, subplot(4,1,4); 4); plot(final); plot(final); axis([1 axis([1 length(fina length(final) l) -1 1.1*max(fin 1.1*max(final)]); al)]); title(’Valor title(’Valores es positivos positivos al cuadrado’); cuadrado’); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

120

B.9 B.9

Ejercicios

Wavelets elets pac packets ets

El análisis con transformada wavelet packet lo haremos únicamente desde la herramienta gráfica. Debemos iniciar por tanto el entorno gráfico del menú principal con la orden wavemenu y elegir “Wavelet Packet 1-D”.

B.9. B.9.1 1

Anal Analiz izan ando do un una a seña señall

Cargar la señal almacenada en el fichero sumlichr.mat que proporciona la “toolbox”. Aparecerá la señal y ahora se debe elegir la wavelet wavelet a utilizar y los niveles de transformada. Tomar la db2 y 4 niveles. Es importante el método de cálculo de la entropía de la señal. En función del método elegido se puede optimizar la descomposición de la señal. Buscar en la ayuda wentropy para más datos. De momento tomaremos la entropía threshold con parámetro 1. Una vez seleccionadas estas opciones podemos analizar al señal. Aparece Aparece un árbol con la descomposició descomposición n de la señal. señal. Investiga Investigarr las distintas distintas opciones de representación y el método para representar cada subseñal. Usualmente no es necesario llegar al último nivel en todas las ramas del árbol. La estructura de descomposición se optimiza aplicando el método de entropía seleccionado. La descomposición que optimiza la entropía se obtiene con “besttree”. Observar que, con la elección de entropía tomada en este caso, obtenemos casi la descomposición inicial.

Ejercicio B.40 Inve Investi stigar gar las demás demás selecci seleccione oness de entro entropía pía que ofrece ofrece la herramienta gráfica prestando atención a los posibles cambios que puedan producirse en el cálculo de la mejor descomposición.

B.9. B.9.2 2

Comp Compri rimi mien endo do la seña señall

Seleccionando la opción de comprimir se pasa a la herramienta de compresión. La “toolbox” nos ofrece un valor para el umbral; este valor puede tomarse o cambiarse por otro. Podemos también observar la energía retenida y el número de ceros que tendríamos en la transformada.

Ejercicio B.41 Cambiar el valor valor del umbral a 0.8938 y observar los resultados de la compresión. Experimentar con otros valores. Una vez decidido el método de “threshold” y sus valores correspondientes se puede comprimir la señal seleccionando Compress. Al salir de la herramienta de

B.10 Transformada wavelet continua

121

compresión el programa nos pregunta si queremos actualizar la señal original con la versión comprimida, esto es, con la reconstrucción que se hace utilizando la señal transformada que hemos modificado.

B.9. B.9.3 3

Redu Reducc cció ión n de ru ruid ido o

Desde la ventana de análisis de señales con “wavelets packets” cambiaremos el método de cálculo de entropía a sure. Usando este método el umbral a poner sería 4.2975. A continuación podemos analizar la señal. Para reducir el ruido es conveniente seleccionar la mejor descomposición ya que la reducción de ruido suele ser más eficiente. Seleccionamos la herramienta de eliminación de ruido. Diferentes opciones pueden cambiarse si las que ofrece el programa no satisfacen nuestras necesidades. Por ejemplo, con valor 3.331 para el “Global threshold” hemos eliminado tanto “ruido” que perdimos parte esencial de la señal. De nuevo, al salir de la herramienta de eliminación de ruido el programa nos pregunta si queremos actualizar la señal original con la nueva versión que acabamos de obtener.

B.10 B.10

Transfor ransformad mada a wa wavelet contin continua ua

A continuación se describen las herramientas de la “Wavelet Toolbox” de MATLAB para realizar un análisis continuo de señales unidimensionales, tanto desde el entorno gráfico como desde la línea de comandos. Describiremos fundamentalmente las siguientes cuestiones.

• Realizar una transformada wavelet continua de la señal. • Representación gráfica de los coeficientes de la transformada continua. • Elección de las escalas a analizar. B.10 B.10.1 .1

Trans ransfo form rmad ada a wavelet elet con contin tinua desd desde e la líne línea a de cocomandos

Para realizar la transformada wavelet continua de una señal 1D sólo se necesita la orden cwt. Veamos su funcionamiento con un ejemplo.

122

Ejercicios

Consideremos la señal generada a partir de la función 2

100π (x−0.2) f ( f (x) = sin(40πx sin(40πx))e−100π + 2

50π (x−0.5) (sin(40πx (sin(40πx)) + 2 cos(16 cos(1600πx)) πx)) e−50π +

(B.1)

2

100π (x−0.8) 2 sin(16 sin(1600πx) πx)e−100π

en el intervalo [0 intervalo [0,, 1] realizando 1] realizando un muestreo con 2048 puntos. Supongamos que la señal ha sido generada y almacenada en la variable s. w = ’mexh’; % Fijamos Wavelet escalas = 1:24; % Escalas C = cwt(s, cwt(s,esc escala alas,w s,w); ); % Calcul Calculamo amos s la transf transform ormada ada

En el argumento de cwt hay que especificar la señal a analizar, las escalas (rango y paso) y el tipo de wavelet a utilizar. El resultado es una matriz que contiene los coeficientes de la transformada. En este ejemplo, la señal está formada por 2048 puntos, por lo que C es una matriz 24 2048 2048,, y cada fila corresponde a los coeficientes para una escala determinada. El comando cwt admite un argumento adicional que cuando está presente produce una representación gráfica de los valores absolutos de los coeficientes de la transformada wavelet continua.

×

C = cwt(s,escal cwt(s,escalas,w,‘ as,w,‘plot’) plot’); ; (véase la figura B.16)

Las normas para elegir las escalas son las que se indican a continuación. Haremos notar que son casi las mismas que en el entorno gráfico que veremos en la siguiente sección, aunque en el caso gráfico no se permite elegir la escala en orden decreciente.

• Todas las escalas deben ser números reales positivos. • El incremento de escala puede ser negativo o positivo (solo positivo cuando estamos con la herramienta gráfica).

puede ser mayor mayor de cierta cierta cantidad cantidad que depende de • La escala máxima no puede la señal a analizar (el entorno gráfico sí nos indica este valor máximo).

Algunos ejemplos válidos para definir las escalas pueden observarse a continuación. Es conveniente definirlos en una variable para su rápido cambio y actualización de los cálculos.

B.10 Transformada wavelet continua

123

3 2 1 0 −1 −2 −3

0

0 .1

0.2

0 .3

0 .4

0 .5

0.6

0 .7

0 .8

0.9

1

Absolute Values Values of Ca,b Coefficients Coefficients for a = 1 2 3 4 5 ... 23 21 19 17 a 15 s e 13 l a c 11 s 9 7 5 3 1 200

400

600

800

1000 1200 time (or space) b

1400

1600

1800

2000

Figura B.16: Señal f ( f (x); escalograma para escalas 1 a 24 con paso de 1.

escalas escalas escalas escalas

= 1:64; % Escalas 1 a 64 de 1 en 1 = 2:2:128; % Escalas 2 a 128 con paso 2 = [2 8 16 32:2:64]; % Escalas 2, 8, 16 y de 32 a 64 con paso 2 = 64:-2:2; % Escalas 64 a 2 con paso decreciente de 2

Ejercicio B.42 Cargar la señal noissin.mat que tiene ya generada la “toolbox”. 1. Represen Representar tar gráficament gráficamentee la señal. 2. Calcular Calcular y represen representar tar la transform transformada ada wavelet wavelet continu continuaa de la señal para las las esca escala lass 1, 2, . . . , 48. 48. 3. Calcular Calcular y represen representar tar la transform transformada ada wavelet wavelet continu continuaa de la señal para las las escal escalas as 2, 4, 4, . . . , 128. 128. 4. Observ Observar como a escalas escalas altas se destaca la periodici p eriodicidad dad de la señal.

B.10.2 B.10.2

Transformad ransformada a wave wavelet let contin continua ua en el ent entorno orno gráfico gráfico

El inicio del uso de la herramien herramienta ta gráfica que la “toolbox” “toolbox” tiene incorporada incorporada es similar a como se explicó en la sección A.2.2 para la transformada discreta. En

124

Ejercicios

primer lugar debemos cargar el menú general de la herramienta gráfica con la orden wavemenu. A continuación debemos elegir la opción “Continuous Wavelet 1-D”. Como ocurría con el entorno gráfico de la “Discrete Wavelet 1-D” la señal debe estar previamente grabada en un fichero *.mat que se carga de la forma habitual.

Ejercicio B.43 Analizar el ejemplo de la sección anterior (señal generada a partir de la función dada en ecuación (B.1)) utilizando la herramienta gráfica. Investigar las distintas maneras de colorear los coeficientes de la transformada (“Coloration Mode”) y observar, por ejemplo, que con una posible opción consiste en asignar colores entre los valores mínimo y máximo de los coeficientes. Otra posibilidad es hacerlo entre cero y el máximo del valor absoluto de los coeficientes, de este modo se obtiene una representación de los valores absolutos de la transformada.

Ejercicio B.44 Repetir el ejercicio B.42 desde la herramienta gráfica. Practicar nuevas posibilidades y opciones de la herramienta gráfica.

Apéndice C

Ejemplos y aplicaciones En este apéndice enumeraremos una serie de ejemplos y aplicaciones en los que el tratamiento con wavelets es de gran utilidad. Los cálculos se pueden hacer desde la herramienta gráfica o también desde la línea de comandos. Estos ejemplos se encuentran planteados en la documentación de la “Wavelet Toolbox for MATLAB” [5] donde, además, podemos encontrar otras aplicaciones.

C.1 C.1

Detec Detecci ción ón de “bre “break akdo down wn points” points”

Pretendemos encontrar cambios abruptos, cambios de forma, discontinuidades en las derivadas derivadas (primera (primera,, segunda, segunda, . . . ), etc. Las líneas generales generales para hacerlo hacerlo son: 1. Para Para encontrar encontrar una discontin discontinuida uidad d en la señal usar la wavel wavelet et de haar. 2. Para Para encon encontra trarr una una ruptur rupturaa en la deriv derivada ada j -ésima -ésima usar usar una wave wavelet let suficientemente regular para que se anulen al menos los primeros j momentos.

C.1. C.1.1 1

Un cam cambi bio o en en la frec frecue uenc ncia ia

Analizar la señal freqbrk.mat y obtener las conclusiones posibles en lo que a puntos de ruptura o cambio de frecuencia se refiere. Justificación. Veamos como el análisis wavelet aporta información sobre el instante en que varía la frecuencia de una señal. Esta información no la puede proporcionar la transformada de Fourier. 125

126

Ejemplos y aplicaciones

La señal pertenece al muestreo de la función g (x) =



sin(0. sin(0.03 x) sin(0. sin(0.3 x)

1 x 500 501 x 1000

≤ ≤ ≤ ≤

El fichero freqbrk_analisis.m puede ser utilizado para analizar la señal. Esta señal también se puede analizar desde la herramienta gráfica aunque incluiremos aquí el fichero por si puede ser útil para posteriores usos. El gráfico que se genera al ejecutar el fichero puede verse en la figura C.1 Original 1 0 −1 0 2 0 −2 0 2 0 −2 0 2 0 −2 0 2 0 −2 0

500 Aproximación 1

500 Aproximación 2

500 Aproximación 3

500 Aproximación 4

500 Aproximación 5

2 0 −2 0

500

Reconstruida

1000

1 0 −1

0

500 Detalle 1

1000

500 Detalle 2

1000

500 Detalle 3

1000

0

500 Detalle 4

1000

0

500 Detalle 5

1000

500

1000

1000

0.5 0 −0.5 0

1000

0.5 0 −0.5 0

1000

1000

1000

1 0 −1 2 0 −2

0.2 0 −0.2 0

Figura C.1: Análisis wavelet completo de la señal freqbrk con la db5 y 5 niveles. Puede observarse que en el detalle D 1 es en el que mejor se observa el cambio de frecuencia.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Fichero Fichero break_anal. break_anal.m m %% Ficher Fichero o para para analiz analizar ar la señal señal freqbr freqbrk.m k.mat at %% Buscam Buscamos os ‘‘brea ‘‘breakdo kdown wn points points’’ ’’ clear; load freqbrk.mat freqbrk.mat; ; s = freqbrk;

C.1 Detección de “breakdown p oints”

127

ls = length length(s) (s); ; w = ’db5’; % Wavelet n = 5; % Niveles [C,L] [C,L] = wavede wavedec(s c(s,n, ,n,w); w); % Descom Descompos posici ición ón for i = 1:n, % Iterador para automatizar A(i,:) = wrcoef(’a’,C wrcoef(’a’,C,L,w,i ,L,w,i); ); % Aproximacion Aproximaciones es D(i,:) D(i,:) = wrcoef wrcoef(’d (’d’,C ’,C,L, ,L,w,i w,i); ); % Detall Detalles es end; sr = wavere waverec(C c(C,L, ,L,w); w); % Recons Reconstru trucci cción ón figure figure; ; % Gráfic Gráficos os subplot(n+1, subplot(n+1,2,1); 2,1); plot(s); plot(s); title(’Orig title(’Original’) inal’); ; subplot(n+1, subplot(n+1,2,2); 2,2); plot(sr); plot(sr); title(’Reco title(’Reconstrui nstruida’); da’); for i = 1:n, % Iterador para automatizar subplot(n+1,2,1+2*i); plot(A(i,:)); [tit_a,err] [tit_a,err] = sprintf(’Apr sprintf(’Aproxima oximación ción %d’,i); %d’,i); title(tit_a); subplot(n+1,2,2+2*i); plot(D(i,:)); [tit_d,err] [tit_d,err] = sprintf(’Det sprintf(’Detalle alle %d’,i); %d’,i); title(tit_d); end; disc disc = find(D find(D(1, (1,:)= :)==ma =max(D x(D(1, (1,:)) :))); ); % Índice Índice con discon discontin tinuid uidad ad sprintf(’Índ sprintf(’Índice ice de la discontinuid discontinuidad: ad: %d’,disc) %d’,disc) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Los detalles D detalles D 1 y D 2 muestran claramente la discontinuidad. Para localizar el punto buscamos la posición del máximo del detalle D detalle D 1 . Por medio de la orden find obtendremos fácilmente que el máximo se alcanza en el índice 500. disc = find(D(1,:)= find(D(1,:)==max(D =max(D(1,:)) (1,:))); );

Las señales de detalles D 3 y D4 contienen el seno de menor frecuencia mientras que el seno de mayor frecuencia queda perfectamente aislado en A5 .

128

C.1. C.1.2 2


Evol Evoluc ució ión n a largo largo tiem tiempo po

Analizar la señal wnoislop.mat y obtener las conclusiones posibles en lo que a puntos de ruptura se refiere y también identificar la evolución a largo tiempo de la señal. Justificación. La señal pertenece a una “rampa” a la que se le ha añadido tanto ruido que no permite distinguir la evolución de la señal principal. Concretamente pertenece al muestreo de la función g (t) =



3t 500

+ ruido 3 + ruido

1 t 499 500 t 1000

≤ ≤ ≤ ≤

Para analiza analizarla rla se recomi recomiend endaa utiliz utilizar ar la wavele avelett ’db3’ y 6 niveles. El gráfico que se genera al ejecutar el fichero evol_anal.m puede verse en la figura C.2. Observar fundamentalmente como el ruido (frecuencias altas) queda aislado en los detalles. Afinando un poco se puede incluso identificar el cambio en la rampa alrededor del valor 500. Original 4 2 0 −2 0

200

400

600

800

1000

800

1000

800

1000

6

Aproximación A 5 0 −5 0

200

400

600 1

6

Detalles D + . . +D 1 0 −1 0

200

400

600

Figura C.2: Análisis wavelet parcial de la señal wnoislop con la db3 y 6 niveles. Puede observarse en la aproximación la evolución en el tiempo de la señal.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

C.1 Detección de “breakdown p oints”

129

%% Fichero Fichero evol_anal.m evol_anal.m %% Ficher Fichero o para para analiz analizar ar la señal señal wnoisl wnoislop. op.m m %% Buscamos la evolución de la señal a largo tiempo clear; load wnoislop.ma wnoislop.mat; t; s = wnoislop; ls = length length(s) (s); ; w = ’db3’; % Wavelet n = 6; % Niveles [C,L] [C,L] = wavede wavedec(s c(s,n, ,n,w); w); % Descom Descompos posici ición ón for i = 1:n, % Iterador para automatizar A(i,:) = wrcoef(’a’,C wrcoef(’a’,C,L,w,i ,L,w,i); ); % Aproximacion Aproximaciones es D(i,:) D(i,:) = wrcoef wrcoef(’d (’d’,C ’,C,L, ,L,w,i w,i); ); % Detall Detalles es end; figure figure; ; % Gráfic Gráficos os subplot(3,1, subplot(3,1,1); 1); plot(s,); plot(s,); title(’Orig title(’Original’) inal’); ; subplot(3,1, subplot(3,1,2); 2); plot(A(6,:)) plot(A(6,:)); ; title(’Aprox title(’Aproximació imación n A^6’); subplot(3,1, subplot(3,1,3); 3); plot(s-A(6,: plot(s-A(6,:)); )); title(’Deta title(’Detalles lles D^1+ . . +D^6’); +D^6’); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

C.1.3 C.1.3

Discon Discontin tinuid uidad ad en en deri deriv vadas adas

Analizar la señal scddvbrk.mat y obtener las conclusiones posibles en lo que a puntos de discontinuidad se refiere (en la señal o en sus derivadas). Justificación. Para detectar discontinuidades en la derivada i-ésima de una señal hay que utilizar una wavelet con al menos i momentos nulos. Utilizaremos una ’db4’ y varios niveles de transformada. La señal considerada pertenece al muestreo de la función g(x) =



exp( 4t2 ) exp( t2 )

− −

−0.5 ≤ x ≤ 0 0 < x ≤ 0.5

que presenta una discontinuidad en su derivada segunda. Desde la herramienta gráfica y ampliando convenientemente se observa que los detalles tienen valores significativos en el centro de la señal y son despreciables en cualquier otro lugar. Esto sugiere que hay información sobre algún

130


tipo de cambio. Si repetimos el cálculo con la wavelet haar o incluso con la db2 se observa que los detalles en la discontinuidad valen cero. Esto es debido a que el análisis debe realizarse con wavelets con al menos 2 momentos nulos. Por comodidad para el lector incluimos el fichero disc_anal.m con el análisis de la señal y los parámetros de ampliación del gráfico necesarios para observar la discontinuidad en la derivada segunda (figura C.3). Original 1 0.5

1 0 −1 1 0 −1

100

200

300

400

500 600 1 Detalle D

700

800

900

1000

100

200

300

400

500 600 2 Detalle D

700

800

900

1000

100

200

300 400 500 600 700 1 Ampliación apropiada de D

800

900

1000

200

300 400 500 600 700 2 Ampliación apropiada de D

800

900

1000

200

300

800

900

1000

−6

1 0 −1

x 10

100 −6

x 10 2 0 −2

100

400

500

600

700

Figura C.3: Análisis wavelet de la señal scddvbrk con la db4 y 2 niveles. Puede observarse en la ampliación de los detalles la discontinuidad (en la segunda derivada) de la señal.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Ficher Fichero o disc_a disc_anal nal.m .m %% Ficher Fichero o para para analiz analizar ar la señal señal scddvb scddvbrk. rk.mat mat %% Buscam Buscamos os discon discontin tinuid uidade ades, s, en la señal señal o sus deriva derivadas das clear; load scddvbrk.ma scddvbrk.mat; t; s = scddvbrk; ls = leng length th(s (s); ); w = ’db4’; % Wavelet n = 2; % Niveles

C.2 Identificando frecuencias puras

131

[C,L] [C,L] = wavede wavedec(s c(s,n, ,n,w); w); % Descom Descompos posici ición ón for i = 1:n, % Iterador para automatizar A(i,:) = wrcoef(’a’,C wrcoef(’a’,C,L,w,i ,L,w,i); ); % Aproximacion Aproximaciones es D(i,:) D(i,:) = wrcoef wrcoef(’d (’d’,C ’,C,L, ,L,w,i w,i); ); % Detall Detalles es end; figure figure; ; % Gráfic Gráficos os subplot(5,1, subplot(5,1,1); 1); plot(s); plot(s); title(’Original’); axis axis([ ([1 1 ls 0.2 0.2 1.1] 1.1]); ); subplot(5,1, subplot(5,1,2); 2); plot(D(1,:)) plot(D(1,:)); ; title(’Detal title(’Detalle le D^1’); axis([1 ls -1 1]); subplot(5,1, subplot(5,1,3); 3); plot(D(2,:)) plot(D(2,:)); ; title(’Detal title(’Detalle le D^2’); axis([1 ls -1 1]); subplot(5,1, subplot(5,1,4); 4); plot(D(1,:)) plot(D(1,:)); ; title(’Ampli title(’Ampliación ación apropiada apropiada de D^1’); D^1’); axis([ axis([1 1 ls -10^(-10^(-6) 6) 10^(-6 10^(-6)]) )]); ; subplot(5,1, subplot(5,1,5); 5); plot(D(2,:)) plot(D(2,:)); ; title(’Ampli title(’Ampliación ación apropiada apropiada de D^2’); D^2’); axis([1 axis([1 ls -.3*10^(-5) -.3*10^(-5) .3*10^(-5)] .3*10^(-5)]); ); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

C.2 C.2

Iden Identific tifican ando do frec frecuen uenci cias as pura purass

Con este ejemplo pretendemos ver cómo el análisis con wavelets nos permite identificar lo que se conoce como una función tipo Fourier, esto es, descomposición de la función como suma de sinusoides de diferentes frecuencias. Analizar la señal sumsin.mat y obtener las conclusiones posibles en lo que se refiere a la descomposición de la señal en suma de sinusoides de diferentes frecuencias. Justificación. Primero diremos que la señal ha sido generada a partir de un muestreo de la función s(t) = sin(3 t) + sin(0. sin(0.3 t) + sin(0. sin(0.003 t). El análisis vamos a realizarlo con la wavelet db3 y 5 niveles de descomposición.

132


Por comodidad usaremos la herramienta gráfica (véase la figura C.4). Signal and Approximation(s)

Coefs, Signal and Detail(s) 5 4 3 2 1

2

s

0 −2

2

1

s

0

a

5

0

−2

−1

0.1 0 −0.1

1

a

4

3

a

2

d

0

4

−1

0

0.5

−2

0

2

−0.5

d

3

0.1

0

d

0 −2

2

−0.1 1

1

1

5

1

2

a

d

0 −1

a

cfs

0

d

0

1

−1 20 0

40 0

60 0

80 0

10 00

−1

20 0

40 0

60 0

80 0

10 0 0

Figura C.4: Análisis wavelet de la señal sumsin con la db3 y 5 niveles. Después de analizar la señal observamos que aparecen varias sinusoides claramente identificadas. Algunas de las que mejor pueden verse las encontramos en aproximación A4 , detalle D4 y detalle D1 . Corresponderían a sinusoides lenta, media y rápida respectivamente. Ampliaciones del detalle D1 nos permiten aislar un periodo y estimar que está compuesto de unas 10 oscilaciones. Este dato puede ser usado para estimar el periodo (respecto de los valores muestreados) con un resultado aproximado de 2. La señal detalle D 3 (y también quizás la D 4 ) está compuesta de las sinusoides de oscilación media. Observamos que de la aproximación A3 a la A4 ha sido eliminada la información de la sinusoide media. Por tanto usaremos las

C.3 Análisis de un caso real

133

aproximaciones de la A1 a la A3 para obtener información del periodo de la sinusoide de oscilación media. Ampliando A1 puede observarse que el periodo es aproximadamente 20. Sólo queda identificar el periodo de la sinusoide de oscilación lenta. Puede ser estimado a partir de la aproximación A4 . Sin dificultad se observa que la distancia entre dos máximos sucesivos es aproximadamente 200.

Ejercicio C.1 Construir (o modificar la anterior) otras señales formadas for sumas de sinusoides e intentar identificar los periodos de las mismas.

C.3 C.3

Aná Análisi lisiss de de un un cas caso o rea reall

En este ejemplo vamos a hacer un análisis más detallado de un caso práctico (referencia [6] de la bibliografía). La fuente original del artículo es Misi Misiti, ti, M.; M.; Misi Misiti ti,, Y.; Oppe Oppenh nhei eim, m, G.; G.; Poggi Poggi,, J.M. J.M., Décompositio en ondelettes et méthodes comparatives: étudy d’une courbe de charge électrique , Revue de Statistique Appliqué vol. XLII (1994), n. 2, pp. 57–77. Los datos pertenecen a un planta eléctrica. Consisten en los valores tomados minuto a minuto de la carga eléctrica consumida durante un periodo de cinco semanas. Se tiene un total de 50400 datos y en el fichero leleccum.mat hay una porción (4320 datos) de ellos. Algunas características de los datos deben ser comentados.

• Los datos finales se han calculado sumando cientos de valores tomados con sensores de medida y por lo tanto se han generado errores al medir.

términos os generale generaless el consum consumoo puede puede atribuir atribuirse se en un 50 % a la in• En términ dustri dustriaa y el otro otro 50 % al consum consumoo indivi individu dual. al.

• Hay más de 10 millones de consumidores individuales. • Los periodos fundamentales son diarios y semanales, ligados a la actividad económica.

• Los patrones de consumo diario también cambian con las horas (temporizadores de calefacción, etc.).

• Si en algún periodo los datos han “desaparecido”, éstos han sido reemplazados por otros.

134


• En algunos periodos de tiempo los sensores han fallado por encima de lo normal y los errores de media se han incrementado (entre datos 2400 y 3400).

La wavelet que vamos a utilizar es db3 y calcularemos transformadas hasta nivel 5 ya que 2 que 2 5 = 32 32 es es cercano a un periodo de media hora. Las dos primeras cuestiones corresponden al análisis de una parte de la señal. Las restantes al análisis de un periodo de tres días (señal completa). 1. Analizar Analizar los datos datos comprendi comprendidos dos entre entre [3600 : 3700] que 3700] que corresponden a una franja horaria diurna entre las 12:00 y las 13:45 aproximadamente. 2. Analizar Analizar los datos datos comprendi comprendidos dos entre entre [1570 : 1720] que 1720] que corresponden a una franja horaria de dos horas y media aproximadamente al final del periodo nocturno. nocturno. 3. Observando Observando la transformada wav wavelet elet trate de identificar donde donde se produjo un fallo en los sensores de medida. Sugerencia: Examine directamente los detalles D1 , D 2 y D 3 despreciando el resto. 4. Reduzca el ruido. Puede hacerlo hacerlo directamente directamente o utilizando la herramienta herramienta gráfica correspondiente. Sugerencia: Examine las aproximaciones de dos días sucesivos, sucesivos, el primero primero sin fallo en los aparatos aparatos de medida y el segundo con fallo. 5. Identifiqu Identifiquee patrones en los detalles, por ejemplo, aquellos aquellos producidos por temporizadores en los sistemas de agua caliente. Sugerencia: Examine detalles D2 , D 3 y D 4 alrededor de los puntos 1350, 1383 y 1415. 6. Identifique Identifique datos (aproximadament (aproximadamentee media hora) que han sido sustituidos con estimaciones (calculadas con splines) en lugar de los valores reales correspondientes. Sugerencia: Examine los datos alrededor del 2870 y utilice las pequeñas variaciones que aparecen en D1 para detectar los datos sustituidos.

C.4 C.4

Dete Detecc cció ión n de auto auto-s -sem emej ejan anza za (est (estru ruct ctur ura a frac frac-tal)

Veamos ahora una aplicación de transformada wavelet continua. El objetivo es mostrar que el análisis con wavelets puede detectar auto-semejanza (estruc-

C.4 Detección de auto-semejanza (estructura fractal)

135

tura fractal) en señales. Desde un punto de vista intuitivo la transformada wavelet continua es como un “índice de parecido” entre la señal y la wavelet. Cuanto mayor es el índice mayor es el parecido y viceversa. Los índices son los coeficientes de transformada. Si una señal es auto-similar a diferentes escalas, entonces los coeficientes (“índ (“índice ice de parecid parecido”) o”) tambié también n serán serán simila similares res a difere diferent ntes es escalas escalas.. Cuando Cuando representamos gráficamente el escalograma de la señal con las escalas en el eje vertical vertical esta auto-semejanza auto-semejanza en la señal genera un “patrón “patrón característico” característico” de las estructuras fractales auto-semejantes.

Ejercicio C.2 En el fichero vonkoch.mat puede encontrarse una señal construida de una manera recursiva. 1. Analizar Analizar la señal con la transform transformada ada wavelet wavelet continu continuaa del tipo coif3. Los niveles de escalas serán desde 2 a 128 con incrementos de 2. Una vez construido el escalograma es útil hacer uso de la herramienta de zoom. 2. Analiza Analizarr la señal señal con otros tipos de wav wavelets elets y observ observar ar si también también se genera en el escalograma patrones propios de estructuras fractales auto similares. 3. Si construimos construimos el escalograma escalograma con un conjunto conjunto de escalas escalas disjunto disjunto con el anterior ¿Qué podemos observar?

Bibliografía [1] E. Aboufadel Aboufadel and Schlicker Schlicker S., Discovering Wavelets , Jonh Wiley & Sons, 1999. [2] M. Frazier Frazier,, An Introduction to Wavelets through Linear Algebra , Springer, 1999. [3] S. Mallat, A wavelet Tour of Signal Processing , Academic Press, 1999. [4] The MathW MathWorks, Getting Started with MATLAB . [5] M. Misiti, Misiti, Y. Misiti, G. Oppenheim, Oppenheim, and J.M. Poggi, Poggi, Wavelet Toolbox for use with MATLAB , The MathWorks. [6] M. Misiti, Y. Misiti, G. Oppenheim, and J.M. Poggi, Décompositio en

ondelettes et méthodes comparatives: Étudy d’une courbe de charge électrique , Revue de Statistique Appliqué 42 (1994), no. 2, 57–77. [7] J.S. Walker, alker, A Primer Primer on Wavelet aveletss and their their Scient Scientific ific Applic Applicati ations ons , Chapman & Hall/CRC, 1999. [8] D.F. Waln Walnut, ut, An Introduction to Wavelet Analysis , Birkhäuser, 2002.

137

Tratamiento de señales digitales mediante wavelets y su uso con MATLAB - Félix Martínez Giménez

Recommend Documents