¿cuanta gente a visto el blog?

Piense en el procesamiento de señales digitales DSP en Python, versión 1.1.4

 Think DSP

Procesamiento Digital de Señales con Python Versión 1.1.4 Allen B. Downey Green Tea Press Needham, Massachusetts

Explicación: Este texto se refiere a un libro titulado "Think DSP" (Piensa en DSP, por sus siglas en inglés), que aborda el tema del Procesamiento Digital de Señales utilizando el lenguaje de programación Python. El autor del libro es Allen B. Downey. La editorial es Green Tea Press y se encuentra ubicada en Needham, Massachusetts.

El procesamiento digital de señales es una disciplina que se ocupa de analizar, modificar y sintetizar señales digitales para extraer información o realizar tareas específicas. Este libro en particular se centra en cómo aplicar técnicas de procesamiento de señales utilizando el lenguaje de programación Python, que es ampliamente utilizado en el campo de la ingeniería.

El libro podría ser útil para estudiantes de sexto semestre de ingeniería, ya que les proporcionaría una introducción al procesamiento digital de señales y les mostraría cómo utilizar Python como herramienta para implementar y experimentar con diferentes técnicas.

Prefacio El procesamiento de señales es uno de mis temas favoritos. Es útil en muchas áreas de la ciencia y la ingeniería, y si comprendes las ideas fundamentales, proporciona perspectivas sobre muchas cosas que vemos en el mundo, especialmente las cosas que oímos.

Pero a menos que hayas estudiado ingeniería eléctrica o mecánica, es probable que no hayas tenido la oportunidad de aprender sobre el procesamiento de señales. El problema es que la mayoría de los libros (y las clases que los utilizan) presentan el material de abajo hacia arriba, comenzando con abstracciones matemáticas como los fasores. Y tienden a ser teóricos, con pocas aplicaciones y poca relevancia aparente.

La premisa de este libro es que si sabes programar, puedes usar esa habilidad para aprender otras cosas y divertirte haciéndolo.

Con un enfoque basado en la programación, puedo presentar las ideas más importantes de inmediato. Al final del primer capítulo, podrás analizar grabaciones de sonido y otras señales, y generar nuevos sonidos. Cada capítulo presenta una nueva técnica y una aplicación que puedes aplicar a señales reales. En cada paso, aprendes cómo usar una técnica primero y luego cómo funciona.

Este enfoque es más práctico y, espero que estés de acuerdo, más divertido.

0.1 ¿Para quién es este libro? Los ejemplos y el código de soporte de este libro están en Python. Debes conocer el núcleo de Python y estar familiarizado con las características orientadas a objetos, al menos utilizando objetos, aunque no necesariamente definiendo los tuyos propios.

Si aún no estás familiarizado con Python, es posible que desees comenzar con mi otro libro, "Think Python", que es una introducción a Python para personas que nunca han programado, o "Learning Python" de Mark Lutz, que podría ser mejor para personas con experiencia en programación.

0.2 Usando el código El código y las muestras de sonido utilizadas en este libro están disponibles en https://github.com/AllenDowney/ThinkDSP. Git es un sistema de control de versiones que te permite realizar un seguimiento de los archivos que conforman un proyecto. Una colección de archivos bajo el control de Git se llama "repositorio". GitHub es un servicio de alojamiento que proporciona almacenamiento para repositorios de Git y una interfaz web conveniente.

La página principal de GitHub para mi repositorio ofrece varias formas de trabajar con el código:

  • Puedes crear una copia de mi repositorio en GitHub presionando el botón "Fork". Si aún no tienes una cuenta en GitHub, deberás crear una. Después de hacer el fork, tendrás tu propio repositorio en GitHub que podrás usar para realizar un seguimiento del código que escribas mientras trabajas en este libro. Luego puedes clonar el repositorio, lo que significa que copias los archivos en tu computadora.

  • O puedes clonar mi repositorio. No necesitas una cuenta en GitHub para hacer esto, pero no podrás enviar tus cambios de vuelta a GitHub.

  • Si no deseas usar Git en absoluto, puedes descargar los archivos en un archivo Zip utilizando el botón en la esquina inferior derecha de la página de GitHub.

Todo el código está escrito para funcionar tanto en Python 2 como en Python 3 sin necesidad de traducción.

Desarrollé este libro utilizando Anaconda de Continuum Analytics, que es una distribución gratuita de Python que incluye todos los paquetes que necesitarás para ejecutar el código (y muchos más). Encontré que Anaconda es fácil de instalar. Por defecto, realiza una instalación a nivel de usuario, no a nivel del sistema, por lo que no necesitas privilegios administrativos. Y admite tanto Python 2 como Python 3. Puedes descargar Anaconda desde https://www.anaconda.com/distribution/.

Si no deseas utilizar Anaconda, necesitarás los siguientes paquetes:

• NumPy para cálculos numéricos básicos, http://www.numpy.org/; • SciPy para cálculos científicos, http://www.scipy.org/; • matplotlib para visualización, http://matplotlib.org/.

Aunque estos paquetes son comúnmente utilizados, no están incluidos en todas las instalaciones de Python y pueden ser difíciles de instalar en algunos entornos. Si tienes problemas para instalarlos, te recomiendo utilizar Anaconda o alguna otra distribución de Python que incluya estos paquetes.

La mayoría de los ejercicios utilizan scripts de Python, pero algunos también utilizan notebooks de Jupyter. Si no has utilizado Jupyter antes, puedes obtener más información en http://jupyter.org.

Hay tres formas en las que puedes trabajar con los notebooks de Jupyter:

  1. Ejecutar Jupyter en tu computadora: Si has instalado Anaconda, es probable que ya tengas Jupyter por defecto. Para verificarlo, inicia el servidor desde la línea de comandos de esta manera: $ jupyter notebook

    Si no está instalado, puedes instalarlo en Anaconda de la siguiente manera: $ conda install jupyter

    Al iniciar el servidor, debería abrir tu navegador web por defecto o crear una nueva pestaña en una ventana de navegador abierta.

  2. Ejecutar Jupyter en Binder: Binder es un servicio que ejecuta Jupyter en una máquina virtual. Si sigues este enlace, http://mybinder.org/repo/AllenDowney/ThinkDSP, deberías obtener una página de inicio de Jupyter con los notebooks de este libro y los datos y scripts de soporte. Puedes ejecutar los scripts y modificarlos para ejecutar tu propio código, pero la máquina virtual en la que se ejecuta es temporal. Cualquier cambio que realices desaparecerá, junto con la máquina virtual, si la dejas inactiva durante más de una hora.

  3. Ver notebooks en nbviewer: Cuando nos refiramos a los notebooks más adelante en el libro, proporcionaremos enlaces a nbviewer, que ofrece una vista estática del código y los resultados. Puedes utilizar estos enlaces para leer los notebooks y escuchar los ejemplos, pero no podrás modificar o ejecutar el código, ni utilizar los widgets interactivos.

¡Buena suerte y diviértete!

Lista de colaboradores Si tienes alguna sugerencia o corrección, por favor envía un correo electrónico a downey@allendowney.com. Si realizo un cambio basado en tus comentarios, te agregaré a la lista de colaboradores (a menos que solicites ser omitido). Si incluyes al menos parte de la oración en la que se encuentra el error, me facilitará la búsqueda. Los números de página y sección también son útiles, pero no tan fáciles de manejar. ¡Gracias!

  • Antes de comenzar a escribir, mis ideas sobre este libro se beneficiaron de conversaciones con Boulos Harb en Google y Aurelio Ramos, anteriormente en Harmonix Music Systems.
  • Durante el semestre de otoño de 2013, Nathan Lintz e Ian Daniher trabajaron conmigo en un proyecto de estudio independiente y me ayudaron con el primer borrador de este libro.
  • En el foro de DSP de Reddit, el usuario anónimo RamjetSoundwave me ayudó a solucionar un problema con mi implementación de Ruido Browniano. Y andodli encontró un error tipográfico.
  • En la primavera de 2015, tuve el placer de enseñar este material junto con el profesor Oscar Mur-Miranda y el profesor Siddhartan Govindasamy. Ambos hicieron muchas sugerencias y correcciones.
  • Silas Gyger corrigió un error aritmético.
  • Giuseppe Masetti envió varias sugerencias muy útiles.
  • Eric Peters envió muchas sugerencias útiles. Un agradecimiento especial a Freesound, que es la fuente de muchas de las muestras de sonido que utilizo en este libro, y a los usuarios de Freesound que subieron esos sonidos. Incluyo algunos de sus archivos de onda en el repositorio de GitHub de este libro, utilizando los nombres de archivo originales, por lo que debería ser fácil encontrar sus fuentes. Desafortunadamente, la mayoría de los usuarios de Freesound no hacen públicos sus nombres reales, por lo que solo puedo agradecerles utilizando sus nombres de usuario. Las muestras utilizadas en este libro fueron contribuidas por los usuarios de Freesound: iluppai, wcfl10, thirsk, docquesting, kleeb, landup, zippi1, themusicalnomad, bcjordan, rockwehrmann, marcgascon7, jcveliz. ¡Gracias a todos ustedes!

Contenidos

Prefacio v 0.1 ¿Para quién es este libro? . . . . . . . . . . . . . . . . . . . . . v 0.2 Uso del código . . . . . . . . . . . . . . . . . . . . . . . . . . . vi 1 Sonidos y señales 1 1.1 Señales periódicas . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Descomposición espectral . . . . . . . . . . . . . . . . . . . . 3 1.3 Señales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.4 Lectura y escritura de ondas . . . . . . . . . . . . . . . . . . 7 1.5 Espectros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.6 Objetos de onda . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.7 Objetos de señal . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.8 Ejercicios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2 Armónicos 13 2.1 Ondas triangulares . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2 Ondas cuadradas . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3 Alias . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.4 Cálculo del espectro . . . . . . . . . . . . . . . . . . . . . . . 19 2.5 Ejercicios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

3 Señales no periódicas 23 3.1 Chirrido lineal . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.2 Chirrido exponencial . . . . . . . . . . . . . . . . . . . . . . . 25 3.3 Espectro de un chirrido . . . . . . . . . . . . . . . . . . . . . 26 3.4 Espectrograma . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.5 El límite de Gabor . . . . . . . . . . . . . . . . . . . . . . . . 28 3.6 Fugas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.7 Ventaneo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.8 Implementación de espectrogramas . . . . . . . . . . . . . . 32 3.9 Ejercicios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4 Ruido 37 4.1 Ruido no correlacionado . . . . . . . . . . . . . . . . . . . . 37 4.2 Espectro integrado . . . . . . . . . . . . . . . . . . . . . . . . 40 4.3 Ruido Browniano . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.4 Ruido Rosa . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 4.5 Ruido Gaussiano . . . . . . . . . . . . . . . . . . . . . . . . . 46 4.6 Ejercicios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 5 Autocorrelación 51 5.1 Correlación . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 5.2 Correlación en serie . . . . . . . . . . . . . . . . . . . . . . . 54 5.3 Autocorrelación . . . . . . . . . . . . . . . . . . . . . . . . . 55 5.4 Autocorrelación de señales periódicas . . . . . . . . . . . . . 56 5.5 Correlación como producto escalar . . . . . . . . . . . . . . 60 5.6 Uso de NumPy . . . . . . . . . . . . . . . . . . . . . . . . . . 60 5.7 Ejercicios . . .

Capítulo 1: Sonidos y señales

Una señal representa una cantidad que varía en el tiempo. Esta definición es bastante abstracta, así que comencemos con un ejemplo concreto: el sonido. El sonido es una variación en la presión del aire. Una señal de sonido representa las variaciones en la presión del aire a lo largo del tiempo.

Un micrófono es un dispositivo que mide estas variaciones y genera una señal eléctrica que representa el sonido. Un altavoz es un dispositivo que toma una señal eléctrica y produce sonido. Los micrófonos y los altavoces se llaman transductores porque transducen, o convierten, señales de una forma a otra.

Este libro trata sobre el procesamiento de señales, que incluye procesos para sintetizar, transformar y analizar señales. Me centraré en las señales de sonido, pero los mismos métodos se aplican a señales electrónicas, vibraciones mecánicas y señales en muchos otros dominios.

También se aplican a señales que varían en el espacio en lugar del tiempo, como la elevación a lo largo de un sendero de senderismo. Y se aplican a señales en más de una dimensión, como una imagen, que se puede pensar como una señal que varía en un espacio bidimensional. O una película, que es una señal que varía en un espacio bidimensional y en el tiempo.

Pero comenzaremos con el sonido simple unidimensional.

El código para este capítulo se encuentra en chap01.ipynb, que está en el repositorio de este libro (ver Sección 0.2). También puedes verlo en http://tinyurl.com/thinkdsp01.


Figura 1.1: Segmento de una grabación de una campana.

1.1 Señales periódicas

Comenzaremos con señales periódicas, que son señales que se repiten a sí mismas después de cierto período de tiempo. Por ejemplo, si golpeas una campana, esta vibra y genera sonido. Si grabas ese sonido y graficas la señal transducida, se vería como la Figura 1.1.

Esta señal se asemeja a una sinusoidal, lo que significa que tiene la misma forma que la función seno trigonométrica.

Puedes ver que esta señal es periódica. Elegí una duración que muestra tres repeticiones completas, también conocidas como ciclos. La duración de cada ciclo, llamada período, es de aproximadamente 2,3 ms.

La frecuencia de una señal es el número de ciclos por segundo, que es el inverso del período. Las unidades de frecuencia son ciclos por segundo, o Hertzios, abreviado "Hz" (Estrictamente hablando, el número de ciclos es una cantidad adimensional, por lo que un Hertzio es realmente "por segundo").

La frecuencia de esta señal es de aproximadamente 439 Hz, ligeramente más baja que los 440 Hz, que es el tono de afinación estándar para la música orquestal. El nombre musical de esta nota es La, o más específicamente, La4. Si no estás familiarizado con la "notación científica del tono", el sufijo numérico indica en qué octava se encuentra la nota. La4 es el La por encima del do central. La5 es una octava más alta. Ver http://en.wikipedia.org/wiki/Scientific_pitch_notation.

Un diapasón genera una sinusoidal porque la vibración de las púas es una forma de movimiento armónico simple. La mayoría de los instrumentos musicales producen señales periódicas, pero la forma de estas señales no es sinusoidal. Por ejemplo, la Figura 1.2 muestra...

Figura 1.2: Segmento de una grabación de un violín.

un segmento de una grabación de un violín interpretando el tercer movimiento del Quinteto de Cuerdas No. 5 en Mi de Boccherini.

Una vez más, podemos ver que la señal es periódica, pero la forma de la señal es más compleja. La forma de una señal periódica se llama forma de onda. La mayoría de los instrumentos musicales producen formas de onda más complejas que una sinusoidal. La forma de la forma de onda determina el timbre musical, que es nuestra percepción de la calidad del sonido. Por lo general, percibimos las formas de onda complejas como ricas, cálidas y más interesantes que las sinusoidales.

1.2 Descomposición espectral

El tema más importante en este libro es la descomposición espectral, que es la idea de que cualquier señal se puede expresar como la suma de sinusoides con diferentes frecuencias.

La idea matemática más importante en este libro es la transformada discreta de Fourier, o DFT, que toma una señal y produce su espectro. El espectro es el conjunto de sinusoides que se suman para producir la señal.

Y el algoritmo más importante en este libro es la transformada rápida de Fourier, o FFT, que es una forma eficiente de calcular la DFT.

Por ejemplo, la Figura 1.3 muestra el espectro de la grabación del violín en la Figura 1.2. El eje x es el rango de frecuencias que componen la señal. El eje y muestra la intensidad o amplitud de cada componente de frecuencia.

Figura 1.3: Espectro de un segmento de la grabación del violín.

El componente de frecuencia más bajo se llama frecuencia fundamental. La frecuencia fundamental de esta señal está cerca de 440 Hz (en realidad un poco más baja o "bemol").

En esta señal, la frecuencia fundamental tiene la amplitud más grande, por lo que también es la frecuencia dominante. Normalmente, el tono percibido de un sonido está determinado por la frecuencia fundamental, incluso si no es dominante.

Las otras puntas en el espectro se encuentran en las frecuencias 880, 1320, 1760 y 2200, que son múltiplos enteros de la fundamental. Estos componentes se llaman armónicos porque son armónicos musicalmente con la fundamental:

  • 880 es la frecuencia de A5, una octava por encima de la fundamental. Una octava es un doble en frecuencia.
  • 1320 es aproximadamente E6, que es una quinta justa por encima de A5. Si no estás familiarizado con intervalos musicales como "quinta justa", consulta https://es.wikipedia.org/wiki/Intervalo_(música).
  • 1760 es A6, dos octavas por encima de la fundamental.
  • 2200 es aproximadamente C#7, que es un tercer mayor por encima de A6.

Estos armónicos componen las notas de un acorde de La mayor, aunque no todas en la misma octava. Algunas de ellas son solo aproximadas debido a que las notas que componen la música occidental se han ajustado para la igualdad de temperamento (ver http://es.wikipedia.org/wiki/Temperamento_igual).

Dado los armónicos y sus amplitudes, puedes reconstruir la señal sumando sinusoides. A continuación, veremos cómo.

1.3 Señales

He creado un módulo de Python llamado thinkdsp.py que contiene clases y funciones para trabajar con señales y espectros. Lo encontrarás en el repositorio de este libro (ver Sección 0.2).

Para representar señales, thinkdsp proporciona una clase llamada Signal, que es la clase principal para varios tipos de señales, incluyendo Sinusoid, que representa tanto señales seno como coseno.

thinkdsp proporciona funciones para crear señales seno y coseno: cos_sig = thinkdsp.CosSignal(freq=440, amp=1.0, offset=0) sin_sig = thinkdsp.SinSignal(freq=880, amp=0.5, offset=0)

freq es la frecuencia en Hz. amp es la amplitud en unidades no especificadas, donde 1.0 se define como la amplitud más grande que podemos grabar o reproducir.

offset es un desplazamiento de fase en radianes. El desplazamiento de fase determina dónde comienza la señal en el período. Por ejemplo, una señal seno con offset=0 comienza en sin 0, que es 0. Con offset=pi/2 comienza en sin π/2, que es 1.

Las señales tienen un método add, por lo que puedes usar el operador + para sumarlas: mix = sin_sig + cos_sig

El resultado es una SumSignal, que representa la suma de dos o más señales.

Una Signal es básicamente una representación en Python de una función matemática. La mayoría de las señales se definen para todos los valores de t, desde menos infinito hasta infinito.

No puedes hacer mucho con una Signal hasta que la evalúes. En este contexto, "evaluar" significa tomar una secuencia de puntos en el tiempo, ts, y calcular los valores correspondientes de la señal, ys. Represento ts e ys utilizando matrices de NumPy y los encapsulo en un objeto llamado Wave.

Un Wave representa una señal evaluada en una secuencia de puntos en el tiempo. Cada punto en el tiempo se llama frame (un término tomado de películas y videos). La propia medición se llama muestra, aunque a veces se usan "frame" y "sample" indistintamente.

Signal proporciona make_wave, que devuelve un nuevo objeto Wave: wave = mix.make_wave(duration=0.5, start=0, framerate=11025)


figura 1.4 Segmento de una mezcla de dos señales sinusoidales.

La duración es la longitud de la Wave en segundos. Start es el tiempo de inicio, también en segundos. La velocidad de fotogramas es el número (entero) de cuadros por segundo, que también es el número de muestras por segundo.

11,025 cuadros por segundo es una de las varias velocidades de fotogramas comúnmente utilizadas en formatos de archivos de audio, incluidos los archivos de audio de forma de onda (WAV) y mp3.

Este ejemplo evalúa la señal desde t=0 hasta t=0.5 en 5,513 cuadros igualmente espaciados (porque 5,513 es la mitad de 11,025). El tiempo entre cuadros, o paso de tiempo, es de 1/11025 segundos, aproximadamente 91 µs.

Wave proporciona un método plot que utiliza pyplot. Puedes trazar la onda de esta manera: wave.plot() pyplot.show() pyplot es parte de matplotlib; está incluido en muchas distribuciones de Python, o es posible que tengas que instalarlo.

En freq=440 hay 220 períodos en 0.5 segundos, por lo que esta trama se vería como un bloque sólido de color. Para acercarse a un número pequeño de períodos, podemos usar segment, que copia un segmento de una Wave y devuelve una nueva onda: period = mix.period segment = wave.segment(start=0, duration=period*3) period es una propiedad de una Signal; devuelve el período en segundos.

start y duration están en segundos. Este ejemplo copia los primeros tres períodos de la mezcla. El resultado es un objeto Wave.

Si trazamos el segmento, se verá como la Figura 1.4. Esta señal contiene dos componentes de frecuencia, por lo que es más complicada que la señal del diapasón, pero menos complicada que el violín.

1.4 Lectura y escritura de Waves

thinkdsp proporciona read_wave, que lee un archivo WAV y devuelve una Wave: violin_wave = thinkdsp.read_wave('input.wav')

Y Wave proporciona write, que escribe un archivo WAV: wave.write(filename='output.wav')

Puedes escuchar la Wave con cualquier reproductor multimedia que reproduzca archivos WAV. En sistemas UNIX, utilizo aplay, que es simple, robusto y está incluido en muchas distribuciones de Linux.

thinkdsp también proporciona play_wave, que ejecuta el reproductor multimedia como un subprocess: thinkdsp.play_wave(filename='output.wav', player='aplay')

Utiliza aplay de forma predeterminada, pero puedes proporcionar el nombre de otro reproductor.

1.5 Espectros

Wave proporciona make_spectrum, que devuelve un Spectrum: spectrum = wave.make_spectrum()

Y Spectrum proporciona plot: spectrum.plot()

Spectrum proporciona tres métodos que modifican el espectro:

  • low_pass aplica un filtro paso bajo, lo que significa que los componentes por encima de una determinada frecuencia de corte se atenúan (es decir, se reducen en magnitud) por un factor.
  • high_pass aplica un filtro paso alto, lo que significa que atenúa los componentes por debajo de la frecuencia de corte.
  • band_stop atenúa los componentes en la banda de frecuencias entre dos frecuencias de corte.

En este ejemplo, se atenúan todas las frecuencias por encima de 600 en un 99%: spectrum.low_pass(cutoff=600, factor=0.01)

1.5 Relaciones entre las clases en thinkdsp.

spectrum.low_pass(cutoff=600, factor=0.01) Un filtro paso bajo elimina los sonidos brillantes de alta frecuencia, por lo que el resultado suena amortiguado y más oscuro. Para escuchar cómo suena, puedes convertir el Spectrum nuevamente en una Wave y luego reproducirla. wave = spectrum.make_wave() wave.play('temp.wav')

El método play escribe la wave en un archivo y luego la reproduce. Si usas Jupyter notebooks, puedes usar make_audio, que crea un widget de audio que reproduce el sonido.

1.6 Objetos Wave

No hay nada muy complicado en thinkdsp.py. La mayoría de las funciones que proporciona son envoltorios simples alrededor de funciones de NumPy y SciPy.

Las clases principales en thinkdsp son Signal, Wave y Spectrum. Dada una Signal, puedes crear una Wave. Dada una Wave, puedes crear un Spectrum y viceversa. Estas relaciones se muestran en la Figura 1.5.

Un objeto Wave contiene tres atributos: ys es una matriz NumPy que contiene los valores de la señal; ts es una matriz de los tiempos en los que se evaluó o muestreó la señal; y framerate es el número de muestras por unidad de tiempo. La unidad de tiempo suele ser segundos, pero no tiene que serlo. En uno de mis ejemplos, es días.

Wave también proporciona tres propiedades de solo lectura: start, end y duration. Si modificas ts, estas propiedades cambian en consecuencia.

Para modificar una wave, puedes acceder a ts y ys directamente. Por ejemplo: wave.ys *= 2 wave.ts += 1

La primera línea escala la wave en un factor de 2, haciéndola más fuerte. La segunda línea desplaza la wave en el tiempo, haciendo que comience 1 segundo más tarde.

Pero Wave proporciona métodos que realizan muchas operaciones comunes. Por ejemplo, las mismas dos transformaciones podrían escribirse así:

1.7 Objetos Signal

Signal es una clase padre que proporciona funciones comunes a todos los tipos de señales, como make_wave. Las clases hijas heredan estos métodos y proporcionan evaluate, que evalúa la señal en una secuencia dada de tiempos.

Por ejemplo, Sinusoid es una clase hija de Signal, con la siguiente definición:

class Sinusoid(Signal): def init(self, freq=440, amp=1.0, offset=0, func=np.sin): Signal.init(self) self.freq = freq self.amp = amp self.offset = offset self.func = func

Los parámetros de init son:

  • freq: frecuencia en ciclos por segundo, o Hz.
  • amp: amplitud. Las unidades de amplitud son arbitrarias, generalmente elegidas de manera que 1.0 corresponda a la entrada máxima desde un micrófono o la salida máxima hacia un altavoz.
  • offset: indica dónde comienza la señal en su período; el desplazamiento se mide en radianes, por razones que explicaré a continuación.
  • func: una función de Python utilizada para evaluar la señal en un punto particular en el tiempo. Por lo general, es np.sin o np.cos, lo que genera una señal sinusoidal o cosinusoidal.

Al igual que muchos métodos init, este simplemente almacena los parámetros para su uso futuro.

Signal proporciona make_wave, que se ve así:

def make_wave(self, duration=1, start=0, framerate=11025): n = round(duration * framerate) ts = start + np.arange(n) / framerate ys = self.evaluate(ts) return Wave(ys, ts, framerate=framerate)

En este método, start y duration son el tiempo de inicio y la duración en segundos. El framerate es el número de fotogramas (muestras) por segundo.

n es el número de muestras, y ts es una matriz NumPy de tiempos de muestra.

Para calcular los ys, make_wave invoca a evaluate, que es proporcionada por Sinusoid:

def evaluate(self, ts): phases = PI2 * self.freq * ts + self.offset ys = self.amp * self.func(phases) return ys

Vamos a desglosar esta función paso a paso:

  1. self.freq es la frecuencia en ciclos por segundo, y cada elemento de ts es un tiempo en segundos, por lo que su producto es el número de ciclos desde el tiempo de inicio.

  2. PI2 es una constante que almacena 2π. Multiplicar por PI2 convierte de ciclos a fase. Puedes pensar en la fase como "ciclos desde el tiempo de inicio" expresados en radianes. Cada ciclo es de 2π radianes.

  3. self.offset es la fase cuando t es ts[0]. Tiene el efecto de desplazar la señal hacia la izquierda o hacia la derecha en el tiempo.

  4. Si self.func es np.sin o np.cos, el resultado es un valor entre -1 y +1.

  5. Multiplicar por self.amp produce una señal que oscila entre -self.amp y +self.amp.

En notación matemática, evaluate se escribe así:

y = A cos(2πf t + φ0)

donde A es la amplitud, f es la frecuencia, t es el tiempo y φ0 es el desplazamiento de fase. Puede parecer que he escrito mucho código para evaluar una expresión simple, pero como veremos, este código proporciona un marco de trabajo para tratar con todo tipo de señales, no solo sinusoides.

1.8 ejercicios

Ejercicio 1.1: Si tienes Jupyter Notebook, carga el archivo chap01.ipynb, léelo y ejecuta los ejemplos. También puedes ver este cuaderno en http://tinyurl.com/thinkdsp01.

Para resolver este ejercicio, sigue estos pasos:

  1. Descarga el código para este libro siguiendo las instrucciones en la Sección 0.2.
  2. Abre Jupyter Notebook.
  3. Navega hasta la ubicación donde descargaste el código del libro.
  4. Busca el archivo "chap01.ipynb" y ábrelo.
  5. Lee cuidadosamente el contenido del cuaderno y ejecuta los ejemplos de código proporcionados.
  6. Si deseas, también puedes acceder al cuaderno en línea utilizando el enlace proporcionado.

Al realizar estos pasos, podrás familiarizarte con el contenido del capítulo 1 y ejecutar los ejemplos de código correspondientes. Esto te ayudará a comprender los conceptos presentados y te preparará para resolver los siguientes ejercicios.

Ejercicio 1.2:

Para resolver este ejercicio, sigue estos pasos:

Ve a http://freesound.org y descarga una muestra de sonido que incluya música, habla u otros sonidos con una frecuencia fundamental bien definida. ver ejemplos de soluciones a los ejercicios

selecciona un segmento de aproximadamente medio segundo en el que la frecuencia fundamental sea constante. Puedes utilizar software de edición de audio para recortar y guardar este segmento por separado. Utiliza la función read_wave de thinkdsp para leer el archivo de audio descargado y crear un objeto Wave. Utiliza el método make_spectrum() del objeto Wave para obtener el espectro de frecuencia de la señal. Utiliza el método plot() del objeto Spectrum para graficar el espectro. Observa la estructura armónica en el espectro y cómo se relaciona con el timbre del sonido. Los armónicos son múltiplos enteros de la frecuencia fundamental y su amplitud relativa determina el timbre o calidad tonal del sonido. Utiliza los métodos high_pass(), low_pass() y band_stop() del objeto Spectrum para filtrar algunos de los armónicos. Convierte el espectro filtrado de nuevo a una señal Wave utilizando el método make_wave() del objeto Spectrum. Escucha la señal Wave resultante utilizando el método play() del objeto Wave. Observa cómo se relaciona el sonido resultante con los cambios que realizaste en el espectro. Al filtrar armónicos específicos, estás alterando la composición tonal y el timbre del sonido. Recuerda que puedes consultar el archivo chap01soln.ipynb del código fuente descargado par

Ejercicio 1.3: Sintetizar una señal compuesta creando objetos SinSignal y CosSignal ysumándolos. Evalúe la señal para obtener una onda y escúchela. Calcule su espectro y grafítelo. ¿Qué sucede si agrega componentes de frecuencia que no son múltiplos de la fundamental?

Para resolver este ejercicio, sigue estos pasos:

  1. Importa los módulos necesarios:

    python
    import numpy as np import thinkdsp from thinkdsp import SinSignal, CosSignal
  2. Crea objetos SinSignal y CosSignal con diferentes frecuencias y amplitudes:

    python
    frequency1 = 440 # Frecuencia fundamental en Hz frequency2 = 660 # Frecuencia no múltiplo de la fundamental en Hz amplitude1 = 1.0 # Amplitud de la señal sinusoidal amplitude2 = 0.5 # Amplitud de la señal cosinusoidal sin_signal = SinSignal(freq=frequency1, amp=amplitude1) cos_signal = CosSignal(freq=frequency2, amp=amplitude2)
  3. Combina las señales sumándolas:

    python
    compound_signal = sin_signal + cos_signal
  4. Evalúa la señal para obtener una forma de onda (Wave):

    python
    duration = 2.0 # Duración de la forma de onda en segundos framerate = 44100 # Tasa de muestreo en fotogramas por segundo wave = compound_signal.make_wave(duration=duration, start=0, framerate=framerate)
  5. Escucha la forma de onda utilizando el método play() del objeto Wave:

    python
    wave.play()
  6. Calcula el espectro de la forma de onda utilizando el método make_spectrum() y luego grafícalo utilizando el método plot() del objeto Spectrum:

    python
    spectrum = wave.make_spectrum() spectrum.plot()
  7. Observa el espectro y nota cómo los componentes de frecuencia que no son múltiplos de la frecuencia fundamental generan picos adicionales en el espectro. Estos componentes de frecuencia no armónicos introducen contenido espectral adicional y pueden afectar la calidad tonal y el timbre del sonido.

Recuerda que puedes ajustar los valores de frecuencia, amplitud, duración y tasa de muestreo según tus preferencias.

Ejercicio 1.4 Escribe una función llamada estiramiento que tome una Onda y un estiramiento factor y acelera o ralentiza la onda modificando ts y framerate. Sugerencia: solo debe tomar dos líneas de código.

Para resolver este ejercicio, puedes seguir estos pasos:

  1. Define una función llamada stretch que tome como argumentos un objeto Wave y un factor de estiramiento:

    python
    def stretch(wave, factor):
  2. Modifica los atributos ts y framerate del objeto Wave según el factor de estiramiento:

    python
    wave.ts /= factor wave.framerate *= factor

El código completo de la función stretch sería el siguiente:

python
def stretch(wave, factor): wave.ts /= factor wave.framerate *= factor

Con esta función, puedes estirar o comprimir una forma de onda modificando los valores de tiempo y la tasa de muestreo. Recuerda que el factor de estiramiento debe ser mayor que 1 para acelerar la forma de onda y menor que 1 para desacelerarla.

Por ejemplo, si tienes una forma de onda llamada wave y deseas estirarla por un factor de 2 para duplicar su duración, puedes llamar a la función de la siguiente manera:

python
stretch(wave, 2)

Esto modificará la forma de onda wave para que tenga una duración duplicada y se reproduzca a una velocidad más lenta.

El Capítulo 1 del libro "Think DSP" cubre los conceptos básicos de señales y espectros. Algunos puntos importantes de este capítulo incluyen:

  1. Introducción a las señales: Se presenta el concepto de señales y se discuten las propiedades comunes de las señales, como la amplitud, frecuencia, periodo y fase.

  2. Representación de señales en Python: Se introduce el módulo thinkdsp.py, que proporciona clases y funciones para trabajar con señales y espectros en Python. Se exploran las clases Signal, Wave y Spectrum, que se utilizan para representar y manipular señales.

  3. Generación de señales: Se muestran ejemplos de cómo crear señales sinusoidales utilizando las clases SinSignal y CosSignal del módulo thinkdsp.py. Se explica cómo ajustar la frecuencia, amplitud y desplazamiento de fase de estas señales.

  4. Creación y manipulación de ondas: Se introduce la clase Wave, que representa una señal evaluada en una secuencia de puntos en el tiempo. Se explican los conceptos de duración, tiempo de inicio y frecuencia de muestreo. Se muestra cómo crear y manipular ondas utilizando el método make_wave y cómo visualizar y reproducir las ondas.

  5. Lectura y escritura de archivos WAV: Se muestra cómo leer y escribir archivos WAV utilizando las funciones read_wave y write del módulo thinkdsp.py. Se explica cómo utilizar reproductores de medios externos para escuchar las ondas.

  6. Espectros: Se introduce el concepto de espectro, que representa la descomposición de una señal en sus componentes de frecuencia. Se muestra cómo crear y visualizar espectros utilizando el método make_spectrum y el método plot de la clase Spectrum.

  7. Filtrado de espectros: Se presentan las funciones de filtrado disponibles en la clase Spectrum, como low_pass, high_pass y band_stop. Se muestra cómo aplicar estos filtros para atenuar o eliminar componentes de frecuencia específicos en un espectro.

  8. Relación entre el espectro y el timbre: Se realiza una conexión entre la estructura armónica de un espectro y el timbre de un sonido. Se muestra cómo los armónicos presentes en el espectro contribuyen a la percepción del timbre de un sonido.

En resumen, este capítulo proporciona una introducción sólida a los conceptos básicos de señales y espectros, y muestra cómo aplicar estos conceptos utilizando Python y el módulo thinkdsp.py. Se presentan ejemplos prácticos y se realizan ejercicios para que los lectores practiquen y refuercen los conceptos aprendidos.

Capítulo 2: Armónicos

En este capítulo presento varias nuevas formas de onda; examinaremos sus espectros para comprender su estructura armónica, que es el conjunto de sinusoides del que están compuestas.

También introduciré uno de los fenómenos más importantes en el procesamiento de señales digitales: el aliasing. Y explicaré un poco más sobre cómo funciona la clase Spectrum.

El código para este capítulo se encuentra en chap02.ipynb, que está en el repositorio de este libro (ver Sección 0.2). También puedes verlo en http://tinyurl.com/thinkdsp02.

2.1 Ondas triangulares

Un sinusoidal contiene solo una componente de frecuencia, por lo que su espectro tiene solo un pico. Formas de onda más complejas, como la grabación de un violín, producen DFT con muchos picos. En esta sección investigamos la relación entre las formas de onda y sus espectros.

Comenzaré con una forma de onda triangular, que es una versión lineal de un sinusoidal. La Figura 2.1 muestra una forma de onda triangular con una frecuencia de 200 Hz.

Para generar una onda triangular, puedes usar thinkdsp.TriangleSignal:

python
class TriangleSignal(Sinusoid): def evaluate(self, ts): cycles = self.freq * ts + self.offset / PI2
Figura 2.1: Segmento de una señal triangular a 200 Hz.
python
frac, _ = np.modf(cycles) 
ys = np.abs(frac - 0.5
ys = normalize(unbias(ys), self.amp) 
return ys

TriangleSignal hereda __ini__t de Sinusoid, por lo que toma los mismos argumentos: freq, amp y offset. La única diferencia está en evaluate. Como vimos antes, ts es la secuencia de tiempos de muestra en los que queremos evaluar la señal.

Hay muchas formas de generar una onda triangular. Los detalles no son importantes, pero así es como funciona evaluate:

  1. cycles es el número de ciclos desde el tiempo de inicio. np.modf divide el número de ciclos en la parte fraccionaria, que se almacena en frac, y la parte entera, que se ignora.
  2. frac es una secuencia que va desde 0 hasta 1 con la frecuencia dada. Restar 0.5 produce valores entre -0.5 y 0.5. Tomar el valor absoluto produce una forma de onda que va de zigzag entre 0.5 y 0.
  3. unbias desplaza la forma de onda hacia abajo para que esté centrada en 0; luego normalize la escala con la amplitud dada, amp.

Aquí está el código que genera la Figura 2.1:

python
signal = thinkdsp.TriangleSignal(200) signal.plot()
Figura 2.2: Espectro de una señal triangular a 200 Hz, mostrado en dos escalas verticales. La versión de la derecha corta el fundamental para mostrar los armónicos de manera más clara.

A continuación, podemos utilizar la Señal para crear una Onda y utilizar la Onda para crear un Espectro:

python
wave = signal.make_wave(duration=0.5, framerate=10000
spectrum = wave.make_spectrum() 
spectrum.plot()

La Figura 2.2 muestra dos vistas del resultado; la vista de la derecha está escalada para mostrar los armónicos de manera más clara. Como era de esperar, el pico más alto está en la frecuencia fundamental, 200 Hz, y hay picos adicionales en las frecuencias armónicas, que son múltiplos enteros de 200.

Pero una sorpresa es que no hay picos en los múltiplos pares: 400, 800, etc. Los armónicos de una onda triangular son todos múltiplos impares de la frecuencia fundamental, en este ejemplo 600, 1000, 1400, etc.

Otra característica de este espectro es la relación entre la amplitud y la frecuencia de los armónicos. Su amplitud disminuye en proporción al cuadrado de la frecuencia. Por ejemplo, la relación de frecuencia de los dos primeros armónicos (200 y 600 Hz) es 3, y la relación de amplitud es aproximadamente 9. La relación de frecuencia de los siguientes dos armónicos (600 y 1000 Hz) es 1.7, y la relación de amplitud es aproximadamente 1.7^2 = 2.9. Esta relación se llama estructura armónica.


Figura 2.3: Segmento de una señal cuadrada a 100 Hz.

2.2 Ondas cuadradas

thinkdsp también proporciona SquareSignal, que representa una señal cuadrada. Aquí está la definición de la clase:

python
class SquareSignal(Sinusoid):
    def evaluate(self, ts):
         cycles = self.freq * ts + self.offset / PI2
         frac, _ = np.modf(cycles)
         ys = self.amp * np.sign(unbias(frac))
        return ys

Al igual que TriangleSignal, SquareSignal hereda init de Sinusoid, por lo que toma los mismos parámetros.

Y el método evaluate es similar. Nuevamente, cycles es el número de ciclos desde el tiempo de inicio, y frac es la parte fraccional, que aumenta de 0 a 1 en cada período.

unbias desplaza frac para que aumente de -0.5 a 0.5, luego np.sign asigna los valores negativos a -1 y los valores positivos a 1. Multiplicar por amp produce una onda cuadrada que salta entre -amp y amp.

La Figura 2.3 muestra tres períodos de una onda cuadrada con frecuencia de 100 Hz, y la Figura 2.4 muestra su espectro.

Al igual que una onda triangular, la onda cuadrada contiene solo armónicos impares, por eso hay picos en 300, 500 y 700 Hz, etc. Pero la amplitud de los armónicos

Figura 2.4: Espectro de una señal cuadrada a 100 Hz.

La amplitud de los armónicos decae más lentamente que en la onda triangular. Específicamente, la amplitud disminuye en proporción a la frecuencia (no al cuadrado de la frecuencia).

Los ejercicios al final de este capítulo te brindan la oportunidad de explorar otras formas de onda y otras estructuras armónicas.

2.3 Aliasing Tengo una confesión. Elegí los ejemplos en la sección anterior cuidadosamente para evitar mostrarte algo confuso. Pero ahora es hora de confundirnos.

La Figura 2.5 muestra el espectro de una onda triangular a 1100 Hz, muestreada a 10,000 fotogramas por segundo. Nuevamente, la vista de la derecha está escalada para mostrar los armónicos.

Los armónicos de esta onda deberían estar en 3300, 5500, 7700 y 9900 Hz. En la figura, hay picos en 1100 y 3300 Hz, como se esperaba, pero el tercer pico está en 4500 Hz, no en 5500 Hz. El cuarto pico está en 2300 Hz, no en 7700 Hz. Y si te fijas bien, el pico que debería estar en 9900 Hz en realidad está en 100 Hz. ¿Qué está pasando?

El problema es que cuando evaluas la señal en puntos discretos en el tiempo, pierdes información sobre lo que sucedió entre las muestras. Para componentes de baja frecuencia, eso no es un problema, porque tienes muchas muestras por período. Pero si muestreas una señal a 5000 Hz con 10,000 fotogramas por segundo, solo tienes dos muestras por período. Resulta ser suficiente, pero si la frecuencia es más alta, no lo es.

Figura 2.5: Espectro de una señal triangular a 1100 Hz muestreada a 10,000 fotogramas por segundo. La vista de la derecha está escalada para mostrar los armónicos.

Para entender por qué esto sucede, generemos señales coseno a 4500 y 5500 Hz, y las muestreemos a 10,000 fotogramas por segundo: framerate = 10000 signal = thinkdsp.CosSignal(4500) duration = signal.period*5 segment = signal.make_wave(duration, framerate=framerate) segment.plot() signal = thinkdsp.CosSignal(5500) segment = signal.make_wave(duration, framerate=framerate) segment.plot()

La Figura 2.6 muestra el resultado. He trazado las señales con líneas grises delgadas y las muestras con líneas verticales, para facilitar la comparación de las dos ondas. El problema debería ser evidente: aunque las señales son diferentes, ¡las ondas son idénticas!

Cuando muestreamos una señal de 5500 Hz a 10,000 fotogramas por segundo, el resultado es indistinguible de una señal de 4500 Hz. Por la misma razón, una señal de 7700 Hz es indistinguible de 2300 Hz y una señal de 9900 Hz es indistinguible de 100 Hz.

Este efecto se llama aliasing porque cuando la señal de alta frecuencia se muestrea, parece ser una señal de baja frecuencia.

En este ejemplo, la frecuencia más alta que podemos medir es de 5000 Hz, que es la mitad de la frecuencia de muestreo. Las frecuencias por encima de 5000 Hz se pliegan por debajo de 5000 Hz,


Figura 2.6: Señales coseno a 4500 y 5500 Hz, muestreadas a 10,000 fotogramas por segundo. Las señales son diferentes, pero las muestras son idénticas.
razón por la cual este umbral a veces se denomina "frecuencia de plegado". A veces también se le llama la frecuencia de Nyquist. Ver http://en.wikipedia.org/wiki/Nyquist_frequency

Esto se debe a que cuando la frecuencia de aliasing cae por debajo de cero, el patrón de plegado continúa. Por ejemplo, el quinto armónico de la señal triangular a 1100 Hz se encuentra en 12,100 Hz. Plegado a 5000 Hz, aparecería en -2100 Hz, pero se pliega nuevamente a 0 Hz, volviendo a 2100 Hz. De hecho, se puede ver un pequeño pico en 2100 Hz en la Figura 2.4, y el siguiente en 4300 Hz.

2.4 Cálculo del espectro

Hemos visto el método make_spectrum de la clase Wave varias veces. Aquí está la implementación (sin incluir algunos detalles que veremos más adelante): from np.fft import rfft, rfftfreq

clase Wave:

def make_spectrum(self): n = len(self.ys)

d = 1 / self.framerate

 hs = rfft(self.ys) 

fs = rfftfreq(n, d) 

return Spectrum(hs, fs, self.framerate)  El parámetro self es un objeto Wave. n es el número de muestras en la onda, y d es el inverso de la velocidad de fotogramas, que es el tiempo entre muestras. np.fft es el módulo NumPy que proporciona funciones relacionadas con Fast Transformada de Fourier (FFT), que es un algoritmo eficiente que calcula la transformada discreta de Fourier (DFT).

El método make_spectrum utiliza la transformada rápida de Fourier (FFT) para calcular el espectro de la señal Wave. El resultado es una instancia de la clase Spectrum, que proporciona métodos para manipular y analizar espectros.

La FFT asume que la duración de la señal Wave es un número entero de periodos. Si no es así, pueden ocurrir cosas extrañas. Si la duración no es un múltiplo entero del período, puedes usar el método pad para aumentar la duración hasta que lo sea.

Aquí tienes un ejemplo que agrega relleno a una Wave con una duración de 0.3 segundos:

python
wave = thinkdsp.SquareSignal(440).make_wave(duration=0.3) spectrum = wave.make_spectrum() spectrum.plot() wave.plot()

La figura superior en la Figura 2.7 muestra el espectro de la Wave original. La figura inferior muestra el mismo espectro después de agregar relleno.

El parámetro self es un objeto Wave. n es el número de muestras en la Wave y d es el inverso de la frecuencia de cuadro, que es el tiempo entre muestras.

np.fft es el módulo de NumPy que proporciona funciones relacionadas con la Transformada Rápida de Fourier (FFT), que es un algoritmo eficiente que calcula la Transformada Discreta de Fourier (DFT).

make_spectrum utiliza rfft, que significa "FFT real", porque la Wave contiene valores reales, no complejos. Más adelante veremos la FFT completa, que puede manejar señales complejas. El resultado de rfft, al que llamo hs, es una matriz de números complejos de NumPy que representa la amplitud y el desplazamiento de fase de cada componente de frecuencia en la señal.

El resultado de rfftfreq, al que llamo fs, es una matriz que contiene las frecuencias correspondientes a las hs.

Para entender los valores en hs, considera estas dos formas de pensar en números complejos:

  • Un número complejo es la suma de una parte real y una parte imaginaria, a menudo escrita como x + iy, donde i es la unidad imaginaria, √(-1). Puedes pensar en x e y como coordenadas cartesianas.
  • Un número complejo también es el producto de una magnitud y un exponencial complejo, Aeiφ, donde A es la magnitud y φ es el ángulo en radianes, también llamado "argumento". Puedes pensar en A y φ como coordenadas polares.

Cada valor en hs corresponde a un componente de frecuencia: su magnitud es proporcional a la amplitud del componente correspondiente; su ángulo es el desplazamiento de fase.

La clase Spectrum proporciona dos propiedades de solo lectura, amps y angles, que devuelven matrices de NumPy que representan las magnitudes y ángulos de las hs. Cuando representamos gráficamente un objeto Spectrum, generalmente graficamos amps versus fs. A veces también es útil graficar angles versus fs.

Aunque puede resultar tentador examinar las partes real e imaginaria de hs, casi nunca necesitarás hacerlo. Te animo a pensar en la DFT como una herramienta para analizar las magnitudes y fases de los componentes de frecuencia de una señ


El método make_spectrum utiliza la transformada rápida de Fourier (FFT, por sus siglas en inglés) para calcular el espectro de la onda. El resultado es una instancia de la clase Spectrum, que proporciona métodos para manipular y analizar espectros.

La FFT asume que la duración de la onda es un número entero de períodos. Si no lo es, pueden ocurrir cosas extrañas. Si la duración no es un múltiplo entero del período, puedes usar el método pad para aumentar la duración hasta que lo sea.

Aquí tienes un ejemplo que añade relleno a una onda con una duración de 0.3 segundos:

python
wave = thinkdsp.SquareSignal(440).make_wave(duration=0.3) spectrum = wave.make_spectrum() spectrum.plot() wave.plot()

La figura superior en la Figura 2.7 muestra el espectro de la onda original. La figura inferior muestra el mismo espectro después de añadir relleno.

El parámetro self es un objeto Wave. n es el número de muestras en la onda, y d es el inverso de la frecuencia de fotogramas, que es el tiempo entre muestras.

np.fft es el módulo NumPy que proporciona funciones relacionadas con la Transformada Rápida de Fourier (FFT), que es un algoritmo eficiente que calcula la Transformada de Fourier Discreta (DFT).

make_spectrum utiliza rfft, que significa "FFT real", porque la onda contiene valores reales, no complejos. Más adelante veremos la FFT completa, que puede manejar señales complejas. El resultado de rfft, que llamo hs, es una matriz NumPy de números complejos que representa la amplitud y el desfase de cada componente de frecuencia en la onda.

El resultado de rfftfreq, que llamo fs, es una matriz que contiene las frecuencias correspondientes a las componentes hs.

Para entender los valores en hs, considera estas dos formas de pensar en los números complejos:

  • Un número complejo es la suma de una parte real y una parte imaginaria, a menudo escrita como x + iy, donde i es la unidad imaginaria, √(-1). Puedes pensar en x e y como coordenadas cartesianas.

  • Un número complejo también es el producto de una magnitud y un exponencial complejo, Aeiφ, donde A es la magnitud y φ es el ángulo en radianes, también llamado "argumento". Puedes pensar en A y φ como coordenadas polares.

Cada valor en hs corresponde a una componente de frecuencia: su magnitud es proporcional a la amplitud de la componente correspondiente; su ángulo es el desfase de fase.

La clase Spectrum proporciona dos propiedades de solo lectura, amps y angles, que devuelven matrices NumPy que representan las magnitudes y los ángulos de hs. Cuando trazamos un objeto Spectrum, generalmente trazamos amps versus fs. A veces también es útil trazar angles versus fs.

Aunque puede ser tentador mirar las partes real e imaginaria de hs, casi nunca necesitarás hacerlo. Te animo a que pienses en la DFT como una representación compleja que captura tanto la amplitud como la fase de las componentes de frecuencia

Un "vector de amplitudes y desfases que están codificados en forma de números complejos".

Para modificar un Spectrum, puedes acceder directamente a los hs. Por ejemplo:

python
spectrum.hs *= 2
spectrum.hs[spectrum.fs > cutoff] = 0

La primera línea multiplica los elementos de hs por 2, lo que duplica las amplitudes de todas las componentes. La segunda línea establece en 0 solo los elementos de hs donde la frecuencia correspondiente excede cierta frecuencia de corte.

Pero Spectrum también proporciona métodos para realizar estas operaciones:

python
spectrum.scale(2) spectrum.low_pass(cutoff)

Puedes leer la documentación de estos métodos y otros en http://greenteapress.com/thinkdsp.html.

En este punto, deberías tener una mejor idea de cómo funcionan las clases Signal, Wave y Spectrum, pero no he explicado cómo funciona la Transformada Rápida de Fourier. Eso requerirá algunos capítulos más.


Ejercicios

Las soluciones a estos ejercicios se encuentran en chap02soln.ipynb.
Ejercicio 2.1 Si usa Jupyter, cargue chap02.ipynb y pruebe los ejemplos. También puede ver el cuaderno en http://tinyurl.com/thinkdsp02.

Ejercicio 2.2 Una señal de diente de sierra tiene una forma de onda que aumenta linealmente desde
-1 a 1, luego cae a -1 y se repite. Ver http://en.wikipedia.org/wiki/
Onda_diente_de_sierra
Escriba una clase llamada SawtoothSignal que amplíe Signal y proporcione
evaluar para evaluar una señal de diente de sierra.
Calcular el espectro de una onda de diente de sierra. ¿Cómo es la estructura armónica?
comparar con ondas triangulares y cuadradas?
Para crear una señal de diente de sierra en Python usando la biblioteca thinkdsp, puede usar la clase SawtoothSignal. Así es como puede generar y trazar una señal de diente de sierra:
import numpy as np 
import matplotlib.pyplot as plt 
from thinkdsp import SawtoothSignal # Create a sawtooth signal with frequency 440 Hz 
signal = SawtoothSignal(freq=440) duration = 1.0 # Duration of the signal in seconds 
wave = signal.make_wave(duration=duration, framerate=44100) # Plot the waveform wave.plot()
 plt.xlabel('Time (s)')
 plt.ylabel('Amplitude')
 plt.title('Sawtooth Signal')
 plt.show()


Una señal de diente de sierra tiene una forma de onda que se eleva linealmente desde -1 hasta 1, luego cae a -1 y se repite. Puedes obtener más información sobre la señal de diente de sierra en http://en.wikipedia.org/wiki/Sawtooth_wave

Para resolver este ejercicio, necesitas escribir una clase llamada SawtoothSignal que extienda la clase Signal y proporcione un método evaluate para evaluar una señal de diente de sierra.

Aquí tienes un ejemplo de cómo puedes implementar la clase SawtoothSignal:

python
class SawtoothSignal(Signal): def evaluate(self, ts): cycles = self.freq * ts + self.offset / PI2 frac, _ = np.modf(cycles) ys = normalize(frac, self.amp) return ys

En el método evaluate, ts es la secuencia de tiempos de muestra en los que deseas evaluar la señal. Utilizando la fórmula para generar una señal de diente de sierra, calculamos el número de ciclos cycles desde el tiempo de inicio. Luego utilizamos np.modf para obtener la parte fraccionaria de cycles. La función normalize se encarga de escalar los valores resultantes a la amplitud deseada.

Una vez que hayas implementado la clase SawtoothSignal, puedes crear una instancia de ella y generar la forma de onda y el espectro correspondiente de la siguiente manera:

python
sawtooth = SawtoothSignal(freq=200) wave = sawtooth.make_wave(duration=0.5, framerate=10000) spectrum = wave.make_spectrum() # Plot the waveform wave.plot() # Plot the spectrum spectrum.plot()

Comparando la estructura armónica de la señal de diente de sierra con las señales de triángulo y cuadrada, encontrarás que la señal de diente de sierra contiene todas las armónicas (tanto pares como impares). En contraste, la señal de triángulo solo contiene armónicas impares, y la señal cuadrada solo contiene armónicas impares. Esto se debe a las diferencias en las formas de onda y cómo se generan las señales.

Ejercicio 2.3 Haz una señal cuadrada a 1100 Hz y haz una onda que muestree
a 10000 fotogramas por segundo. Si trazas el espectro, puedes ver que la mayoría
de los armónicos tienen alias. Cuando escuchas la ola, ¿puedes oír el
armónicos alias? 

Exercise 2.3 - Crea una señal cuadrada a 1100 Hz y genera una forma de onda que la muestre a 10000 cuadros por segundo. Si trazas el espectro, podrás ver que la mayoría de las armónicas están aliadas. Cuando escuches la forma de onda, ¿puedes oír las armónicas aliadas?

Para resolver este ejercicio, puedes utilizar la clase SquareSignal proporcionada en el libro y generar una forma de onda y un espectro a partir de ella. Aquí tienes un ejemplo de cómo hacerlo:

python
square = SquareSignal(freq=1100
wave = square.make_wave(duration=0.5, framerate=10000
spectrum = wave.make_spectrum() # Plot the spectrum 
spectrum.plot() # Listen to the wave wave.make_audio()

Al trazar el espectro, podrás observar que hay armónicas aliadas presentes. Esto significa que hay componentes de frecuencia que no corresponden a las frecuencias reales de la señal original, sino que han sido doblados o plegados debido al aliasing.

Cuando escuchas la forma de onda generada, es posible que puedas percibir los efectos del aliasing como tonos adicionales o distorsiones en el sonido. Dependiendo de la calidad de reproducción de audio y de tus habilidades auditivas, es posible que puedas identificar algunos de los armónicos aliados en la señal. Sin embargo, la percepción de los armónicos aliados puede variar de una persona a otra y dependerá de varios factores, como la precisión de la reproducción de audio y la sensibilidad auditiva individual.


Ejercicio 2.4 Si tienes un objeto de espectro, espectro e imprime los primeros
valores de espectro.fs, verá que comienzan en cero. Entonces espectro.hs[0]
es la magnitud del componente con frecuencia 0. Pero, ¿qué significa eso?
¿significar?
Prueba este experimento:
1. Hacer una señal triangular con frecuencia 440 y hacer una Onda con duración 0.01 segundos. Trazar la forma de onda.
2. Cree un objeto Spectrum e imprima espectro.hs[0]. ¿Cuál es la amplitud y la fase de este componente?
3. Configure el espectro.hs[0] = 100. Haga una onda a partir del espectro modificado
y trazarlo. ¿Qué efecto tiene esta operación sobre la forma de onda?
Ejercicio 2.5 Escribe una función que tome un Espectro como parámetro y
lo modifica dividiendo cada elemento de hs por la frecuencia correspondiente
de fs. Sugerencia: dado que la división por cero no está definida, es posible que desee establecer
espectro.hs[0] = 0.
Pruebe su función usando una onda cuadrada, triangular o de diente de sierra.
1. Calcular el Espectro y trazarlo.
2. Modifique el Espectro usando su función y vuelva a trazarlo.
3. Haz una onda del espectro modificado y escúchala. Que efecto
¿Tiene esta operación en la señal?
Ejercicio 2.6 Las ondas triangulares y cuadradas solo tienen armónicos impares; la onda de diente de sierra tiene armónicos pares e impares. Los armónicos del cuadrado.
y las ondas de diente de sierra caen en proporción a 1/f; los armónicos de la onda triangular caen como 1/f 2
. ¿Puedes encontrar una forma de onda que tenga pares e impares?
armónicos que caen como 1/f 2
?
Sugerencia: hay dos formas de abordar esto: podría construir el
señal que desea sumando sinusoides, o podría comenzar con una señal que
es similar a lo que quieres y modifícalo.

Exercise 2.4 - Si tienes un objeto de espectro spectrum y imprimes los primeros valores de spectrum.fs, verás que comienzan en cero. Por lo tanto, spectrum.hs[0] es la magnitud del componente con frecuencia 0. Pero, ¿qué significa eso?

Realiza el siguiente experimento:

  1. Crea una señal triangular con frecuencia de 440 Hz y genera una forma de onda con una duración de 0.01 segundos. Grafica la forma de onda.
  2. Crea un objeto Spectrum y imprime spectrum.hs[0]. ¿Cuál es la amplitud y fase de este componente?
  3. Establece spectrum.hs[0] = 100. Crea una forma de onda a partir del Spectrum modificado y grafícala. ¿Qué efecto tiene esta operación en la forma de onda?

Aquí tienes un ejemplo de cómo realizar este experimento:

python
triangle = TriangleSignal(freq=440
wave = triangle.make_wave(duration=0.01) # Plot the waveform 
wave.plot() # Create the spectrum 
spectrum = wave.make_spectrum() # Print the amplitude and phase of the component at frequency 0 
print(spectrum.hs[0]) # Modify the amplitude of the component at frequency 0 spectrum.hs[0] = 100 # Create a new wave from the modified spectrum 
modified_wave = spectrum.make_wave() # Plot the modified waveform modified_wave.plot()

Al realizar este experimento, descubrirás que el componente con frecuencia 0 en el espectro corresponde a la componente de corriente continua (DC) de la señal. Tener una amplitud y fase en el componente de frecuencia 0 significa que hay una componente de corriente continua en la señal original.

Al establecer spectrum.hs[0] = 100, estás aumentando la amplitud de la componente de frecuencia 0 en el espectro. Al crear una nueva forma de onda a partir del espectro modificado y graficarla, notarás que la forma de onda se ve afectada al tener un componente de corriente continua más grande. En este caso, la amplitud de la forma de onda se incrementará en 100 unidades en relación con la forma de onda original.

Ten en cuenta que, en la práctica, la mayoría de las señales de audio no tienen una componente de corriente continua significativa, por lo que el componente en spectrum.hs[0] suele ser cercano a cero

Ejercicio 2.5 Escribe una función que tome un Espectro como parámetro y
lo modifica dividiendo cada elemento de hs por la frecuencia correspondiente
de fs. Sugerencia: dado que la división por cero no está definida, es posible que desee establecer
espectro.hs[0] = 0.
Pruebe su función usando una onda cuadrada, triangular o de diente de sierra.
1. Calcular el Espectro y trazarlo.
2. Modifique el Espectro usando su función y vuelva a trazarlo.
3. Haz una onda del espectro modificado y escúchala. Que efecto
¿Tiene esta operación en la señal?

Exercise 2.5 - Escribe una función que tome un objeto Spectrum como parámetro y lo modifique dividiendo cada elemento de hs por la frecuencia correspondiente de fs. Pista: como la división por cero no está definida, puedes establecer spectrum.hs[0] = 0.

Prueba tu función usando una señal cuadrada, triangular o diente de sierra.

  1. Calcula el espectro y grafícalo.
  2. Modifica el espectro usando tu función y grafícalo nuevamente.
  3. Crea una forma de onda a partir del espectro modificado y escúchala. ¿Qué efecto tiene esta operación en la señal?

Aquí tienes un ejemplo de cómo implementar esta función:

python
def modify_spectrum(spectrum):
 spectrum.hs /= spectrum.fs spectrum.hs[0] = 0 # Set the first element to 0 to avoid division by zero # Create a square signal 
square = SquareSignal(freq=440) wave = square.make_wave(duration=1.0) # Compute and plot the spectrum 
spectrum = wave.make_spectrum() spectrum.plot() # Modify the spectrum modify_spectrum(spectrum) # Plot the modified spectrum 
spectrum.plot() # Create a wave from the modified spectrum 
modified_wave = spectrum.make_wave() # Listen to the modified wave modified_wave.make_audio()

Al ejecutar este código, verás que la función modify_spectrum divide cada elemento en hs por la frecuencia correspondiente en fs y establece el primer elemento en 0 para evitar la división por cero.

Después de modificar el espectro y graficarlo, notarás que los componentes de frecuencia más altos tienen amplitudes más bajas en comparación con los componentes de frecuencia más bajas.

Luego, al crear una forma de onda a partir del espectro modificado y escucharla, notarás que la señal resultante tiene un cambio en su timbre. En general, la operación de dividir cada elemento por su frecuencia reduce la amplitud de los componentes de frecuencia más altos, lo que puede resultar en un sonido más suave o menos agudo.


Ejercicio 2.6 Las ondas triangulares y cuadradas solo tienen armónicos impares; la onda de diente de sierra tiene armónicos pares e impares. Los armónicos del cuadrado.
y las ondas de diente de sierra caen en proporción a 1/f; los armónicos de la onda triangular caen como 1/f 2
. ¿Puedes encontrar una forma de onda que tenga pares e impares?
armónicos que caen como 1/f 2
?
Sugerencia: hay dos formas de abordar esto: podría construir el
señal que desea sumando sinusoides, o podría comenzar con una señal que
es similar a lo que quieres y modifícalo

Exercise 2.6 - La onda triangular y la onda cuadrada tienen armónicos impares únicamente; la onda diente de sierra tiene tanto armónicos pares como impares. Los armónicos de la onda cuadrada y la onda diente de sierra disminuyen en proporción a 1/f; los armónicos de la onda triangular disminuyen como 1/f^2. ¿Puedes encontrar una forma de onda que tenga armónicos pares e impares que disminuyan como 1/f^2?

Pista: Hay dos formas en las que podrías abordar esto: puedes construir la señal que deseas sumando sinusoides, o puedes empezar con una señal que sea similar a lo que deseas y modificarla.

Aquí tienes un ejemplo de cómo podrías construir una forma de onda que tenga armónicos pares e impares que disminuyan como 1/f^2:

python
import numpy as np 
import matplotlib.pyplot as plt # Define the parameters 
duration = 1.0 # Duration of the waveform in seconds 
framerate = 44100 # Frame rate in frames per second 
frequency = 440 # Fundamental frequency of the waveform # Create the time array 
t = np.linspace(0, duration, int(duration * framerate), endpoint=False) # Initialize the waveform as an array of zeros 
waveform = np.zeros_like(t) # Add up sinusoids with frequencies that decrease as 1/f^2 
for n in range(1, 11): 
 f = frequency / (n ** 2) # Compute the frequency of the n-th harmonic 
 waveform += np.sin(2 * np.pi * f * t) # Normalize the waveform to have maximum amplitude of 1
waveform /= np.max(np.abs(waveform)) # Plot the waveform 
plt.plot(t, waveform) 
plt.xlabel('Time (s)'
plt.ylabel('Amplitude') plt.title('Custom Waveform') plt.show()

En este ejemplo, se construye una forma de onda sumando sinusoides con frecuencias que disminuyen como 1/f^2, donde f es la frecuencia fundamental. El rango de sumatoria se limita a los primeros 10 armónicos para simplificar el ejemplo, pero puedes ajustar este valor según tus necesidades.

Al graficar la forma de onda resultante, verás que tiene tanto armónicos pares como impares y que su amplitud disminuye en proporción a 1/f^2. Puedes experimentar con diferentes valores de frecuencia fundamental y rango de sumatoria para obtener diferentes formas de onda con esta característica.

Resumen del Capítulo 2:

El Capítulo 2 del libro "Think DSP" aborda los conceptos fundamentales de las señales y su representación en el dominio de la frecuencia utilizando la Transformada de Fourier. A lo largo del capítulo, se introducen las clases Signal, Wave y Spectrum para representar y manipular señales.

Las partes más importantes del Capítulo 2 son las siguientes:

  1. Señales: Se explican las diferentes características y tipos de señales, como las señales sinusoidales, triangulares, cuadradas y diente de sierra. Se muestra cómo generar y visualizar estas señales utilizando la librería thinkdsp.

  2. Espectros: Se introduce el concepto de espectro, que representa la distribución de frecuencias en una señal. Se demuestra cómo calcular el espectro de una señal utilizando la Transformada de Fourier.

  3. Aliasing: Se explica el fenómeno del aliasing, que ocurre cuando una señal se muestrea a una frecuencia insuficiente y las frecuencias altas se pliegan hacia frecuencias más bajas. Se ilustra este efecto y se discuten sus implicaciones.

  4. Cálculo del espectro: Se presenta la función make_spectrum() que calcula el espectro de una señal. Se explica cómo acceder y modificar los componentes del espectro, como las amplitudes y las fases.

  5. Ejercicios: Se proponen una serie de ejercicios prácticos para aplicar los conceptos aprendidos, como la implementación de una señal de onda diente de sierra, la modificación de espectros y la exploración de diferentes formas de onda.

En general, el Capítulo 2 sienta las bases para comprender cómo se representan y manipulan las señales en el dominio de la frecuencia. Proporciona una introducción sólida a los conceptos clave necesarios para el análisis y procesamiento de señales digitales.

Capítulo 3: Señales no periódicas

Las señales con las que hemos trabajado hasta ahora son periódicas, lo que significa que se repiten infinitamente. También significa que los componentes de frecuencia que contienen no cambian con el tiempo. En este capítulo, consideramos señales no periódicas, cuyos componentes de frecuencia sí cambian a lo largo del tiempo. En otras palabras, prácticamente todas las señales de sonido.

Este capítulo también presenta espectrogramas, una forma común de visualizar señales no periódicas.

El código correspondiente a este capítulo se encuentra en el archivo chap03.ipynb, que está en el repositorio de este libro (consulte la Sección 0.2). También puedes verlo en http://tinyurl.com/thinkdsp03.

3.1 Chirp lineal

Comenzaremos con un "chirp", que es una señal con una frecuencia variable. thinkdsp proporciona una clase llamada Chirp que genera un senoide que varía linealmente a través de un rango de frecuencias.

Aquí tienes un ejemplo que varía desde 220 hasta 880 Hz, lo que equivale a dos octavas desde A3 hasta A5:

python
signal = thinkdsp.Chirp(start=220, end=880) wave = signal.make_wave()

La Figura 3.1 muestra segmentos de esta onda cerca del inicio, en el medio y al final. Es evidente que la frecuencia está aumentando.

Antes de continuar, veamos cómo se implementa la clase Chirp. Aquí está la definición de la clase:

Figura 3.1: Forma de onda de un chirp cerca del inicio, en el medio y al final.

python
class Chirp(Signal): 
def __init__(self, start=440, end=880, amp=1.0):
 self.start = start self.end = end self.amp = amp def evaluate(self, ts):
 freqs = np.linspace(self.start, self.end, len(ts)) 
 dts = np.diff(ts, prepend=0
 dphis = PI2 * freqs * dts phases = np.cumsum(dphis) 
 ys = self.amp * np.cos(phases) 
return ys

start y end son las frecuencias, en Hz, al inicio y al final del chirp. amp es la amplitud.

Aquí está la función que evalúa la señal:

python
def evaluate(self, ts):
 freqs = np.linspace(self.start, self.end, len(ts)) 
 dts = np.diff(ts, prepend=0
 dphis = PI2 * freqs * dts 
 phases = np.cumsum(dphis) 
 ys = self.amp * np.cos(phases)
return ys

ts es la secuencia de puntos en el tiempo donde se debe evaluar la señal. Para mantener esta función simple, asumo que los puntos están equiespaciados.

Para calcular la frecuencia en cada punto en el tiempo, uso np.linspace, que devuelve un array de NumPy con n valores entre start y end.

np.diff calcula la diferencia entre elementos adyacentes de ts, devolviendo la longitud de cada intervalo en segundos. Si los elementos de ts están equiespaciados, los dts son todos iguales.

El siguiente paso es determinar cuánto cambia la fase durante cada intervalo. En la Sección 1.7 vimos que cuando la frecuencia es constante, la fase, φ, aumenta linealmente con el tiempo:

φ = 2πf t

Cuando la frecuencia es una función del tiempo, el cambio en la fase durante un intervalo de tiempo corto, ∆t, es:

∆φ = 2πf(t)∆t

En Python, dado que freqs contiene f(t) y dts contiene los intervalos de tiempo, podemos escribir:

dphis = PI2 * freqs * dts

Ahora, dado que dphis contiene los cambios en la fase, podemos obtener la fase total en cada paso de tiempo sumando los cambios:

phases = np.cumsum(dphis) phases = np.insert(phases, 0, 0)

np.cumsum calcula la suma acumulativa, que es casi lo que queremos, pero no comienza en 0. Por eso uso np.insert para agregar un 0 al principio.

El resultado es un array de NumPy donde el i-ésimo elemento contiene la suma de los primeros i términos de dphis; es decir, la fase total al final del i-ésimo intervalo.

Finalmente, np.cos calcula la amplitud de la onda como función de la fase (recuerda que la fase se expresa en radianes).

Si conoces cálculo, es posible que notes que el límite cuando ∆t se vuelve pequeño es:

dφ = 2πf(t)dt

Dividiendo por dt obtenemos:

dφ dt = 2πf(t)

En otras palabras, la frecuencia es la derivada de la fase. Recíprocamente, la fase es la integral de la frecuencia. Cuando usamos cumsum para pasar de frecuencia a fase, estábamos aproximando la integración.

3.2 Chirp exponencial

Cuando escuchas este chirp, es posible que notes que el tono aumenta rápidamente al principio y luego se desacelera. El chirp abarca dos octavas, pero solo toma 2/3 de segundo para abarcar la primera octava, y el doble de tiempo para abarcar la segunda.

La razón es que nuestra percepción del tono depende del logaritmo de la frecuencia. Como resultado, el intervalo que escuchamos entre dos notas depende de la razón de sus frecuencias, no de la diferencia. "Intervalo" es el término musical para la diferencia percibida entre dos tonos.

Por ejemplo, una octava es un intervalo donde la razón de dos tonos es 2. Entonces, el intervalo de 220 a 440 es una octava y el intervalo de 440 a 880 también es una octava. La diferencia en frecuencia es mayor, pero la razón es la misma.

Como resultado, si la frecuencia aumenta linealmente, como en un chirp lineal, el tono percibido aumenta de manera logarítmica.

Si deseas que el tono percibido aumente linealmente, la frecuencia debe aumentar exponencialmente. Una señal con esa forma se llama chirp exponencial.

Aquí está el código que crea uno:

python
class ExpoChirp(Chirp):
def evaluate(self, ts):
 start, end = np.log10(self.start), np.log10(self.end) 
 freqs = np.logspace(start, end, len(ts)-1
return self._evaluate(ts, freqs)

En lugar de np.linspace, esta versión de evaluate usa np.logspace, que crea una serie de frecuencias cuyos logaritmos están igualmente espaciados, lo que significa que aumentan exponencialmente.

Eso es todo; todo lo demás es igual que en el Chirp. Aquí está el código que crea uno:

python
signal = thinkdsp.ExpoChirp(start=220, end=880) wave = signal.make_wave(duration=1)

Puedes escuchar estos ejemplos en el archivo chap03.ipynb y comparar los chirps lineales y exponenciales.

3.3 Espectro de un chirp

¿Qué crees que sucede si calculamos el espectro de un chirp? Aquí tienes un ejemplo que construye un chirp de una octava y un segundo de duración y su espectro:

python
signal = thinkdsp.Chirp(start=220, end=440) wave = signal.make_wave(duration=1) spectrum = wave.make_spectrum()
Figura 3.2: Espectro de un chirp de una octava y un segundo de duración.

La Figura 3.2 muestra el resultado. El espectro tiene componentes en cada frecuencia de 220 a 440 Hz, con variaciones que se asemejan un poco al Ojo de Sauron (ver http://en.wikipedia.org/wiki/Sauron).

El espectro es aproximadamente plano entre 220 y 440 Hz, lo que indica que la señal pasa igual cantidad de tiempo en cada frecuencia dentro de este rango. Basándonos en esa observación, deberías poder adivinar cómo se ve el espectro de un chirp exponencial.

El espectro proporciona pistas sobre la estructura de la señal, pero oculta la relación entre frecuencia y tiempo. Por ejemplo, no podemos determinar viendo este espectro si la frecuencia aumentó, disminuyó o ambas cosas.

3.4 Espectrograma

Para recuperar la relación entre frecuencia y tiempo, podemos dividir el chirp en segmentos y graficar el espectro de cada segmento. El resultado se llama transformada de Fourier de tiempo corto (STFT, por sus siglas en inglés).

Existen varias formas de visualizar una STFT, pero la más común es un espectrograma, que muestra el tiempo en el eje x y la frecuencia en el eje y. Cada columna en el espectrograma muestra el espectro de un segmento corto, utilizando color o escala de grises para representar la amplitud.

Como ejemplo, calcularé el espectrograma de este chirp:

python
signal = thinkdsp.Chirp(start=220, end=440) wave = signal.make_wave(duration=1, framerate=11025)
Figura 3.3: Espectrograma de un chirp de una octava y un segundo de duración.

Wave proporciona la función make_spectrogram, que devuelve un objeto Spectrogram:

python
spectrogram = wave.make_spectrogram(seg_length=512) spectrogram.plot(high=700)

seg_length es el número de muestras en cada segmento. Elegí 512 porque la FFT es más eficiente cuando el número de muestras es una potencia de 2.

La Figura 3.3 muestra el resultado. El eje x muestra el tiempo de 0 a 1 segundo. El eje y muestra la frecuencia de 0 a 700 Hz. Corté la parte superior del espectrograma; el rango completo llega hasta 5512.5 Hz, que es la mitad de la velocidad de muestreo.

El espectrograma muestra claramente que la frecuencia aumenta linealmente con el tiempo. De manera similar, en el espectrograma de un chirp exponencial, podemos ver la forma de la curva exponencial.

Sin embargo, debes notar que el pico en cada columna está difuminado en 2-3 celdas. Este difuminado refleja la resolución limitada del espectrograma.

3.5 El límite de Gabor

La resolución temporal del espectrograma es la duración de los segmentos, que corresponde al ancho de las celdas en el espectrograma. Dado que cada segmento tiene 512 fotogramas y hay 11,025 fotogramas por segundo, la duración de cada segmento es de aproximadamente 0.046 segundos.

La resolución de frecuencia es el rango de frecuencia entre elementos en el espectro, que corresponde a la altura de las celdas. Con 512 fotogramas, tenemos

Para obtener 256 componentes de frecuencia en un rango de 0 a 5512.5 Hz, el rango entre componentes es de 21.6 Hz.

Más generalmente, si n es la longitud del segmento, el espectro contiene n/2 componentes. Si la velocidad de muestreo es r, la frecuencia máxima en el espectro es r/2. Por lo tanto, la resolución temporal es n/r y la resolución de frecuencia es r/n, que es r/n.

Idealmente, nos gustaría que la resolución temporal fuera pequeña para poder ver cambios rápidos en la frecuencia. Y nos gustaría que la resolución de frecuencia fuera pequeña para poder ver cambios pequeños en la frecuencia. Pero no se puede tener ambas cosas. Observa que la resolución temporal, n/r, es el inverso de la resolución de frecuencia, r/n. Por lo tanto, si una se hace más pequeña, la otra se hace más grande.

Por ejemplo, si duplicas la longitud del segmento, reduces a la mitad la resolución de frecuencia (lo cual es bueno), pero duplicas la resolución temporal (lo cual es malo). Incluso aumentar la velocidad de muestreo no ayuda. Obtienes más muestras, pero el rango de frecuencias también aumenta al mismo tiempo.

Este compromiso se llama límite de Gabor y es una limitación fundamental de este tipo de análisis tiempo-frecuencia.

3.6 Leakage

Para explicar cómo funciona make_spectrogram, tengo que explicar el uso de ventanas (windowing). Y para explicar las ventanas, debo mostrarte el problema que intentan solucionar, que es la fuga (leakage).

La Transformada Discreta de Fourier (DFT), que utilizamos para calcular los espectros, trata a las ondas como si fueran periódicas; es decir, asume que el segmento finito en el que opera es un periodo completo de una señal infinita que se repite a lo largo de todo el tiempo. En la práctica, esta suposición a menudo es falsa, lo que crea problemas.

Un problema común es la discontinuidad al comienzo y al final del segmento. Debido a que la DFT asume que la señal es periódica, conecta implícitamente el final del segmento con el comienzo para formar un bucle. Si el final no se conecta suavemente con el comienzo, la discontinuidad crea componentes de frecuencia adicionales en el segmento que no están en la señal.

Como ejemplo, comencemos con una señal senoidal que contiene solo un componente de frecuencia de 440 Hz.

Figura 3.4: Espectro de un segmento periódico de una sinusoide (izquierda), un segmento no periódico (centro), un segmento no periódico con ventana (derecha). 
 Para ilustrar el fenómeno de fuga (leakage), utilizaremos la señal senoidal de frecuencia 440 Hz. Si seleccionamos un segmento que sea un múltiplo entero del período de la señal, el final del segmento se conectará suavemente con el principio y la DFT se comportará correctamente.
python
signal = thinkdsp.SinSignal(freq=440
duration = signal.period * 30 
wave = signal.make_wave(duration) 
spectrum = wave.make_spectrum()

La Figura 3.4 (izquierda) muestra el resultado. Como se esperaba, hay un único pico en 440 Hz.

Pero si la duración no es un múltiplo del período, ocurren cosas no deseadas. Con duration = signal.period * 30.25, la señal comienza en 0 y termina en 1.

La Figura 3.4 (centro) muestra el espectro de este segmento. Nuevamente, el pico está en 440 Hz, pero ahora hay componentes adicionales distribuidos desde 240 hasta 640 Hz. Esta distribución se llama fuga espectral (spectral leakage), porque parte de la energía que está en la frecuencia fundamental se filtra hacia otras frecuencias.

En este ejemplo, la fuga ocurre porque estamos utilizando la DFT en un segmento que se vuelve discontinuo cuando lo tratamos como periódico.

Podemos reducir la fuga (leakage) suavizando la discontinuidad entre el principio y el final del segmento, y una forma de hacerlo es mediante el uso de ventanas (windowing).

3.7 Windowing

Podemos reducir las fugas suavizando la discontinuidad entre el principio y el final del segmento, y una forma de hacerlo es creando ventanas.

Una "ventana" es una función diseñada para transformar un segmento no periódico en algo que pueda parecer periódico. La Figura 3.5 (arriba) muestra un segmento donde el final no se conecta suavemente con el principio.

La Figura 3.5 (centro) muestra una "ventana Hamming", una de las funciones de ventana más comunes. Ninguna función de ventana es perfecta, pero algunas se pueden mostrar como óptimas para diferentes aplicaciones, y Hamming es una buena ventana de uso general.

La Figura 3.5 (abajo) muestra el resultado de multiplicar la ventana por la señal original. Donde la ventana es cercana a 1, la señal no se altera. Donde la ventana es cercana a 0, la señal se atenúa. Debido a que la ventana se estrecha en ambos extremos, el final del segmento se conecta suavemente con el principio.

La Figura 3.4 (derecha) muestra el espectro de la señal con ventana. El uso de la ventana ha reducido sustancialmente la fuga, pero no completamente.
Figura 3.6: Ventanas de Hamming superpuestas.

Aquí tienes el código. Wave proporciona el método window, que aplica una ventana de Hamming:

python
#class Wave: def window(self, window): self.ys *= window

Y NumPy proporciona la función hamming, que calcula una ventana de Hamming con una longitud dada:

python
window = np.hamming(len(wave)) wave.window(window)

NumPy proporciona funciones para calcular otras funciones de ventana, incluyendo bartlett, blackman, hanning y kaiser. Uno de los ejercicios al final de este capítulo te pide que experimentes con estas otras ventanas.

Ahora que entendemos el uso de ventanas, podemos entender la implementación del espectrograma. Aquí está el método de la clase Wave que calcula espectrogramas:

python
#class Wave: def make_spectrogram(self, seg_length): window = np.hamming(seg_length) i, j = 0, seg_length step = seg_length / 2 spec_map = {}

El método make_spectrogram divide la forma de onda en segmentos solapados utilizando una ventana de Hamming. El tamaño de cada segmento se especifica mediante el parámetro seg_length. Luego, se itera sobre los segmentos con un paso de seg_length/2, aplicando la ventana de Hamming a cada segmento. El resultado es un mapa de espectrograma que contiene los segmentos y sus espectros correspondientes.

while j < len(self.ys):
 segmento = self.slice(i, j) segmento.window(window) 
 t = (segmento.start + segmento.end) / 2 spec_map[t] = segmento.make_spectrum() i += step j += step return Spectrogram(spec_map, seg_length)

Esta es la función más larga del libro, así que si puedes entender esto, puedes entender cualquier cosa.

El parámetro self es un objeto de tipo Wave. seg_length es el número de muestras en cada segmento. window es una ventana de Hamming con la misma longitud que los segmentos. i y j son los índices de corte que seleccionan los segmentos de la forma de onda. step es el desplazamiento entre segmentos. Dado que step es la mitad de seg_length, los segmentos se superponen en la mitad. La Figura 3.6 muestra cómo se ven estas ventanas superpuestas. spec_map es un diccionario que mapea una marca de tiempo a un espectro. Dentro del bucle while, seleccionamos una porción de la forma de onda y aplicamos la ventana; luego construimos un objeto Spectrum y lo agregamos a spec_map. El tiempo nominal de cada segmento, t, es el punto medio. Luego avanzamos i y j, y continuamos siempre que j no supere el final de la Wave. Finalmente, el método construye y devuelve un Spectrogram. Aquí está la definición de la clase Spectrogram:

python
class Spectrogram(object): def __init__(self, spec_map, seg_length): self.spec_map = spec_map self.seg_length = seg_length

Al igual que muchos métodos __init__, este simplemente guarda los parámetros como atributos. Spectrogram proporciona el método plot, que genera una representación en pseudocolor con el tiempo en el eje x y la frecuencia en el eje y.

Y así es como se implementan los espectrogramas.

Ejercicio 3.1: Ejecuta y escucha los ejemplos en chap03.ipynb, que se encuentra en el repositorio de este libro y también está disponible en http://tinyurl.com/thinkdsp03. En el ejemplo de "leakage", intenta reemplazar la ventana de Hamming con una de las otras ventanas proporcionadas por NumPy, y observa qué efecto tienen en el "leakage". Puedes consultar http://docs.scipy.org/doc/numpy/reference/routines.window.html para más información sobre las ventanas proporcionadas por NumPy.

Ejercicio 3.2: Escribe una clase llamada "SawtoothChirp" que extienda la clase "Chirp" y sobrescriba el método "evaluate" para generar una forma de onda de sierra cuya frecuencia aumente (o disminuya) linealmente. Sugerencia: combina las funciones "evaluate" de las clases "Chirp" y "SawtoothSignal". Dibuja un esquema de cómo crees que se verá el espectrograma de esta señal y luego plótalo. El efecto del aliasing debería ser visualmente evidente y si escuchas cuidadosamente, podrás oírlo.

Ejercicio 3.3: Crea una señal de sierra que varíe en frecuencia desde 2500 Hz hasta 3000 Hz y úsala para generar una forma de onda con una duración de 1 segundo y una frecuencia de muestreo de 20 kHz. Dibuja un esquema de cómo crees que se verá el espectro y luego plotea el espectro para verificar si estabas en lo correcto.

Ejercicio 3.4: En terminología musical, un "glissando" es una nota que se desliza desde un tono hasta otro, similar a un chirp. Encuentra o crea una grabación de un glissando y traza el espectrograma de los primeros segundos. Una sugerencia: el comienzo de la Rhapsody in Blue de George Gershwin tiene un famoso glissando de clarinete, que puedes descargar de http://archive.org/details/rhapblue11924.

Ejercicio 3.5: Un trombonista puede realizar un glissando extendiendo el tubo del trombón mientras sopla continuamente. A medida que el tubo se extiende, la longitud total del tubo aumenta y la frecuencia resultante es inversamente proporcional a la longitud. Suponiendo que el trombonista mueve el tubo a una velocidad constante, ¿cómo varía la frecuencia con el tiempo? Escribe una clase llamada "TromboneGliss" que extienda la clase "Chirp" y proporciona el método "evaluate". Crea una forma de onda que simule un glissando de trombón desde C3 hasta F3 y de regreso a C3. C3 es 262 Hz; F3 es 349 Hz. Plotea el espectrograma de la forma de onda resultante. ¿Un glissando de trombón se parece más a un chirp lineal o exponencial?

Ejercicio 3.6: Haz o encuentra una grabación de una serie de sonidos de vocales y observa el espectrograma. ¿Puedes identificar diferentes vocales?


El Capítulo 3 del libro "Think DSP" aborda el análisis espectral y la transformada de Fourier. A continuación se presenta un resumen y algunas cosas relevantes del capítulo:

  1. El análisis espectral es una herramienta fundamental para comprender la estructura y las características de las señales. Permite descomponer una señal en sus componentes de frecuencia y estudiar su contenido espectral.

  2. La transformada de Fourier es una técnica matemática utilizada para realizar el análisis espectral de las señales. Permite descomponer una señal en una combinación de ondas sinusoidales de diferentes frecuencias y amplitudes.

  3. La DFT (Transformada Discreta de Fourier) es una implementación computacional de la transformada de Fourier que se utiliza para analizar señales discretas.

  4. El espectro de una señal es la representación gráfica de su contenido espectral, que muestra la amplitud de las diferentes frecuencias presentes en la señal.

  5. La resolución temporal y la resolución espectral son dos conceptos importantes en el análisis espectral. Existe un compromiso entre la resolución temporal y la resolución espectral: si se aumenta la resolución temporal, se reduce la resolución espectral y viceversa.

  6. El problema de la fuga es una limitación en el análisis espectral que ocurre cuando se aplica la DFT a segmentos de señales no periódicas. Esto puede causar la aparición de componentes de frecuencia adicionales en el espectro que no están presentes en la señal original.

  7. La técnica de ventana (windowing) se utiliza para reducir la fuga espectral. Consiste en multiplicar el segmento de la señal por una función ventana que suaviza la discontinuidad en los extremos del segmento.

  8. Existen diferentes funciones ventana que se pueden utilizar, como la ventana de Hamming, la ventana de Bartlett, la ventana de Blackman, la ventana de Hanning, entre otras. Cada una tiene propiedades y aplicaciones específicas.

  9. Un espectrograma es una representación visual del espectro de una señal en función del tiempo. Se construye dividiendo la señal en segmentos solapados, aplicando una función ventana a cada segmento y calculando el espectro de cada segmento. El espectrograma proporciona información sobre cómo cambia el contenido espectral de la señal a lo largo del tiempo.

  10. El capítulo presenta ejercicios prácticos que permiten aplicar los conceptos aprendidos, como la exploración de diferentes ventanas, la generación de señales de chirp y sierra, el análisis de glissandos y la identificación de vocales en espectrogramas.

En resumen, el Capítulo 3 introduce los conceptos y técnicas fundamentales del análisis espectral y la transformada de Fourier, proporcionando las bases para comprender y analizar las características de las señales en el dominio de la frecuencia.

Capítulo 4: Ruido

En inglés, "ruido" se refiere a un sonido no deseado o desagradable. En el contexto del procesamiento de señales, tiene dos sentidos diferentes:

  1. Como en inglés, puede significar una señal no deseada de cualquier tipo. Si dos señales interfieren entre sí, cada una consideraría a la otra como ruido.

  2. "Ruido" también se refiere a una señal que contiene componentes en muchas frecuencias, por lo que carece de la estructura armónica de las señales periódicas que vimos en capítulos anteriores.

Este capítulo trata sobre el segundo tipo de ruido.

El código correspondiente a este capítulo se encuentra en "chap04.ipynb", que está en el repositorio del libro (consulte la Sección 0.2). También puedes verlo en http://tinyurl.com/thinkdsp04.

4.1 Ruido no correlacionado

La forma más sencilla de entender el ruido es generándolo, y la forma más simple de generarlo es mediante el ruido uniforme no correlacionado (UU noise). "Uniforme" significa que la señal contiene valores aleatorios de una distribución uniforme, es decir, cada valor en el rango es igualmente probable. "No correlacionado" significa que los valores son independientes, es decir, conocer un valor no proporciona información sobre los demás.

Aquí hay una clase que representa el ruido UU:

Figura 4.1: Forma de onda del ruido uniforme no correlacionado.

La clase "UncorrelatedUniformNoise" hereda de "_Noise", que a su vez hereda de "Signal".

Como es habitual, la función "evaluate" toma los tiempos "ts" en los que se evaluará la señal. Utiliza "np.random.uniform", que genera valores de una distribución uniforme. En este ejemplo, los valores están en el rango entre -amp y amp.

El siguiente ejemplo genera ruido UU con una duración de 0.5 segundos y una frecuencia de muestreo de 11,025 muestras por segundo.

python
signal = thinkdsp.UncorrelatedUniformNoise() 
wave = signal.make_wave(duration=0.5, framerate=11025)

Si reproduces esta onda, suena como la estática que se escucha al sintonizar una radio entre canales. La Figura 4.1 muestra cómo se ve la forma de onda. Como era de esperar, parece bastante aleatoria.

Ahora echemos un vistazo al espectro:

python
spectrum = wave.make_spectrum() spectrum.plot_power()

"Spectrum.plot_power" es similar a "Spectrum.plot", excepto que muestra la potencia en lugar de la amplitud. La potencia es el cuadrado de la amplitud. Estoy cambiando de amplitud a potencia en este capítulo porque es más convencional en el contexto del ruido.



Figura 4.2: Espectro de potencia del ruido uniforme no correlacionado.

La Figura 4.2 muestra el resultado. Al igual que la señal, el espectro parece bastante aleatorio. De hecho, es aleatorio, pero tenemos que ser más precisos sobre la palabra "aleatorio".

Hay al menos tres cosas que nos gustaría saber sobre una señal de ruido o su espectro:

  • Distribución: La distribución de una señal aleatoria es el conjunto de valores posibles y sus probabilidades. Por ejemplo, en la señal de ruido uniforme, el conjunto de valores está en el rango de -1 a 1, y todos los valores tienen la misma probabilidad. Una alternativa es el ruido gaussiano, donde el conjunto de valores está en el rango desde negativo hasta positivo infinito, pero los valores cercanos a 0 son los más probables, con una probabilidad que disminuye según la curva gaussiana o en forma de "campana".
  • Correlación: ¿Cada valor en la señal es independiente de los demás, o existen dependencias entre ellos? En el ruido UU, los valores son independientes. Una alternativa es el ruido browniano, donde cada valor es la suma del valor anterior y un "paso" aleatorio. Entonces, si el valor de la señal es alto en un momento particular en el tiempo, esperamos que siga siendo alto, y si es bajo, esperamos que siga siendo bajo.
  • Relación entre potencia y frecuencia: En el espectro del ruido UU, la potencia en todas las frecuencias se extrae de la misma distribución; es decir, la potencia promedio es la misma para todas las frecuencias. Una alternativa es el ruido rosa, donde la potencia está inversamente relacionada con la frecuencia; es decir, la potencia en la frecuencia f se extrae de una distribución cuya media es proporcional a 1/f.




  • Figura 4.3: Espectro integrado de ruido uniforme no correlacionado.
    • 4.2 Espectro integrado

      Para el ruido UU podemos ver más la relación entre potencia y frecuencia claramente mirando el espectro integrado, que es una función de la frecuencia, f, que muestra la potencia acumulada en el espectro hasta f. Spectrum proporciona un método que calcula el IntegratedSpectrum:

      python
      def make_integrated_spectrum(self):

      self.power es un arreglo de NumPy que contiene la potencia para cada frecuencia. np.cumsum calcula la suma acumulativa de las potencias. Al dividir por el último elemento se normaliza el espectro integrado para que vaya de 0 a 1.

      El resultado es un objeto de tipo IntegratedSpectrum. Aquí está la definición de la clase:

      python
      class IntegratedSpectrum(object):

      Al igual que Spectrum, IntegratedSpectrum proporciona plot_power, por lo que podemos calcular y graficar el espectro integrado de la siguiente manera:

      python
      integ = spectrum.make_integrated_spectrum()

      El resultado, mostrado en la Figura 4.3, es una línea recta, lo que indica que la potencia en todas las frecuencias es constante, en promedio. El ruido con potencia igual en todas las frecuencias se conoce como ruido blanco.



    •  integ.plot_power()
    •  cs = np.cumsum(self.power) cs /= cs[-1return IntegratedSpectrum(cs, self.fs)
    • def __init__(self, cs, fs): self.cs = cs self.fs = fs
    • Spectrum proporciona un método que calcula el Espectro Integrado:

    • Para el ruido UU podemos ver la relación entre la potencia y la frecuencia de manera más clara al observar el espectro integrado, que es una función de la frecuencia, f, que muestra la potencia acumulada en el espectro hasta f.



    • Figura 4.4: Forma de onda de ruido Browniano.

      frecuencias se llama ruido blanco por analogía con la luz, porque una mezcla igual de luz en todas las frecuencias visibles es blanca



    • 4.3 Brownian noise

    • Se llama "Browniano" por analogía con el movimiento Browniano, en el que una partícula suspendida en un fluido se mueve aparentemente al azar debido a interacciones invisibles con el fluido. El movimiento Browniano se describe a menudo utilizando una caminata aleatoria, que es un modelo matemático de una trayectoria donde la distancia entre pasos está caracterizada por una distribución aleatoria.

      Ruido browniano, en el que cada valor es la suma del valor anterior y un “paso” aleatorio En una caminata aleatoria unidimensional, la partícula se mueve hacia arriba o hacia abajo en una cantidad aleatoria en cada paso de tiempo. La ubicación de la partícula en cualquier punto en el tiempo es la suma de todos los pasos anteriores.

      Esta observación sugiere una forma de generar ruido Browniano: generar pasos aleatorios no correlacionados y luego sumarlos. Aquí hay una definición de clase que implementa este algoritmo:

      python
      class BrownianNoise(_Noise): def evaluate(self, ts): dys = np.random.uniform(-11len(ts))

      El ruido blanco es aquel en el que la potencia en todas las frecuencias es constante. Por otro lado, el ruido Browniano es una alternativa al ruido blanco en el que cada valor es la suma del valor anterior y un "paso" aleatorio.ys = np.cumsum(dys)

    •  ys = normalize(unbias(ys), self.amp) 
    • return ys 

Figura 4.5: Espectro de ruido browniano en escala lineal (izquierda) y escala logarítmica (derecha). discord

evaluate utiliza np.random.uniform para generar una señal no correlacionada y np.cumsum para calcular su suma acumulada. Dado que la suma probablemente escape del rango de -1 a 1, debemos usar unbias para desplazar la media a 0 y normalize para obtener la amplitud máxima deseada. Aquí está el código que genera un objeto BrownianNoise y traza la forma de onda:

python
signal = thinkdsp.BrownianNoise() wave = signal.make_wave(duration=0.5, framerate=11025) wave.plot()

La Figura 4.4 muestra el resultado. La forma de onda se desplaza hacia arriba y hacia abajo, pero hay una clara correlación entre los valores sucesivos. Cuando la amplitud es alta, tiende a mantenerse alta y viceversa. Si trazas el espectro del ruido Browniano en una escala lineal, como se muestra en la Figura 4.5 (izquierda), no se ve mucho. Casi toda la potencia se encuentra en las frecuencias más bajas; los componentes de frecuencia más alta no son visibles. Para ver la forma del espectro de manera más clara, podemos trazar la potencia y la frecuencia en una escala log-log. Aquí está el código:

python
import matplotlib.pyplot as plt spectrum = wave.make_spectrum()

La siguiente línea de código traza el espectro de potencia en una escala log-log y calcula la pendiente estimada:

python
spectrum.plot_power(linewidth=1, alpha=0.5) plt.xscale('log') plt.yscale('log')

El resultado se muestra en la Figura 4.5 (derecha). La relación entre la potencia y la frecuencia es ruidosa, pero aproximadamente lineal. Spectrum proporciona el método estimate_slope, que utiliza SciPy para calcular un ajuste de mínimos cuadrados al espectro de potencia:

python
def estimate_slope(self): x = np.log(self.fs[1:]) y = np.log(self.power[1:]) t = scipy.stats.linregress(x, y) return t

Descarta el primer componente del espectro porque este componente corresponde a f = 0 y log 0 no está definido. estimate_slope devuelve el resultado de scipy.stats.linregress, que es un objeto que contiene la pendiente estimada, el intercepto, el coeficiente de determinación (R2), el valor p y el error estándar. Para nuestros propósitos, solo necesitamos la pendiente. Para el ruido Browniano, la pendiente del espectro de potencia es -2 (veremos por qué en el Capítulo 9), por lo que podemos escribir esta relación: log P = k - 2 log f donde P es la potencia, f es la frecuencia y k es el intercepto de la línea, que no es importante para nuestros propósitos. Exponenciando ambos lados, obtenemos: P = K / f^2 donde K es e^k, pero sigue sin ser importante. Lo más relevante es que la potencia es proporcional a 1/f^2, lo cual es característico del ruido Browniano. El ruido Browniano también se llama ruido rojo, por la misma razón por la que el ruido blanco se llama "blanco". Si se combina luz visible con potencia proporcional a 1/f^2, la mayor parte de la potencia estaría en el extremo de baja frecuencia del espectro, que es rojo. A veces también se llama "ruido marrón" al ruido Browniano, pero creo que eso es confuso, así que no lo usaré

Diferencias de Ruido

El ruido Browniano, Blanco, Rosa y Rojo son términos utilizados para describir diferentes características y propiedades de los procesos de ruido en señales.

  1. Ruido Blanco: El ruido blanco es un tipo de ruido en el que la potencia es igual en todas las frecuencias dentro de un rango determinado. Esto significa que todas las frecuencias tienen la misma contribución de energía. En el dominio del tiempo, el ruido blanco se ve como una señal aleatoria que no muestra ninguna correlación entre los valores adyacentes. En el espectro de frecuencia, el ruido blanco tiene una distribución plana y constante de potencia. El ruido blanco se llama así por analogía con la luz blanca, que contiene una mezcla igual de todas las frecuencias visibles.

  2. Ruido Rosa: El ruido rosa, también conocido como ruido 1/f, es un tipo de ruido en el que la potencia disminuye a medida que la frecuencia aumenta. En otras palabras, la densidad espectral de potencia es inversamente proporcional a la frecuencia. Esto significa que las frecuencias más bajas tienen una mayor contribución de potencia en comparación con las frecuencias más altas. En el dominio del tiempo, el ruido rosa muestra una cierta correlación entre los valores adyacentes, lo que resulta en una decaída suave en la potencia a medida que aumenta la frecuencia. El ruido rosa se llama así porque se asemeja al color rosa en el espectro de luz visible.

  3. Ruido Rojo: El ruido rojo es similar al ruido rosa en el sentido de que la potencia disminuye a medida que aumenta la frecuencia. Sin embargo, a diferencia del ruido rosa, el ruido rojo tiene una caída más pronunciada en la potencia a medida que aumenta la frecuencia. Esto significa que las frecuencias más bajas tienen una contribución aún mayor de potencia en comparación con las frecuencias más altas. El ruido rojo se llama así porque se asemeja al color rojo en el espectro de luz visible.

  4. Ruido Browniano: El ruido Browniano es un tipo de ruido en el que cada valor en la señal es la suma del valor anterior y un paso aleatorio. Esto resulta en una señal que muestra una correlación clara entre los valores adyacentes. En el espectro de frecuencia, el ruido Browniano tiene una relación de potencia inversamente proporcional al cuadrado de la frecuencia (1/f^2). Esto significa que la potencia disminuye más rápidamente a medida que aumenta la frecuencia en comparación con el ruido rosa o rojo. El ruido Browniano se llama así por su relación con el movimiento Browniano, que es un fenómeno físico donde las partículas suspendidas en un fluido se mueven aparentemente al azar debido a las interacciones con el fluido.

En resumen, la diferencia entre estos tipos de ruido radica en la distribución de potencia en el espectro de frecuencia y la presencia o ausencia de correlación entre los valores adyacentes en el dominio del tiempo. Cada tipo de ruido tiene características distintivas que los hacen adecuados para diferentes aplicaciones y fenómenos en diversas áreas, como la música, las telecomunicaciones, la física y la ingeniería.

4.4 Ruido Rosa

En el caso del ruido rojo, la relación entre la frecuencia y la potencia es:

P = K/f^2

No hay nada especial en el exponente 2. En general, podemos generar ruido con cualquier exponente β:

P = K/f^β

Cuando β = 0, la potencia es constante en todas las frecuencias, lo que resulta en ruido blanco. Cuando β = 2, el resultado es ruido rojo. Cuando β está entre 0 y 2, el resultado se encuentra entre el ruido blanco y el ruido rojo, por lo que se llama ruido rosa.

Existen varias formas de generar ruido rosa. La más sencilla es generar ruido blanco y luego aplicar un filtro paso bajo con el exponente deseado.

thinkdsp proporciona una clase que representa una señal de ruido rosa:

class PinkNoise(_Noise): def init(self, amp=1.0, beta=1.0): self.amp = amp self.beta = beta

amp es la amplitud deseada de la señal y beta es el exponente deseado.

PinkNoise proporciona make_wave, que genera una Wave

Figura 4.7: Espectro de ruido blanco, rosa y rojo en una escala log-log.
python
def make_wave(self, duration=1, start=0, framerate=11025): signal = UncorrelatedUniformNoise() wave = signal.make_wave(duration, start, framerate) spectrum = wave.make_spectrum() spectrum.pink_filter(beta=self.beta) wave2 = spectrum.make_wave() wave2.unbias() wave2.normalize(self.amp) return wave2

duration es la duración de la onda en segundos. start es el tiempo de inicio de la onda; se incluye para que make_wave tenga la misma interfaz para todos los tipos de señales, pero para ruido aleatorio, el tiempo de inicio no es relevante. Y framerate es el número de muestras por segundo.

make_wave crea una onda de ruido blanco, calcula su espectro, aplica un filtro con el exponente deseado y luego convierte el espectro filtrado de nuevo en una onda. Luego, se realiza el proceso de desplazamiento y normalización de la onda.

Spectrum proporciona pink_filter:

python
def pink_filter(self, beta=1.0): denom = self.fs ** (beta/2.0) denom[0] = 1 self.hs /= denom

pink_filter divide cada elemento del espectro por f^(β/2). Dado que la potencia es el cuadrado de la amplitud, esta operación divide la potencia en cada componente por este factor.

Figura 4.7: Espectro de ruido blanco, rosa y rojo en una escala log-log.

python
def make_wave(self, duration=1, start=0, framerate=11025): signal = UncorrelatedUniformNoise() wave = signal.make_wave(duration, start, framerate) spectrum = wave.make_spectrum() spectrum.pink_filter(beta=self.beta) wave2 = spectrum.make_wave() wave2.unbias() wave2.normalize(self.amp) return wave2

duration es la duración de la onda en segundos. start es el tiempo de inicio de la onda; se incluye para que make_wave tenga la misma interfaz para todos los tipos de señales, pero para ruido aleatorio, el tiempo de inicio no es relevante. Y framerate es el número de muestras por segundo.

make_wave crea una onda de ruido blanco, calcula su espectro, aplica un filtro con el exponente deseado y luego convierte el espectro filtrado de nuevo en una onda. Luego, se realiza el proceso de desplazamiento y normalización de la onda.

Spectrum proporciona pink_filter:

python
def pink_filter(self, beta=1.0): denom = self.fs ** (beta/2.0) denom[0] = 1 self.hs /= denom

pink_filter divide cada elemento del espectro por f^(β/2). Dado que la potencia es el cuadrado de la amplitud, esta operación divide la potencia en cada componente por este factor.

Figura 4.8: Gráfico de probabilidad normal para las partes real e imaginaria del espectro del ruido gaussiano.

En el gráfico se muestra el espectro del ruido gaussiano, representando tanto la parte real como la imaginaria.

En la figura 4.6 se muestra la forma de onda resultante del ruido rosa. Al igual que el ruido browniano, fluctúa al azar de manera que sugiere una correlación entre los valores sucesivos, pero al menos visualmente, parece más aleatorio. En el próximo capítulo volveremos a esta observación y seré más preciso sobre lo que quiero decir con "correlación" y "más aleatorio".

Finalmente, la figura 4.7 muestra un espectro para el ruido blanco, rosa y rojo en la misma escala log-log. La relación entre el exponente β y la pendiente del espectro es evidente en esta figura.

4.5 Ruido gaussiano

Comenzamos con el ruido uniforme no correlacionado (UU) y mostramos que, debido a que su espectro tiene una potencia igual en todas las frecuencias, en promedio, el ruido UU es blanco. Sin embargo, cuando las personas hablan de "ruido blanco", no siempre se refieren al ruido UU. De hecho, más a menudo se refieren al ruido gaussiano no correlacionado (UG).

thinkdsp proporciona una implementación de ruido UG:

python
class UncorrelatedGaussianNoise(_Noise):
def evaluate(self, ts):

En este caso, el ruido se genera utilizando la distribución gaussiana. Cada valor en la señal es una muestra generada independientemente a partir de una distribución gaussiana con media cero y desviación estándar especificada.

Este tipo de ruido se utiliza con frecuencia en aplicaciones prácticas y teóricas debido a sus propiedades estadísticas bien conocidas y su relación con otros fenómenos naturales.

La función np.random.normal devuelve una matriz de NumPy con valores generados a partir de una distribución gaussiana, en este caso con media 0 y desviación estándar self.amp. En teoría, el rango de valores puede ser desde menos infinito hasta más infinito, pero se espera que aproximadamente el 99% de los valores estén entre -3 y 3.

El ruido UG (no correlacionado gaussiano) es similar en muchos aspectos al ruido UU (no correlacionado uniforme). El espectro tiene una potencia igual en todas las frecuencias, en promedio, por lo que el ruido UG también es blanco. Y tiene otra propiedad interesante: el espectro del ruido UG también es ruido UG. Más precisamente, las partes real e imaginaria del espectro son valores gaussianos no correlacionados.

Para probar esta afirmación, podemos generar el espectro del ruido UG y luego generar un "gráfico de probabilidad normal", que es una forma gráfica de probar si una distribución es gaussiana.

python
signal = thinkdsp.UncorrelatedGaussianNoise() wave = signal.make_wave(duration=0.5, framerate=11025) spectrum = wave.make_spectrum() thinkstats2.NormalProbabilityPlot(spectrum.real) thinkstats2.NormalProbabilityPlot(spectrum.imag)

NormalProbabilityPlot está proporcionado por thinkstats2, que está incluido en el repositorio de este libro. Si no estás familiarizado con los gráficos de probabilidad normal, puedes leer más al respecto en el Capítulo 5 de "Think Stats" en http://thinkstats2.com.

La Figura 4.8 muestra los resultados. Las líneas grises muestran un ajuste de modelo lineal a los datos; las líneas oscuras muestran los datos.

Una línea recta en un gráfico de probabilidad normal indica que los datos provienen de una distribución gaussiana. A excepción de algunas variaciones aleatorias en los extremos, estas líneas son rectas, lo que indica que el espectro del ruido UG es ruido UG.

El espectro del ruido UU también es ruido UG, al menos aproximadamente. De hecho, según el Teorema del Límite Central, el espectro de casi cualquier ruido no correlacionado es aproximadamente gaussiano, siempre que la distribución tenga media y desviación estándar finitas, y el número de muestras sea grande.

4.6 Ejercicios 

Las soluciones a estos ejercicios se encuentran en el archivo chap04soln.ipynb.

Ejercicio 4.1: "A Soft Murmur" es un sitio web que reproduce una mezcla de fuentes de ruido natural, como lluvia, olas, viento, etc. En http://asoftmurmur.com/about/ puedes encontrar su lista de grabaciones, la mayoría de las cuales se encuentran en http://freesound.org.

Descarga algunos de estos archivos y calcula el espectro de cada señal. ¿Se parece el espectro de potencia a ruido blanco, ruido rosa o ruido Browniano? ¿Cómo varía el espectro con el tiempo?

Ejercicio 4.2: En una señal de ruido, la mezcla de frecuencias cambia con el tiempo. A largo plazo, esperamos que la potencia en todas las frecuencias sea igual, pero en cualquier muestra, la potencia en cada frecuencia es aleatoria.

Para estimar la potencia promedio a largo plazo en cada frecuencia, podemos dividir una señal larga en segmentos, calcular el espectro de potencia para cada segmento y luego calcular el promedio entre los segmentos. Puedes obtener más información sobre este algoritmo en http://en.wikipedia.org/wiki/Bartlett’s_method.

Implementa el método de Bartlett y úsalo para estimar el espectro de potencia de una onda de ruido. Pista: revisa la implementación de make_spectrogram.

Ejercicio 4.3: En http://www.coindesk.com puedes descargar el precio diario de un BitCoin en un archivo CSV. Lee este archivo y calcula el espectro de los precios de BitCoin en función del tiempo. ¿Se asemeja a ruido blanco, ruido rosa o ruido Browniano?

Ejercicio 4.4: Un contador Geiger es un dispositivo que detecta radiación. Cuando una partícula ionizante golpea el detector, emite una ráfaga de corriente. La salida total en un momento dado se puede modelar como ruido de Poisson no correlacionado (UP, por sus siglas en inglés), donde cada muestra es una cantidad aleatoria de una distribución de Poisson, que corresponde al número de partículas detectadas durante un intervalo.

Escribe una clase llamada UncorrelatedPoissonNoise que herede de thinkdsp._Noise y proporcione el método evaluate. Debería utilizar np.random.poisson para generar valores aleatorios a partir de una distribución de Poisson. El parámetro de esta función, lam, es el número promedio de partículas durante cada intervalo. Puedes usar el atributo amp para especificar lam. Por ejemplo, si la velocidad de fotogramas es de 10 kHz y amp es 0.001, esperamos alrededor de 10 "clics" por segundo.

Genera aproximadamente un segundo de ruido UP y escúchalo. Para valores bajos de amp, como 0.001, debería sonar como un contador Geiger. Para valores más altos, debería sonar como ruido blanco. Calcula y traza el espectro de potencia para ver si se parece a ruido blanco.

Ejercicio 4.5: El algoritmo de generación de ruido rosa en este capítulo es conceptualmente sencillo pero computacionalmente costoso. Existen alternativas más eficientes, como el algoritmo de Voss-McCartney. Investiga este método, impleméntalo, calcula el espectro del resultado y confirma que tiene la relación deseada entre potencia y frecuencia.

El Capítulo 4 del libro "Think DSP" aborda el tema del ruido y su caracterización en el dominio de la señal. A continuación, se resumen algunas de las partes más relevantes del capítulo:

  1. Ruido blanco: Es un tipo de ruido en el que la potencia es constante en todas las frecuencias. Se puede generar utilizando una señal de ruido uniforme no correlacionado (UU noise).

  2. Espectro de potencia: El espectro de potencia representa cómo se distribuye la energía del ruido en diferentes frecuencias. Puede ser visualizado mediante un gráfico de log-log.

  3. Ruido Browniano: Es un tipo de ruido en el que la potencia es proporcional a 1/f^2, donde f es la frecuencia. Tiene una relación lineal en una escala log-log en el espectro de potencia. Se puede generar aplicando un filtro rosa (pink filter) al ruido blanco.

  4. Ruido rosa: Es un tipo de ruido entre el ruido blanco y el ruido Browniano. La potencia disminuye a medida que la frecuencia aumenta, pero no de manera tan pronunciada como en el ruido Browniano.

  5. Ruido rojo: Es otro término utilizado para referirse al ruido Browniano debido a que la mayor parte de su potencia se encuentra en el extremo de baja frecuencia del espectro, que se asocia con el color rojo.

  6. Generación de ruido rosa: Se puede generar ruido rosa aplicando un filtro de paso bajo con un exponente deseado al ruido blanco.

  7. Análisis del espectro de ruido: Se puede utilizar la pendiente del espectro de potencia para determinar el tipo de ruido presente. Por ejemplo, una pendiente de -2 indica ruido Browniano.

  8. Estimación del espectro de potencia a largo plazo: Para estimar la potencia promedio a largo plazo en cada frecuencia, se puede dividir una señal larga en segmentos, calcular el espectro de potencia para cada segmento y luego promediar los resultados.

El capítulo también incluye una serie de ejercicios que permiten practicar los conceptos presentados, como analizar el espectro de grabaciones de ruido natural, implementar el método de Bartlett para estimar el espectro de potencia y analizar el espectro de los precios diarios de Bitcoin.

En general, el capítulo proporciona una comprensión básica de los diferentes tipos de ruido y cómo se pueden analizar y generar en el dominio de la señal.

 5 Autocorrelacion.

En el capítulo anterior caractericé al ruido blanco como "no correlacionado", lo que significa que cada valor es independiente de los demás, y al ruido Browniano como "correlacionado", porque cada valor depende del valor anterior. En este capítulo definiré estos términos de manera más precisa y presentaré la función de autocorrelación, que es una herramienta útil para el análisis de señales.

El código de este capítulo se encuentra en el archivo "chap05.ipynb", el cual está disponible en el repositorio de este libro (consulte la Sección 0.2). También puedes verlo en el siguiente enlace: http://tinyurl.com/thinkdsp05.

5.1 Correlación

En general, la correlación entre variables significa que si conoces el valor de una, tienes cierta información sobre la otra. Hay varias formas de cuantificar la correlación, pero la más común es el coeficiente de correlación de Pearson, generalmente denotado como ρ. Para dos variables, x e y, que contienen cada una N valores, la fórmula para calcular ρ es la siguiente:


ρ =(Sumai(xi − µx)(yi − µy))/Nσxσy

Donde µx y µy son las medias de x e y, y σx y σy son sus desviaciones estándar.

La correlación de Pearson siempre está entre -1 y +1 (incluyendo ambos). Si ρ es positivo, decimos que la correlación es positiva, lo que significa que cuando una variable es alta, la otra tiende a ser alta. Si ρ es negativo, la correlación es negativa, lo que significa que cuando una variable es alta, la otra tiende a ser baja

Figura 5.1: Dos ondas senoidales que difieren en una fase de 1 radian; su coeficiente de correlación es 0.54.

La magnitud de ρ indica la fuerza de la correlación. Si ρ es 1 o -1, las variables están perfectamente correlacionadas, lo que significa que si conoces una, puedes hacer una predicción perfecta sobre la otra. Si ρ está cerca de cero, probablemente la correlación sea débil, por lo que si conoces una variable, no te dice mucho acerca de las demás.

Digo "probablemente débil" porque también es posible que exista una relación no lineal que no se capture mediante el coeficiente de correlación. Las relaciones no lineales a menudo son importantes en estadística, pero menos relevantes para el procesamiento de señales, por lo que no diré más al respecto aquí.

Python proporciona varias formas de calcular correlaciones. np.corrcoef toma cualquier número de variables y calcula una matriz de correlación que incluye las correlaciones entre cada par de variables.

Presentaré un ejemplo con solo dos variables. Primero, defino una función que construye ondas senoidales con diferentes desplazamientos de fase:

python
def make_sine(offset): signal = thinkdsp.SinSignal(freq=440, offset=offset) wave = signal.make_wave(duration=0.5, framerate=10000) return wave

A continuación, instancio dos ondas con diferentes desplazamientos:

python
wave1 = make_sine(offset=0) wave2 = make_sine(offset=1)

La Figura 5.1 muestra cómo se ven los primeros períodos de estas ondas. Cuando una onda es alta, la otra también tiende a ser alta, por lo que esperamos que estén correlacionadas.

Figura 5.2: La correlación de dos ondas senoidales en función del desplazamiento de fase entre ellas. El resultado es un coseno.

python
corr_matrix = np.corrcoef(wave1.ys, wave2.ys, ddof=0) [[ 1. 0.54] [ 0.54 1. ]]

La opción ddof=0 indica que corrcoef debe dividir por N, como en la ecuación anterior, en lugar de usar el valor predeterminado, N - 1.

El resultado es una matriz de correlación: el primer elemento es la correlación de wave1 consigo misma, que siempre es 1. De manera similar, el último elemento es la correlación de wave2 consigo misma.

Los elementos fuera de la diagonal contienen el valor que nos interesa, que es la correlación entre wave1 y wave2. El valor 0.54 indica que la fuerza de la correlación es moderada.

A medida que el desplazamiento de fase aumenta, esta correlación disminuye hasta que las ondas están desfasadas 180 grados, lo que da como resultado una correlación de -1. Luego, aumenta hasta que el desplazamiento difiere en 360 grados. En ese punto, hemos completado un ciclo completo y la correlación es 1.

La Figura 5.2 muestra la relación entre la correlación y el desplazamiento de fase para una onda senoidal. La forma de esa curva debería resultar familiar; es un coseno.

thinkdsp proporciona una interfaz sencilla para calcular la correlación entre ondas:

python
wave1.corr(wave2) 0.54

5.2 Correlación serial

Las señales a menudo representan mediciones de cantidades que varían en el tiempo. Por ejemplo, las señales de sonido con las que hemos trabajado representan mediciones de voltaje (o corriente), que corresponden a los cambios en la presión del aire que percibimos como sonido.

Mediciones como estas casi siempre tienen correlación serial, que es la correlación entre cada elemento y el siguiente (o el anterior). Para calcular la correlación serial, podemos desplazar una señal y luego calcular la correlación de la versión desplazada con la original.

python
def serial_corr(wave, lag=1):
 n = len(wave) 
 y1 = wave.ys[lag:] 
 y2 = wave.ys[:n-lag] 
 corr = np.corrcoef(y1, y2, ddof=0)[0, 1
return corr

serial_corr toma un objeto Wave y lag, que es el número entero de lugares para desplazar la onda. Calcula la correlación de la onda con una versión desplazada de sí misma.

Podemos probar esta función con las señales de ruido del capítulo anterior. Esperamos que el ruido UU no tenga correlación, según la forma en que se genera (sin mencionar el nombre):

python
signal = thinkdsp.UncorrelatedGaussianNoise() 
wave = signal.make_wave(duration=0.5, framerate=11025) serial_corr(wave)

Cuando ejecuté este ejemplo, obtuve 0.006, lo que indica una correlación serial muy pequeña. Es posible que obtengas un valor diferente cuando lo ejecutes, pero debería ser comparativamente pequeño.

En una señal de ruido Browniano, cada valor es la suma del valor anterior y un "paso" aleatorio, por lo que esperamos una fuerte correlación serial:

python
signal = thinkdsp.BrownianNoise()
 wave = signal.make_wave(duration=0.5, framerate=11025) serial_corr(wave)

Efectivamente, obtuve un resultado mayor que 0.999.

Dado que el ruido rosa está en cierto sentido entre el ruido Browniano y el ruido UU, podríamos esperar una correlación intermedia:

python
signal = thinkdsp.PinkNoise(beta=1)
 wave = signal.make_wave(duration=0.5, framerate=11025) serial_corr(wave)

Figure 5.3: Correlación serial para ruido rosa con diferentes parámetros.

Con el parámetro β = 1, obtuve una correlación serial de 0.851. A medida que variamos el parámetro desde β = 0, que es ruido no correlacionado, hasta β = 2, que es ruido Browniano, la correlación serial varía desde 0 hasta casi 1, como se muestra en la Figura 5.3.

5.3 Autocorrelación

En la sección anterior calculamos la correlación entre cada valor y el siguiente, por lo que desplazamos los elementos del arreglo en 1. Pero también podemos calcular fácilmente correlaciones seriales con diferentes desplazamientos.

Podemos pensar en serial_corr como una función que mapea cada valor de desplazamiento a la correlación correspondiente, y podemos evaluar esa función iterando a través de los valores de desplazamiento:

python
def autocorr(wave): 
 lags = range(len(wave.ys)//2
 corrs = [serial_corr(wave, lag) for lag in lags] 
return lags, corrs

autocorr toma un objeto Wave y devuelve la función de autocorrelación como un par de secuencias: lags es una secuencia de enteros desde 0 hasta la mitad de la longitud de la onda; corrs es la secuencia de correlaciones seriales para cada desplazamiento.

La Figura 5.4 muestra las funciones de autocorrelación para ruido rosa con tres valores de β. Para valores bajos de β, la señal tiene una menor correlación y la función de autocorrelación cae rápidamente a cero. Para valores más altos, la correlación serial es más fuerte y disminuye más lentamente. Con β = 1.7, la correlación serial es fuerte.



Figure 5.4: Funciones de autocorrelación para ruido rosa con diferentes parámetros.

Incluso para desplazamientos largos, las funciones de autocorrelación del ruido rosa muestran una correlación significativa; este fenómeno se conoce como dependencia a largo plazo, porque indica que cada valor en la señal depende de muchos valores anteriores.

5.4 Autocorrelación de señales periódicas

La autocorrelación del ruido rosa tiene propiedades matemáticas interesantes, pero aplicaciones limitadas. La autocorrelación de señales periódicas es más útil.

Como ejemplo, descargué de freesound.org una grabación de alguien cantando un trino; el repositorio de este libro incluye el archivo: 28042_bcjordan__voicedownbew.wav. Puedes usar el cuaderno Jupyter de este capítulo, chap05.ipynb, para reproducirlo.

La Figura 5.5 muestra el espectrograma de esta onda. La frecuencia fundamental y algunos de los armónicos se muestran claramente. El trino comienza cerca de 500 Hz y desciende hasta aproximadamente 300 Hz, aproximadamente desde C5 hasta E4.

Para estimar la altura en un punto particular en el tiempo, podríamos usar el espectro, pero no funciona muy bien. Para ver por qué no, tomaré un segmento corto de la onda y trazaré su espectro:

python
duration = 0.01 segment = wave.segment(start=0.2, duration=duration) spectrum = segment.make_spectrum() spectrum.plot(high=1000)

Figura 5.5: Espectrograma de un chirrido vocal

Figura 5.6: Espectro de un segmento de un chirrido vocal

Figure 5.7: Dos segmentos de un trino, uno comenzando 0.0023 segundos después del otro.

Este segmento comienza en 0.2 segundos y dura 0.01 segundos. La Figura 5.6 muestra su espectro. Hay un pico claro cerca de 400 Hz, pero es difícil identificar la altura con precisión. La longitud del segmento es de 441 muestras con una frecuencia de muestreo de 44100 Hz, por lo que la resolución de frecuencia es de 100 Hz (ver Sección 3.5). Esto significa que la estimación de la altura podría tener un margen de error de 50 Hz; en términos musicales, el rango de 350 Hz a 450 Hz equivale aproximadamente a 5 semitonos, ¡lo cual es una diferencia considerable!

Podríamos obtener una mejor resolución de frecuencia tomando un segmento más largo, pero dado que la altura está cambiando con el tiempo, también obtendríamos "desenfoque de movimiento"; es decir, el pico se dispersaría entre la altura de inicio y la altura final del segmento, como vimos en la Sección 3.3.

Podemos estimar la altura con mayor precisión utilizando la autocorrelación. Si una señal es periódica, esperamos que la autocorrelación se dispare cuando el desfase sea igual al periodo. Para mostrar por qué esto funciona, trazaré dos segmentos de la misma grabación.

python
import matplotlib.pyplot as plt def plot_shifted(wave, offset=0.001, start=0.2): segment1 = wave.segment(start=start, duration=0.01) segment1.plot(linewidth=2, alpha=0.8) segment2 = wave.segment(start=start-offset, duration=0.01) segment2.shift(offset) segment2.plot(linewidth=2, alpha=0.4)

Figure 5.8: Función de autocorrelación para un segmento de un trino.
python
corr = segment1.corr(segment2) text = r'$\rho =$ %.2g' % corr plt.text(segment1.start+0.0005, -0.8, text) plt.xlabel('Tiempo (s)')

Una de las segmentos comienza en 0.2 segundos; el otro comienza 0.0023 segundos después. La Figura 5.7 muestra el resultado. Los segmentos son similares y su correlación es de 0.99.

Este resultado sugiere que el período está cerca de 0.0023 segundos, lo que corresponde a una frecuencia de 435 Hz.

Para este ejemplo, estimé el período mediante prueba y error. Para automatizar el proceso, podemos usar la función de autocorrelación.

python
lags, corrs = autocorr(segment) 
plt.plot(lags, corrs)

La Figura 5.8 muestra la función de autocorrelación para el segmento que comienza en t = 0.2 segundos. El primer pico ocurre en el desfase 101. Podemos calcular la frecuencia que corresponde a ese período de la siguiente manera:

python
period = lag / segment.framerate
 frequency = 1 / period

La frecuencia fundamental estimada es de 437 Hz. Para evaluar la precisión de la estimación, podemos realizar el mismo cálculo con los desfases 100 y 102, que corresponden a las frecuencias 432 y 441 Hz. La precisión de frecuencia utilizando la autocorrelación es inferior a 10 Hz, en comparación con los 100 Hz utilizando el espectro. En términos musicales, el error esperado es de aproximadamente 30 centésimas (un tercio de un semitono).


En el contexto del procesamiento de señales, la definición de correlación a menudo se simplifica. Para señales no sesgadas con una media de 0 y señales normalizadas con una desviación estándar de 1, el coeficiente de correlación ρ se puede simplificar a:

ρ = 1/N * Σ(xi * yi)

Y se simplifica aún más a:

r = Σ(xi * yi)

Esta definición de correlación no está estandarizada, por lo que no necesariamente se encuentra entre -1 y 1. Sin embargo, tiene otras propiedades útiles.

En términos de vectores, esta fórmula se puede reconocer como el producto escalar de dos vectores, x · y. El producto escalar indica el grado de similitud entre las señales. Si las señales están normalizadas de manera que sus desviaciones estándar sean 1, entonces el producto escalar, x · y, es igual al coseno del ángulo θ entre los vectores:

x · y = cos θ

Esta relación explica por qué la Figura 5.2 muestra una curva coseno.

5.6 Usando Numpy

N.umPy proporciona una función llamada correlate que se puede utilizar para calcular la correlación o autocorrelación de señales. Podemos usarla para calcular la autocorrelación del segmento de la sección anterior de la siguiente manera:

python
corrs2 = np.correlate(segment.ys, segment.ys, mode='same')

El arreglo resultante corrs2 contiene los valores de autocorrelación para diferentes desplazamientos (lags).

Figura 5.9: Función de autocorrelación calculada con np.correlate.

La opción 'mode' le indica a 'correlate' qué rango de desplazamiento (lag) utilizar. Con el valor 'same', el rango va desde -N/2 hasta N/2, donde N es la longitud del arreglo de la señal.

La Figura 5.9 muestra el resultado. Es simétrica debido a que las dos señales son idénticas, por lo que un desplazamiento negativo en una tiene el mismo efecto que un desplazamiento positivo en la otra.

Para comparar con los resultados obtenidos con 'autocorr', podemos seleccionar la segunda mitad de los valores:

N = len(corrs2) half = corrs2[N//2:]

Si comparas la Figura 5.9 con la Figura 5.8, notarás que las correlaciones calculadas por np.correlate disminuyen a medida que aumentan los desplazamientos. Esto se debe a que np.correlate utiliza la definición de correlación no estandarizada; a medida que el desplazamiento aumenta, la superposición entre las dos señales disminuye, por lo que la magnitud de las correlaciones también disminuye.

Podemos corregir esto dividiendo los valores por las longitudes correspondientes:

lengths = range(N, N//2, -1) half /= lengths

Finalmente, podemos estandarizar los resultados para que la correlación con desplazamiento 0 sea igual a 1:

half /= half[0]

Con estos ajustes, los resultados obtenidos con 'autocorr' y np.correlate son casi idénticos. Aún pueden diferir en un 1-2%. La razón no es importante, pero si tienes curiosidad: 'autocorr' estandariza las correlaciones de manera independiente para cada desplazamiento, mientras que en np.correlate las estandarizamos al final.

Lo más importante es que ahora sabes qué es la autocorrelación, cómo utilizarla para estimar el período fundamental de una señal y dos formas de calcularla.


Ejercicio 5.1: Para estimar el tono del trino vocal para diferentes tiempos de inicio, puedes seguir estos pasos: Abre el cuaderno Jupyter chap05.ipynb. Busca la sección que te permite calcular autocorrelaciones para diferentes retardos. Utiliza la interacción proporcionada para ajustar el tiempo de inicio del trino vocal. Ejecuta el código para calcular la autocorrelación y observa los resultados. Repite los pasos 3 y 4 para diferentes tiempos de inicio y estimar el tono del trino vocal. Ejercicio 5.2: Para encapsular el código en una función llamada estimar_fundamental y realizar un seguimiento del tono de un sonido grabado, puedes seguir estos pasos: Define una función llamada estimar_fundamental que tome como entrada el sonido grabado. Dentro de la función, implementa el código del ejemplo en chap05.ipynb que utiliza la autocorrelación para estimar la frecuencia fundamental del sonido. Devuelve la frecuencia fundamental estimada. Utiliza la función estimar_fundamental para realizar un seguimiento del tono del sonido grabado en diferentes intervalos de tiempo. Superpone las estimaciones de tono en un espectrograma del sonido grabado para visualizar los resultados. Ejercicio 5.3: Para calcular la autocorrelación de los precios de Bitcoin utilizando los datos históricos de precios, sigue estos pasos: Obtén los datos históricos de precios de Bitcoin. Calcula la autocorrelación de los precios de Bitcoin utilizando la función de autocorrelación. Grafica la función de autocorrelación y observa si disminuye rápidamente. Analiza los resultados de la autocorrelación para determinar si hay evidencia de comportamiento periódico en los precios de Bitcoin. Ejercicio 5.4: Para completar el Ejercicio 5.4, sigue estos pasos: Localiza el cuaderno Jupyter saxophone.ipynb en el repositorio de este libro. Abre el cuaderno y lee su contenido. Ejecuta los ejemplos proporcionados en el cuaderno para explorar la autocorrelación, la percepción del tono y el fenómeno del tono fundamental ausente. Selecciona un segmento diferente de la grabación en el cuaderno y vuelve a ejecutar los ejemplos para observar cualquier diferencia en los resultados. Mira el video de Vi Hart "What is up with Noises? (The Science and Mathematics of Sound, Frequency, and Pitch)" para obtener una comprensión más profunda del fenómeno del tono fundamental ausente y la percepción del tono.

Resumen del Capítulo 5:

El Capítulo 5 del libro aborda el tema de la autocorrelación y su aplicación en el análisis de señales. Se explora cómo calcular la autocorrelación de una señal y cómo interpretar los resultados obtenidos. Además, se comparan diferentes métodos para estimar la frecuencia fundamental de una señal periódica, como el uso de espectrogramas y la autocorrelación.

El capítulo comienza explicando la autocorrelación y su relación con la correlación de Pearson. Se muestran diferentes ejemplos y gráficas para ilustrar cómo la autocorrelación puede revelar información sobre la estructura y periodicidad de una señal.

Se presenta el concepto de autocorrelación serial y se muestra cómo calcular la función de autocorrelación para diferentes retardos. Se discute la noción de correlación a partir del producto punto entre dos señales y cómo se relaciona con la similitud entre las señales.

A continuación, se explora el uso de la autocorrelación para estimar la frecuencia fundamental de una señal periódica. Se muestra cómo calcular la autocorrelación de segmentos de la señal y encontrar el pico correspondiente al periodo fundamental. Se compara esta técnica con el análisis espectral utilizando espectrogramas.

El capítulo concluye introduciendo la función np.correlate de la biblioteca NumPy, que proporciona una forma conveniente de calcular la autocorrelación de una señal.

Puntos importantes:

  1. La autocorrelación es una medida de la similitud entre una señal y una versión desplazada de sí misma.
  2. La autocorrelación puede revelar patrones y periodicidad en una señal.
  3. La autocorrelación se puede utilizar para estimar la frecuencia fundamental de una señal periódica.
  4. La autocorrelación puede calcularse utilizando el producto punto entre dos señales.
  5. La función np.correlate de NumPy proporciona una forma eficiente de calcular la autocorrelación.
  6. La autocorrelación puede complementar el análisis espectral en el estudio de señales.

Es importante tener en cuenta que la autocorrelación es una herramienta poderosa en el análisis de señales y puede tener aplicaciones en diversos campos, como procesamiento de señales, análisis de audio, reconocimiento de patrones y más. Comprender cómo calcular y interpretar la autocorrelación puede ampliar nuestras capacidades para analizar y comprender las señales en nuestro entorno.

No hay comentarios.:

Publicar un comentario