1. Experimento, espacio muestral, evento

Referencia: Leon-García 2.1 p21, Gubner 1.1 p.7

Experimento

Un experimento aleatorio es un experimento en que el resultado varia de forma no predecible cuando se repite bajo las mismas condiciones.

Un experimento aleatorio se especifica dado el procedimiento experimental y el grupo de uno o mas medidas y observaciones.

  • Un solo procedimiento, pero con diferentes en las observaciones
  • puede tener mas de una medida u observación
  • las medidas pueden ser valores contínuos
  • pueden ser secuenciales, o una secuencia de experimentos más simples o subexperimentos

Espacio Muestral

Se define como resultado o muestra ζ (zeta griega) de un experimento aleatorio como el resultado que no puede ser descompuesto a partir de otros resultados, o una combinación de ellos.

El espacio muestral S de un experimento aleatorio se define como el grupo de todos los resultados posibles, que se especifican como:

  • Lista de todos los resultados a={0,1,2,3}
  • Mediante una propiedad de los resultados A={x: 0≤x≤3, x es entero}

Un espacio muestral discreto es contable, ejemplo las caras de un dado.
un espacio muestral contínuo es contable infinito o no contable, ejemplo los valores de voltaje de descarga de un capacitor, o el tiempo hasta la primera falla de un equipo.

Algunas veces es conveniente incluir en el espacio muestral el resultado imposible. Ejemplo, incluir que el espacio muestral es un número real positivo para el tiempo antes de la falla de un dispositivo, aunque conocemos que no tiene tiempo de operación infinito.

Eventos

Normalmente no es de interés un resultado, muestra u observación específica, el interés se enfoca más en que ocurra un evento, lo que requiere usar subgrupos de S. Por ejemplo, que el voltaje de una señal sea negativo {señal: -∞ < señal < 0}

Dos eventos de interés especial son: evento de certeza S, y el evento imposible o nulo ∅, que no resultados o nunca ocurre.

Eventos Clase

La teoría de probabilidades usa también lo que se conoce como una clase F con los eventos de interés, que son aquellos a los que se les asigna valores de probabilidad.

Axiomas de probabilidad

Referencia: Leon-García 1.3.3 p.23, Gubner 1.4 p.36, Parsen p.9/pdf.30, Ross p.4

La teoría de probabilidad inicia con un grupo de axiomas que indican las propiedades que los valores de probabilidad deben satisfacer.

Supone que se definieron:

  1. el experimento aleatorio con todos los posibles valores en el espacio muestral S
  2. las clases y subgrupos llamados eventos
  3. cada evento A tiene asignado un valor P[A]

con lo que se satisfacen los AXIOMAS:

  • Axioma I    :  0 ≤ P[A] ≤ 1
  • Axioma II   : P[S] = 1
  • Axioma III  : Si A∩B=∅,  P[A∪B] = P[A] + P[B]

Lo que se expresa como, que la probabilidad siempre es no negativa, que la probabilidad del espacio muestral es 1, que incluye la probabilidad P[∅] que a veces se le denomina «imposible» e igual a 0.

Resultados de los axiomas

Unión de eventos finitos y separados

P\left[ \bigcup\limits_{n=1}^{N} A_n \right] = \sum\limits_{n=1}^{N} P[A_n]

Probabilidad de un complemento

P[A^c] = 1 - P[A]

Evento imposible

P[\varnothing] = 0

Monotonicidad

A \subset B ,\text{ implica } P[A] \leq P[B]

Inclusión – Exclusión

P[A \cup B] = P[A] + P[B] - P[A \cap B]

Propiedades de límites

P \left[ \bigcup\limits_{n=1}^{\infty} A_n \right] = \lim \limits_{N \rightarrow \infty} P\left[ \bigcup\limits_{n=1}^{N} A_n \right] P \left[ \bigcap\limits_{n=1}^{\infty} A_n \right] = \lim \limits_{N \rightarrow \infty} P\left[ \bigcap\limits_{n=1}^{N} A_n \right]

Union bound/countable subadditivity

P \left[ \bigcup\limits_{n=1}^{\infty} A_n \right] \leq \sum \limits_{n=1}^{\infty} P[ A_n]

2. Regularidad estadística y Frecuencia Relativa

Para que un modelo sea útil, se debe poder realizar predicciones del comportamiento del sistema en el futuro, es decir debe tener un comportamiento regular.

Muchos modelos de probabilidad en ingeniería se fundamentan en los promedios obtenidos en las secuencias de resultados, así las muchas repeticiones o intentos del experimento aleatorio son consistentes y tienden aproximadamente al mismo valor.

Esta es la propiedad de regularidad estadítica.

En el experimento de la urna y pelotas con S={0, 1, 2}, se pueden procesar las muestras del archivo: aleatoriosurna100.txt

muestras = [1,1,1,1,1,2,1,2,2,2,2,0,0,0,0,0,2,1,2,2,
            0,2,0,2,1,2,1,2,2,0,2,0,1,2,0,2,2,0,2,1,
            1,1,1,1,0,2,1,2,1,0,1,0,1,2,0,2,2,2,1,0,
            2,2,0,2,0,2,1,1,0,0,2,1,0,1,0,0,2,1,2,0,
            1,1,2,0,0,0,1,0,2,2,2,0,0,0,1,2,1,1,0,0]

Algoritmo en Python

Para seleccionar los valores únicos de las muestras y contar el número de veces de cada uno se puede usar la instrucción de numpy np.unique

los valores unicos:   [ 0.  1.  2.]
la cuenta de unicos:  [33 31 36]
Nk:
[[  0.  33.]
 [  1.  31.]
 [  2.  36.]]

Instrucciones en Python

# Muestras de una urna con pelotas S=[0,1,2]
# propuesta: edelros@espol.edu.ec
import numpy as np

# INGRESO
# archivotxt=input('nombre del archivo.txt: ')
archivotxt = 'aleatoriosurna100.txt'

# PROCEDIMIENTO
muestras = np.loadtxt(archivotxt)
n = len(muestras)

# cuenta los únicos
[unicos,cuenta] = np.unique(muestras,return_counts=True)

# une en tabla N, y aplica T transpuesta
N = np.asarray((unicos, cuenta)).T 

# SALIDA 
print('los eventos unicos:  ', unicos)
print('la cuenta de unicos: ', cuenta)
print('Nk:')
print(N)

El número de veces que aparece la pelota 0 en el experimento es N0(n)=33,la pelota N1(n)=31 y N2(n)=36.


Frecuencia relativa

La frecuencia relativa de un resultado k, se define como:

f_k(n)=\frac {N_k(n)}{n}

y lo que la regularidad estadística quiere decir es que fk(n) varía siempre acercandose a una constante pk cuando n tiende a ser muy grande.

lim_{n\to\infty} f_k(n) = p_k
# frecuencia relativa
fk  = cuenta/n
fkn = np.asarray((unicos, fk)).T
print(fkn)
[[ 0.    0.33]
 [ 1.    0.31]
 [ 2.    0.36]]

Tarea

Repita el experimento de la urna con n=1000 intentos, y determine las frecuencias relativas de cada pelota. Use el algoritmo de la hoja «modelo de probabilidad» para generar el archivo de aleatorios


Ejercicios

Referencia: Gubner 1.1 p.18

Se lanza una moneda 100 veces, y se anota la frecuencia relativa de las caras (0) y sellos (1). La experiencia nos indica que la frecuencia relativa debería estar alrededor de 0.5. Realice el experimento simulado y compare los resultados.

 caras de moneda: 
[0 0 1 1 1 0 0 1 0 1 0 0 0 0 0 0 0 0 1 0 
 1 0 1 1 1 0 0 1 1 0 1 1 1 0 1 0 0 1 1 1 
 0 0 1 0 1 0 1 0 0 1 0 0 1 0 0 0 0 1 0 1
 1 1 1 0 1 0 0 0 0 0 0 0 1 0 1 0 0 1 1 1
 0 1 1 0 0 1 0 0 0 1 0 0 1 1 1 0 0 1 1 1]
frecuencia relativa:
[[ 0.    0.55]
 [ 1.    0.45]]

# Simulador de lanza monedas [0,1]
# propuesta: edelros@espol.edu.ec
import numpy as np
import matplotlib.pyplot as plt

# INGRESO 
intentos = 100
eventos  = 2

# PROCEDIMIENTO
# lanza moneda
caramoneda = np.random.randint(eventos,size=intentos)
# cuenta los únicos
[unicos, cuenta] = np.unique(caramoneda,return_counts=True)

# frecuencia relativa
fk  = cuenta/intentos

# une en tabla fkn, y aplica T transpuesta
fkn = np.asarray((unicos, fk)).T 

# SALIDA
print(' caras de moneda: ')
print(caramoneda)
print('frecuencia relativa:')
print(fkn)

# grafica pmf
plt.plot(caramoneda,'o--')
plt.xlabel('intentos')
plt.ylabel('caramoneda')
plt.margins(0.1)
plt.show()

Ejercicio

Referencia: Gubner 1.2 p.18, Leon-García pdf/p.21

Un experimento por si mismo es lanzar una moneda 100 veces y guardar las frecuencias relativas de cada caras.

Dado que el numero de caras puede estar en el rango entre 0 y 100, hay 101 posibilidades de respuesta.

Realice el experimento anterior 1000 veces y almacene la frecuencia relativa para cada experimento. Muestre los resultados.

# Simulador de lanza monedas [0,1]
# Ejercicio: Gubner 1.2 p.18, Leon-García pdf/p.21
# propuesta: edelros@espol.edu.ec
import numpy as np
import matplotlib.pyplot as plt

# INGRESO 
experimentos = 1000
intentos = 100
eventos  = 2

# PROCEDIMIENTO
# Para cada experimento
resultado = np.zeros(shape=(experimentos,eventos),dtype=int)
for i in range(0,experimentos):
    
    # lanza moneda
    caramoneda = np.random.randint(eventos,size=intentos)
    n=len(caramoneda)
    [unicos, cuenta] = np.unique(caramoneda,return_counts=True)
    
    # resultado de cada experimento
    for j in range(0,eventos):
        resultado[i][j] = cuenta[j]

# Frecuencias relativas por evento
s = np.arange(eventos)
fr_cuenta = []
fr_unicos = []
for evento in s:
    unevento = resultado[:,evento]
    [unicos,cuenta] = np.unique(unevento,return_counts=True)
    fr_cuenta.append(cuenta)
    fr_unicos.append(unicos)

# SALIDA - gráfica
for evento in s:
    plt.plot(fr_unicos[evento]/100,fr_cuenta[evento],label=str(evento))
plt.xlabel('frecuencia relativa')
plt.legend()
plt.show()

Tarea

Escriba sus observaciones en el ejercicio anterior, cambiando el número de intentos en el rango de [10, 100000].

Referencia: Leon-García 1.3 p.5, Gubner p.3

4. Modelos de Probabilidad

Referencia: Leon-García 1.3 p.4

Muchos modelos de señales muestran variaciones impredecibles y un comportamiento aleatorio o probabilístico.

En un experimento aleatorio, repetido sobre las mismas condiciones, entrega resultados que varian de forma poco predecible. Si consideramos realizar un modelos determinístico, encontramos que no es el más apropiado al predecir siempre el mismo resultado.

Experimento aleatorio clásico

http://www.gifsanimados.org/data/media/994/bingo-imagen-animada-0002.gif

En un ánfora se colocan tres pelotas de igual tamaño numeradas con 0, 1 y 2, se agita el ánfora para que no exista un «orden» predeterminado.

Al sacar una pelota cualquiera, no se podría predecir cual será, luego se anota el número de la pelota y se la devuelve al ánfora.

El resultado del experimento, cada vez que se realiza un intento, es un número que pertenece al grupo S={0,1,2}, que formalmente se lo denomina espacio muestral

El resultado del experimento para 20 intentos, se obtiene un vector con los valores siguientes:

modeloprob01.png

Simulación en Python

Para simular el experimento con python, se puede usar el generador de enteros aleatorios np.random.randint().

# Simulador de Urna con pelotas 0,1,2 
# propuesta: edelros@espol.edu.ec
import numpy as np
import matplotlib.pyplot as plt

# INGRESO 
intentos = 20
pelotas  = 3

# PROCEDIMIENTO
aleatorios = np.random.randint(pelotas,size=intentos)

# SALIDA
print(aleatorios)

# grafica
plt.plot(aleatorios,'o-')
plt.xlabel('intentos')
plt.ylabel('aleatorios')
plt.margins(0.1)
plt.show()

# SALIDA - ARCHIVO
archivotxt = 'aleatoriosurna.txt'
np.savetxt(archivo, aleatorios, fmt='%d ')
print(' ... archivo: '+archivotxt+' guardado ...')

Revise la aplicación codificar a morse

Ahora, revise el resultado obtenido al codificar un mensaje cualquiera en morse, y realice sus observaciones:

codigo morse :
'. ... .--. --- .-.. .. -- .--. ..- .-.. ... .- 
 -. -.. --- .-.. .- ... --- -.-. .. . -.. .- -..
 -.. . .-.. -.-. --- -. --- -.-. .. -- .. . -. -
 --- ' 
  1. ¿Cuál sería el conjunto posible de resultados?
  2. defina el espacio muestral S

Señal Componentes

Si diseña un equipo para transmitir el tono «la» de 440Hz, debe tomar en cuenta que no siempre es posible obtener una expresión matemática «sencilla» de una señal.

Escuche los siguientes sonidos, vea sus gráficas y realice sus observaciones:

440Hz_44100Hz_16bit_05sec.wav

440Hz_piano.wav

440Hz_violin_A4.wav

senal440Hz_pianoviolin

        • Proponga algunas expresiones matemáticas para cada una de las graficas mostradas.
        • ¿Que tienen en común cada una de las señales?
        • ¿Que tienen de diferente?

Algoritmo en Python

Para obtener la gráfica de los sonidos, se leen los datos y frecuencias de muestreo de cada uno de los archivos.wav

# Señales analógicas
import numpy as np
import matplotlib.pyplot as plt
import scipy.io.wavfile as waves

# INGRESO
archivo1 = '440Hz_44100Hz_16bit_05sec.wav'
[muestreo1, sonido1] = waves.read(archivo1)

archivo2 = '440Hz_piano.wav'
[muestreo2, sonido2] = waves.read(archivo2)

archivo3 = '440Hz_violin_A4.wav'
[muestreo3, sonido3] = waves.read(archivo3)

# SALIDA - Observacion intermedia
print('archivo, frecuencia , dimensiones:')
print('sonido1: ', muestreo1, np.shape(sonido1))
print('sonido2: ', muestreo2, np.shape(sonido2))
print('sonido3: ', muestreo3, np.shape(sonido3))

con lo que se obtiene:

archivo, frecuencia , dimensiones:
sonido1:  44100 (220500,)
sonido2:  11025 (16538,)
sonido3:  22050 (98160,)

Nota: si el archivo.wav tiene etiquetas, se descartan para el ejercicio. por lo que se presenta una advertencia «WavFileWarning».

c:\python34\lib\site-packages\scipy\io\wavfile.py:179: WavFileWarning: Chunk (non-data) not understood, skipping it.
  WavFileWarning)

Se observa que:

  • todos los «sonidos» son monofónicos,
  • las frecuencias de muestreo de cada archivo son diferentes
  • las dimensiones de cada «sonido» o cantidad de muestras indican una duración diferente

Se recomienda escuchar nuevamente cada sonido y decidir el intervalo de observación para las gráficas. Se propone usar una ventana de observación o fragmento desde de los 1.2 segundos por 11 milisegundos.

Se extraen los fragmentos de cada sonido y se determinan los tiempos a los que corresponden basados en cada dt, con lo que se grafican los fragmentos.

# PROCEDIMIENTO
# ventana de observación
inicia  = 1.2
termina = 1.211
canal   = 0

# Extraen los fragmentos de sonido
dt1 = 1/muestreo1
t1  = np.arange(inicia,termina,dt1)
muestras   = len(t1)
fragmento1 = sonido1[int(inicia/dt1):int(inicia/dt1)+muestras]

dt2 = 1/muestreo2
t2  = np.arange(inicia,termina,dt2)
muestras   = len(t2)
fragmento2 = sonido2[int(inicia/dt2):int(inicia/dt2)+muestras]

dt3 = 1/muestreo3
t3  = np.arange(inicia,termina,dt3)
muestras  = len(t3)
fragmento3 = sonido3[int(inicia/dt3):int(inicia/dt3)+muestras]

# SALIDA
# grafica
plt.subplot(311)
plt.plot(t1,fragmento1)
plt.ylabel('sonido1(t)')

plt.subplot(312)
plt.plot(t2,fragmento2)
plt.ylabel('sonido2(t)')

plt.subplot(313)
plt.plot(t3,fragmento3)
plt.ylabel('sonido3(t)')
plt.xlabel('t segundos')

plt.show()

senal440Hz_pianoviolin
Una vez obtenidas las gráficas, realice observaciones adicionales a cada una de las señales de sonido:
1.
2.
3.

Pregunta:
Para la nota musical «la», ejecutada en un violin tradicional tradicional de madera, ¿tendrá la misma forma de señal siempre?.

Considere las siguientes situaciones:
a. Tocada por diferentes músicos
b. en diferentes lugares
c. en diferentes días


Tarea

Buscar la nota «La» en otros instrumentos en formato.wav. De ser necesario realizar observaciones adicionales.

Referencia: Lathi 1.3.2 pdf/p.63, Estándar ISO16, Frecuencias de afinación de un piano

Señal Determinística

Referencia: Lathi 1.3.2 pdf/p.63, Estándar ISO16, Frecuencias de afinación de un piano

Señal determinística, conocida como una señal cuya descripción física es completamente conocida por su forma matemática o gráfica.

Ejemplo: Señal 440 Hz

Es una señal que puede ser escrita de forma determinística como:

x(t)= cos(2\pi f t)

Siendo f=440 Hz

La señal de 440 Hz es usada como la frecuencia referencia para la afinación de todos los instrumentos musicales desde 1936, adoptada por ISO en 1955 y reafirmado en 1975 como ISO16.

señal 440Hz

Algoritmo en Python – Graficar en dominio del tiempo

El audio del archivo.wav se puede escuchar con windows media player. Para visualizar la forma de señal del archivo.wav, se puede usar python para leer el archivo y graficar la señal en el dominio del tiempo.

El archivo.wav debe encontrarse en el mismo directorio que el archivo .py de python.

# Señales analógicas
import numpy as np
import matplotlib.pyplot as plt
import scipy.io.wavfile as waves

# INGRESO
# archivo = input('archivo.wav a leer: ')
archivo = '440Hz_44100Hz_16bit_05sec.wav'
[muestreo, sonido] = waves.read(archivo)

# SALIDA - Observación intermedia
print('frecuencia de muestreo: ', muestreo)
print('dimensiones de matriz: ', np.shape(sonido))
frecuencia de muestreo:  44100
dimensiones de matriz:  (220500,)

Nota: si el archivo.wav tiene etiquetas, se descartan para el ejercicio. por lo que se presenta una advertencia «WavFileWarning».

c:\python34\lib\site-packages\scipy\io\wavfile.py:179: WavFileWarning: Chunk (non-data) not understood, skipping it. WavFileWarning)

Existen 220500 muestras en el archivo del ejemplo, que en una sola gráfica puede resultar muy denso y poco observable. La gráfica se limitará a una ventana de tiempo del archivo definida en el rango de tiempo por inicio y termina, con valores en el rango de milisegundos.

El rango de la ventana puede cambiarse si es de interés usar otro rango de tiempo.

De las dimensiones de la matriz (220500,), se determina que el audio es de tipo monofónico al tener una sola dimensión. El audio estéreo tiene dimensiones de (n,2).

# PROCEDIMIENTO
# ventana de observación
inicia  = 0
termina = 0.005
canal   = 0

# selecciona datos
dt = 1/muestreo
t  = np.arange(inicia,termina,dt)
muestras = len(t)
fragmento = sonido[int(inicia/dt):int(inicia/dt)+muestras]

# SALIDA
plt.plot(t,fragmento)
plt.ylabel('sonido(t)')
plt.xlabel('t segundos')
plt.show()

señal 440Hz

Aunque la señal es simple de ver, y oir, la riqueza del sonido de un instrumento corresponde a otros componentes de frecuencia que los hace únicos al oido estando en la misma nota musical.

Sigma-Delta Decodificador con Python

Referencia: Leon-Couch, 3–8 Modulación Delta, p.192 ; Delta-sigma_modulation Wikipedia, Sigma-Delta – Modulación

La señal de entrada para el decodificador es la señal codificada en los archivo.txt del ejemplo anterior:

deltasigma_datos.txt, deltasigma_parametros.txt

Los archivos para ser leidos deben copiarse al directorio donde se encuentra el algoritmo.py.

Los parametros obtenidos del archivo.txt son:

  • ΔY = deltaY = datos[1]
  • Δt = datos[0]
  • muestras en el rango de observación: k como tamaño del arreglo yentrada.

DeltaSigma_Decodificador

% matplotlib inline
# Modulacion Sigma-Delta Decodificador
#  entrada y[n], salida: x[n]
# propuesta:edelros@espol.edu.ec
import numpy as np
import matplotlib.pyplot as plt
import scipy.io.wavfile as waves

# INGRESO - Lectura de datos desde archivos
# archivoparam = intput('archivo.txt de parametros: ')
# archivodatos = input('archivo.txt de datos: ')
# arhivoaudio = input('archivo.wav a grabar: ')

archivoparam = 'deltasigma_parametros.txt'
archivodatos = 'deltasigma_datos.txt'
archivoaudio = 'sigmadeltaaudio.wav'

param = np.loadtxt(archivoparam,dtype=float)
yentrada = np.loadtxt(archivodatos,dtype=int)

Decodificar la señal +1 y -1 en el vector yentrada consiste en acumular los valores de la secuencia en xdigital usando escalones de tamaño deltaY.

Se procede de forma semejante para los tiempos td usados para el eje de las abscisas usando pasos deltaT.

# PROCEDIMIENTO
deltaT = param[0] # Tamaño delta en eje tiempo
deltaY = param[1] # Tamaño delta en eje Y
k = len(yentrada) # número de muestras
xdigital = np.zeros(k, dtype='int16')
punto = np.zeros(k, dtype=int)      # número de muestras
td = np.zeros(k, dtype=float)    # tiempo muestreo

# DECOdifica Sigma-Delta
xdigital[0] = yentrada[0]
punto[0] = 0
td[0] = 0
for i in range(1,k):
    punto[i] = i
    td[i] = deltaT*i
    xdigital[i] = xdigital[i-1]+yentrada[i]*deltaY

El resultado puede ser observado en los vectores, en una gráfica de los vectores y en un archivoaudio.wav.

# SALIDA
print('entrada: ')
print(yentrada)
print('salida: ')
print(xdigital)
entrada: 
[ 0  1 -1  1 -1  1  1 -1  1 -1  1  1 -1  1 -1  1 -1  1  1 -1  1 -1  1 -1  1 -1  1 -1  1 -1  1 -1  1 -1 -1  1 -1  1 -1  1 -1 -1  1 -1  1 -1 -1  1 -1  1 -1 -1  1 -1  1 -1 -1  1 -1  1 -1 -1  1 -1  1 -1  1 -1 -1  1 -1  1 -1  1 -1  1 -1  1 -1  1 -1  1 -1  1  1 -1  1 -1  1]
salida: 
[     0   6801      0   6801      0   6801  13602   6801  13602   6801   13602  20403  13602  20403  13602  20403  13602  20403  27204  20403 27204  20403  27204  20403  27204  20403  27204  20403  27204  20403  27204  20403  27204  20403  13602  20403  13602  20403  13602  20403  13602   6801  13602   6801  13602   6801      0   6801      0   6801  0  -6801   0  -6801      0  -6801 -13602  -6801 -13602  -6801  -13602 -20403 -13602 -20403 -13602 -20403 -13602 -20403 -27204 -20403  -27204 -20403 -27204 -20403 -27204 -20403 -27204 -20403 -27204 -20403  -27204 -20403 -27204 -20403 -13602 -20403 -13602 -20403 -13602]
# Graficar
plt.figure(1)       # define la grafica
plt.suptitle('Decodifica Delta-Sigma')

plt.subplot(211)    # grafica de 2x1 y subgrafica 1
plt.ylabel('entrada: y[n]')
plt.axis((0,k-1,-1.1,1.1)) # Ajusta ejes
plt.step(punto,yentrada, where='post')

plt.subplot(212)    # grafica de 2x1 y subgrafica 2
plt.ylabel('salida: x[t]')
xmax = np.max(xdigital)+0.1*np.max(xdigital)# rango en el eje y
xmin = np.min(xdigital)-0.1*np.max(xdigital)
plt.axis((0,td[k-1],xmin,xmax)) # Ajusta ejes
plt.step(td,xdigital, where='post', color='m')

plt.show()

Archivo de Audio del decodificador

Para crear el archivo de audio que permita escuchar el resultado del decodificador, se utiliza una instrucción de scipy que require:

  • el nombre del archivoaudio.wav: ‘sigmadeltaaudioruido.wav’
  • la frecuencia de muestreo del sonido: muestreo
  • el arreglo con la señal digital reconstruida: xdigital

El archivo.wav creado puede ejecutarse con windows media player:
sigmadeltaaudio440Hz_1s.wav

# Salida
# Archivo de audio.wav
muestreo = int(1/deltaT)
waves.write(archivoaudio, muestreo, xdigital)
print(' ... ' + archivoaudio + ' ...')
... sigmadeltaaudio.wav ...

windowsmediaplayer

Sigma-Delta Codificador con Python

Referencia: Leon-Couch, 3–8 Modulación Delta, p.192; Delta-sigma_modulation, Wikipedia ; Sigma-Delta – Modulación

La señal de entrada es el archivo de sonido:

440Hz_44100Hz_16bit_05sec.wav

que es tipo monofónica. La modulación delta-sigma se aplica usando los parámetros:

  • ΔY = deltaY = 0.3 * max(sonido)
  • rango de observación [0, 0.002] segundos para la gráfica
  • muestras en el rango de observación: k

Nota: En el caso de sonido estéreo, con dimensión kx2, se usa solo un canal (sonido[:,0])

Hay que tomar en cuenta que para graficar, la misma se satura con un número grande de muestras . Para escuchar el sonido el tiempo de ‘termina’ será de al menos 1 segundo.
codifica Delta Sigma

Algoritmo en Python

A continuación se describe el algoritmo para codificar en Python, que el el bloque de ingreso especifica el rango de observación en segundos.

% matplotlib inline
# Modulacion Delta Sigma - Codificador
# entrada x(t), salida: y[n]
# propuesta:edelros@espol.edu.ec
import numpy as np
import matplotlib.pyplot as plt
import scipy.io.wavfile as waves

# INGRESO 
# archivo = input('archivo de sonido:' )
archivo = '440Hz_44100Hz_16bit_05sec.wav'
muestreo, sonido = waves.read(archivo)

# rango de observación en segundos
inicia = 0
termina = 0.002
canal = 0   #Usar un canal en caso de estereo
c:\python34\lib\site-packages\scipy\io\wavfile.py:179: WavFileWarning: Chunk (non-data) not understood, skipping it. WavFileWarning)

Si el archivo.wav contiene etiquetas, el proceso de lectura emite una advertencia (WavFileWarning) que no afecta al proceso de codificación.

Se extrae solo una porción del sonido que contiene las k muestras en el intervalo de observación y se envía graficar.

La señal xdigital y ysalida se calcula con las diferencias entre el valor de cada muestra analógica y el valor xdigital anterior; se toma en cuenta solo el signo para acumular y generar el vector ysalida.

# PROCEDIMIENTO
# Codificar Sigma-Delta
deltaY = 0.1*np.max(sonido) 
deltaT = 1/muestreo 

# Usar un canal en caso de estereo
canales=sonido.shape
cuantos=len(canales)
canal = 0   
if (cuantos==1): # Monofónico
    uncanal=sonido[:]  
if (cuantos==2): # Estéreo
    uncanal=sonido[:,canal] 

# Extrae solo una porcion del sonido
donde = int(inicia/deltaT)
# tiempo muestreo de la señal analógica
t = np.arange(inicia,termina,deltaT) 
k = len(t) 

muestra = np.copy(uncanal[donde:donde+k])

# Señal Digital
xdigital = np.zeros(k, dtype=float) 
ysalida = np.zeros(k, dtype=int) 

for i in range(1,k):
    diferencia = muestra[i]-xdigital[i-1]
    if (diferencia>0):
        signo = 1
    else:
        signo = -1
    xdigital[i] = xdigital[i-1]+signo*deltaY
    ysalida[i] = signo
    
parametros=np.array([deltaT,deltaY,k])

Los resultados de pueden observar de dos maneras, en los archivos.txt de parámetros y datos, o en una gráfica.

# SALIDA
print('parámetros:[deltaT, deltaY, k]')
print(parametros)
print('datos:')
print(ysalida)
np.savetxt('deltasigma_parametros.txt',parametros)
np.savetxt('deltasigma_datos.txt',ysalida,fmt='%i')
print('... archivos.txt guardados ...')
parámetros:[deltaT, deltaY, k]
[  2.26757370e-05   6.80100000e+03   8.90000000e+01]
datos:
[ 0  1 -1  1 -1  1  1 -1  1 -1  1  1 -1  1 -1  1 -1  1  1 -1  1 -1  1 -1  1 -1  1 -1  1 -1  1 -1  1 -1 -1  1 -1  1 -1  1 -1 -1  1 -1  1 -1 -1  1 -1  1 -1 -1  1 -1  1 -1 -1  1 -1  1 -1 -1  1 -1  1 -1  1 -1 -1  1 -1  1 -1  1 -1  1 -1  1 -1  1 -1  1 -1  1  1 -1  1 -1  1]
... archivos.txt guardados ...

Gráficas

para graficar y no saturar el gráfico, se recomienda observar solo una porción de los datos determinada por ‘verdesde’ y ‘verhasta’ muestras desde [0,k-1]

verdesde=0
verhasta=90

# Graficar
plt.figure(1)       # define la grafica
plt.suptitle('Codificador Delta-Sigma')

plt.subplot(211)    # grafica de 2x1 y subgrafica 1
plt.ylabel('x(t), x[n]')
plt.plot(t[verdesde:verhasta],muestra[verdesde:verhasta], 'g')
plt.step(t[verdesde:verhasta],xdigital[verdesde:verhasta], where='post',color='m') # Puntos x[n]

plt.subplot(212)    # grafica de 2x1 y subgrafica 2
plt.ylabel('y[n]')
#plt.plot(ysalida, 'bo')     # Puntos y[n]
plt.axis((verdesde,verhasta,-1.1,1.1))
puntos=np.arange(verdesde,verhasta,1)     #posicion eje x para escalon
plt.step(puntos[verdesde:verhasta],ysalida[verdesde:verhasta], where='post')

plt.show()

codifica Delta Sigma
Archivos Resultantes: deltasigma_datos.txt, deltasigma_parametros.txt

4. Morse Decodificador de sonido

Morse – Decodificador de un mensaje Morse desde el sonido en archivo.wav

Para analizar el sonido de un mensaje Morse, se procede de forma semejante al decodificador de un tono morse, revisando la duración de cada símbolo: punto o espacio tienen duración de una marca, la raya tiene una duración de 3 marcas, los simbolos se separan por una pausa con un intervalo de una marca.

telegrafomarca

Se divide el sonido del archivo.wav en partes o ventanas, usando la funcion separaventanas().

ventanas de sonido

Cada ventana se analiza con la transformada de Fourier para encontrar la frecuencia del sonido mas fuerte o marcas de puntos '.', '-', ó ' '. El resultado de análisis de las ventanas es un vector con la marca de frecuencia más fuerte para cada una de las ventanas. (función marcasdeventanas())

La frecuencia o ‘tono frecuente’ se usa como una referencia para diferenciar si existe una señal para un punto o raya Morse. El tono frecuente se obtiene mediante la «moda» en el vector marcas usando scipy.stats.itemfreq() y usando la fila donde donde está la mayor cuenta usando np.argmax().

Para facilitar el análisis de un tono en una marca, se realiza una conversión a una secuencia binaria, usando '1' ó ,'0' si se encuentra el ‘tono frecuente’ en cada ventana.

Para el ejemplo, las ventanas se convierten al vector:
[1 1 1 1 1 1 0 0]

Al contar la duracion de cada intevalo, se determina si el sonido corresponde a un punto, una raya o una espacio:

Ejemplo:
. ... .--. --- .-..

Algoritmo en python

El sonido se obtiene de un archivo.wav del generador de tonos para un punto '.' o raya '-', del cual se leen los datos de sonido, y frecuencia de muestreo.

Para el ejemplo, el sonido del mensaje morse: morsetonoESPOL.wav

Las funciones para cada uno de los procesos descritos se encuentran a continuación

# Código Morse -  Determina un mensaje desde sonido.wav
# propuesta: edelros@espol.edu.ec
import numpy as np
import matplotlib.pyplot as plt
import scipy.io.wavfile as waves
import scipy.fftpack as fourier
import scipy.stats as stats

def separaventanas(sonido, muestreo, duracion = 0.04, partes = 8 ):
    # Divide el sonido en partes o ventanas para estudio
    # duracion = 0.04   # segundos de un punto
    # partes = 8        # divisiones de un punto
    # Extraer ventanas para espectro por partes
    tventana = duracion/partes
    mventana = int(muestreo * tventana) # muestras de una ventana

    # Ajuste de muestras de ventanas para matriz
    anchosonido  = int(len(sonido)/mventana)*mventana 
    sonidoajuste = np.resize(sonido,anchosonido)
    ventanas = np.reshape(sonidoajuste,(-1,mventana))  
    # -1 indica que calcule las filas
    return(ventanas)

def marcasdeventanas(ventanas, muestreo):
    # Analiza con FFT todas las ventanas del sonido
    filas, columnas = np.shape(ventanas)
    marcas = np.zeros(filas,dtype=int)

    # Espectro de Fourier de cada ventana
    # frecuencias para eje frq = fourier.fftfreq(mventana, 1/muestreo) 
    frq = fourier.fftfreq(columnas, 1/muestreo)  
    for f in range(0,filas,1):
        xf = fourier.fft(ventanas[f])
        xf = np.abs(xf)        # magnitud de xf
        tono = np.argmax(xf) # tono, frecuencia mayor 
        marcas[f] = frq[tono]
    return(marcas)

def secuenciabinaria(marcas): 
    # Busca el tono frecuente mayor que cero
    tonosventana = stats.itemfreq(marcas)
    donde = np.argmax(tonosventana[1:,1])
    tonopunto = tonosventana[donde+1,0]

    # Convierte los tonos a secuencia binaria
    secuencia = ''
    for valor in marcas:
        if (valor==tonopunto):
            secuencia = secuencia + '1'
        else:
            secuencia = secuencia + '0'
    return(secuencia)

def duracionmarcas(secuencia):
    # Duración de cada marca
    conteo = []
    caracter = '1'
    if (len(secuencia)>0):
        caracter = secuencia[0]
    i = 0
    k = 0
    while (i<len(secuencia)):
        if (secuencia[i]==caracter):
            k = k + 1
        else:
            conteo.append([int(caracter),k])
            if (caracter=='1'):
                caracter = '0'
            else:
                caracter = '1'
            k = 1
        i = i+1
    conteo = np.array(conteo)
    return(conteo)

def marcasmorse(conteo):
    # Determina la base de un tono
    veces = stats.itemfreq(conteo[:,1])
    donde = np.argmax(veces[:,1])
    base  = veces[donde,0]

    # Genera el código Morse
    tolera = 0.2 # Tolerancia en relacion
    bajo   = 1-tolera
    alto   = 1+tolera
    morse  = ''
    for j in range(0,len(conteo)):
        relacion = conteo[j,1]/base
        simbolo  = conteo[j,0]
        if (simbolo==1):
            if (relacion>(1*bajo) and relacion<(1*alto)):
                morse = morse + '.'
            if  (relacion>(3*bajo) and relacion<(3*alto)):
                morse = morse + '-'
        if  simbolo==0 :
            if (relacion>(3*bajo) and relacion<(3*alto)):
                morse = morse + ' '
            if (relacion>(7*bajo) and relacion<(7*alto)):
                morse=morse+'   '
    return(morse)

# INGRESO 
# archivo = input('nombre del archivo: ')
archivo = 'morsetonoESPOL.wav'
muestreo, sonido = waves.read(archivo)

# PROCEDIMIENTO

ventanas  = separaventanas(sonido, muestreo)
marcas    = marcasdeventanas(ventanas, muestreo)
secuencia = secuenciabinaria(marcas)
conteo    = duracionmarcas(secuencia)
morse     = marcasmorse(conteo)

# SALIDA
print('codigo en morse: ')
print(morse)
codigo en morse: 
. ... .--. --- .-..   .. -- .--. ..- .-.. ... .- -. -.. ---   .-.. .-   ... --- -.-. .. . -.. .- -..   -.. . .-..   -.-. --- -. --- -.-. .. -- .. . -. - ---

Revisión de valores en el procedimiento

Se presentan algunos de los valores intermedios para observar el proceso de detección de símbolos Morse.

# Salida /Observación intermedia
print('ventanas analizadas: ',len(marcas))
print('marcas: ', marcas)
ventanas analizadas:  3287
marcas:  [400 400 400 ...,   0   0   0]

Para facilitar el proceso de detección morse se convierten las marcas a una secuencia binaria, usando como referencia la frecuencia del tono de un punto determinada en la sección anterior.

# Salida /Observación intermedia
if len(secuencia)<100:
    print(secuencia)
else:
    print(secuencia[0:200] + ' ... ')
11111111000000000000000000000000111111110000000011111111000000001111111110000000000000000000000011111111000000001111111111111111111111111000000011111111111111111111111110000000111111111000000000000000 ... 

Cambiar la secuencia a simbolos morse, consiste en determinar la relación entre las veces que aparecen los ‘1’s y ‘0’s, el resultado se puede guardar en un arreglo.

La base de cuántas marcas corresponden a un punto de estima como la «moda» de las veces en que se repite un uno o un cero.

print('simbolo, cuenta:')
print(conteo)
simbolo, cuenta:
[[ 1  8]
 [ 0 24]
 [ 1  8]
 [ 0  8]
 [ 1  8]
 [ 0  8]
 [ 1  9]
 [ 0 23]
 [ 1  8]
 [ 0  8]
 [ 1 25]
 [ 0  7]
 [ 1 25]
...
# Salida
print('codigo en morse: ')
print(morse)
codigo en morse: 
. ... .--. --- .-..   .. -- .--. ..- .-.. ... .- -. -.. ---   .-.. .-   ... --- -.-. .. . -.. .- -..   -.. . .-..   -.-. --- -. --- -.-. .. -- .. . -. - ---

Con el resultado puede usar el decodificador morse de los temas anteriores.

Referencia: Código Morse Wikipedia, Recommendation ITU-R M.1677-1 (10/2009) International Morse code, Leon-Couch, 5–9 Señalización Pasabanda Modulada Binaria (OOK)

3. Morse Decodificador de untono

Referencia: Código Morse Wikipedia, Recommendation ITU-R M.1677-1 (10/2009) International Morse code, Leon-Couch, 5–9 Señalización Pasabanda Modulada Binaria (OOK)

Morse – Decodificador de sonido de Un simbolo

Analisis de archivo.wav en python

Para analizar el sonido de un símbolo Morse, se inicia por recordar la duración de cada símbolo:

Para la detección de un tono, se revisa la frecuencia activa en varios intervalos de tiempo o ventanas.

telegrafomarca

Al contar el número de intevalos (marcas) que aparece un tono vs. el número de intervalos de pausa, se podrá determinar si el sonido corresponde a un punto, una raya o un espacio. Es el proceso contrario a lo realizado para crear el tono, contanto las marcas, la relacion entre marcas y pausas indicará el símbolo recibido.

Algoritmo en python

El sonido para el análisis se obtiene de un archivo.wav, creado con el generador de tonos para un punto '.' o raya '-'. La lectura del archivo proporciona los datos de sonido, y frecuencia de muestreo.

Para el análisis, se supondrá que la ventana tiene la mitad de la duración de un punto, para el ejercicio es 0.04 segundos. Así, una ventana tiene la mitad de muestras que una marca '.' o pausa, obteniendo al menos dos ventanas por cada punto o pausa y tres ventanas por cada raya. Se cuentan las ventanas que tienen sonido o pausa y se comparan para determinar el símbolo.

Que existan al menos dos ventanas por cada punto tienen relación con el muestreo de Nyquist.

Lo expuesto se prueba primero analizando solo una ventana:

# Código Morse -  Determina un simbolo a partir de UN sonido
# propuesta: edelros@espol.edu.ec
import numpy as np
import matplotlib.pyplot as plt
import scipy.io.wavfile as waves
import scipy.fftpack as fourier

# INGRESO 
# archivo = input('nombre del archivo: ')
# archivo = 'morsetonopunto.wav'
archivo = 'morsetonoraya.wav'
muestreo, sonido = waves.read(archivo)

# PROCEDIMIENTO
duracion = 0.04   # segundos de un punto
partes   = 2        # divisiones de un punto

# Extraer ventanas para espectro por partes
tventana = duracion/partes
mventana = int(muestreo * tventana) # muestras de una ventana
# Observacion intermedia
unaventana = sonido[0:mventana]

# SALIDA  GRAFICA
print('muestreo:', muestreo)
print('muestras en sonido', len(sonido))
print('intervalo de ventana: ', tventana)
print('muestras por ventana: ', mventana)

plt.figure(1)
plt.subplot(211)
plt.plot(sonido, label='sonido')
plt.xlabel('muestras')
plt.ylabel('sonido y ventanas')
k = 1
for i in range(0,int(len(sonido)/mventana)):
    plt.axvline(x=k*mventana, color='r', linestyle='--')
    k = k + 1
plt.subplot(212)
plt.plot(unaventana)
plt.xlabel('muestras')
plt.ylabel('unaventana')
plt.margins(0)
plt.show()
muestreo: 11025
muestras en sonido 1764
intervalo de ventana: 0.02
muestras por ventana: 220

Análisis de una ventana

Para determinar si existe un tono en una ventana, se analiza la existencia de una señal observada en el dominio de la frecuencia.

En el dominio de la frecuencia se requiere usar la «Transformada de Fourier» de la señal en la ventana. La Transformada rápida de Fourirer fft() se encuentra disponible en las librerias «scipy» y el rango de frecuencias para la gráfica se obtiene con fftfreq().

En el resultado se busca la frecuencia (np.argmax()) con mayor magnitud (np.abs()) y se determina si corresponde al «tono» morse esperado.

# Espectro de Fourier de una ventana [0:nventana]
xf = fourier.fft(unaventana)
xf = np.abs(xf)  # magnitud de xf
tono = np.argmax(xf)  # tono en donde

# frecuencias para eje
frq = fourier.fftfreq(mventana, 1/muestreo ) 

# SALIDA  GRAFICA /Observacion intermedia
print('frecuencia tono (Hz): ',frq[tono])

plt.figure(2)
plt.xlabel('Frecuencia Hz')
plt.ylabel('Magnitud')
plt.stem(frq,xf)
plt.show()

frecuencia tono (Hz):  451.022727273

Analisis del archivo.wav de un símbolo Morse

Para analizar todo el archivo de sonido, los datos se segmentan por «ventanas» en una matriz, cada fila corresponde a las muestras de una ventana de tiempo.

La matriz debe ser rectangular, por lo que es necesario confirmar que el número de muestras por ventana sean iguales. En caso de no ser así, se realiza un ajuste del ancho del sonido para que el numero de muestras en el archivo sea múltiplo del número de muestras por ventana.

Observar el resultado de crear varias ventanas de tiempo, es semejante a tener un cuaderno en el que cada hoja tiene una imagen con los valores de la señal en cada ventana.

# Ajuste de muestras de ventanas para matriz
antes = len(sonido)
anchosonido = int(len(sonido)/mventana)*mventana 
sonido   = np.resize(sonido,anchosonido)
ventanas = np.reshape(sonido,(-1,mventana))  
# -1 indica que calcule las filas

# SALIDA
print('muestras en sonido: ', antes)
print('muestras recortadas a: ', anchosonido)
print('ventanas: ')
print(ventanas)

muestras en sonido:  1764
muestras recortadas a:  1760
ventanas: 
[[     0   6504  12602 ..., -22161 -24942 -26163]
 [-25748 -23722 -20212 ...,   9594   3241  -3315]
 [ -9663 -15408 -20188 ...,  25762  26158  24919]
 ..., 
 [-15377  -9629  -3278 ...,  -9733 -15468 -20236]
 [-23738 -25755 -26161 ...,      0      0      0]
 [     0      0      0 ...,      0      0      0]]

Para observar el resultado se usa una gráfica 3d de superficie.

# Observacion intermedia
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

nfilas,ncolumnas = np.shape(ventanas)
x = np.arange(0,ncolumnas) # valores para cada eje
y = np.arange(0,nfilas)
X, Y = np.meshgrid(x, y) # valores para punto.
Z = ventanas             # Matriz

# Salida, gráfico 3d
fig  = plt.figure(3,figsize=plt.figaspect(0.5))
ax   = fig.add_subplot(1, 2, 1, projection='3d')
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1,
                       cmap=cm.coolwarm, linewidth=0,
                       antialiased=False)
plt.xlabel('muestras')
plt.ylabel('ventana')
plt.show()

El proceso de analisis se repite para cada ventana (una fila) en la matriz, guardando cada resultado en un vector «marcas«.

# Analiza con FFT todas las ventanas del sonido
filas, columnas = np.shape(ventanas)
marcas = np.zeros(filas,dtype=float)

# Espectro de Fourier de cada ventana
# frecuencias para eje
frq = fourier.fftfreq(mventana, 1/muestreo)  
for fila in range(0,filas,1):
    xf = fourier.fft(ventanas[fila])
    xf = np.abs(xf)        # magnitud de xf
    tono = np.argmax(xf)   # tono, frecuencia mayor 
    marcas[fila] = frq[tono]

Observamos el resultado en una gráfica 2D.

# SALIDA  GRAFICA  /Observacion intermedia
print(marcas)
plt.figure(4)
plt.stem(marcas)
plt.xlabel('número de marca')
plt.ylabel('frecuencia del tono mas fuerte.')
plt.show()
[ 451.02272727  451.02272727  451.02272727  451.02272727  451.02272727  451.02272727    0.            0.        ]

Para determinar si es punto '.', raya'-' o espacio ' ', se cuentan los noceros y los ceros del arreglo de marcas, la relación entre ellos determina el símbolo que representa el sonido del archivo.wav.

# determina el simbolo
noceros = np.count_nonzero(marcas)
ceros  = len(marcas)-noceros
relacion = int(noceros/ceros)
if (relacion==3):
    simbolo = '-'
if (relacion==1):
    simbolo = '.'
if (relacion==0):
    simbolo = 'espacio'

# SALIDA
print('relacion noceros/ceros:', relacion)
print('simbolo morse detectado: ', simbolo)
relacion noceros/ceros: 3
simbolo morse detectado:  -

Nota: Observe que cualquier frecuencia en una ventana marca un punto '.'. El detalle es tratado en «Morse – Decodificador del sonido de un mensaje» en : # Busca frecuencia del tono punto.


Programa resumido

El resultado del programa resumido se muestra a continuación:

# Código Morse -  Determina un simbolo a partir de UN sonido
# propuesta: edelros@espol.edu.ec
import numpy as np
import matplotlib.pyplot as plt
import scipy.io.wavfile as waves
import scipy.fftpack as fourier

# INGRESO 
# archivo = input('nombre del archivo: ')
# archivo = 'morsetonopunto.wav'
archivo = 'morsetonoraya.wav'
muestreo, sonido = waves.read(archivo)

# PROCEDIMIENTO
duracion = 0.04   # segundos de un punto
partes   = 2      # divisiones de un punto

# Extraer ventanas para espectro por partes
tventana = duracion/partes
mventana = int(muestreo * tventana) # muestras de una ventana

# Ajuste de muestras de ventanas para matriz
antes = len(sonido)
anchosonido = int(len(sonido)/mventana)*mventana 
sonido   = np.resize(sonido,anchosonido)
ventanas = np.reshape(sonido,(-1,mventana))  
# -1 indica que calcule las filas

# Analiza con FFT todas las ventanas del sonido
filas, columnas = np.shape(ventanas)
marcas = np.zeros(filas,dtype=float)

# Espectro de Fourier de cada ventana
# frecuencias para eje
frq = fourier.fftfreq(mventana, 1/muestreo)  
for f in range(0,filas,1):
    xf = fourier.fft(ventanas[f])
    xf = np.abs(xf)        # magnitud de xf
    tono = np.argmax(xf) # tono, frecuencia mayor 
    marcas[f] = frq[tono] 

# determina el simbolo
noceros  = np.count_nonzero(marcas)
ceros    = len(marcas)-noceros
relacion = int(noceros/ceros)
if (relacion==3):
    simbolo = '-'
if (relacion==1):
    simbolo = '.'
if (relacion==0):
    simbolo = 'espacio'

# SALIDA
print('relacion de marcas:', relacion)
print('simbolo morse detectado: ', simbolo)
relacion de marcas: 3
simbolo morse detectado:  -