Correlación con Python – Ejercicio Estación Meteorológica

Referencia: Ross 8.2 p334

De una estación meteorológica se obtiene un archivo.csv con los datos de los sensores disponibles.

2021_10_EstMeteorologica.csv

 

Nota: las instrucciones para la lectura del archivo se describen en:

Archivos.csv con Python – Ejercicio con gráfica de temperatura y Humedad

Para observar la correlación entre las medidas obtenidas de varios sensores, se considera usar las variables:

sensor_EM = ["TEMP","Humidity","Solar_rad","Bar_press."]

La correlación entre las variables se realiza seleccionando las columnas desde una lista, para luego realizar la operación desde una de las funciones de Pandas-Python

# Matriz de Correlación
mediciones = tablaEM[sensor_EM].copy()
correlacion_matriz = mediciones.corr()

obteniendo como resultados:

Matriz de correlación:
                TEMP  Humidity  Solar_rad  Bar_press.
TEMP        1.000000 -0.957730   0.743379   -0.387187
Humidity   -0.957730  1.000000  -0.765766    0.447845
Solar_rad   0.743379 -0.765766   1.000000   -0.093236
Bar_press. -0.387187  0.447845  -0.093236    1.000000

Del resultado se observa que el coeficiente de correlación entre las mediciones de «TEMP» y «Humidity» es alta. Para observar el caso se usa la gráfica entre las dos variables comparando los datos en el tiempo por fechas.
Cuando la temperatura aumenta la Humedad disminuye, lo opuesto tambien es observable.

Comparando solo entre las dos variables:

Para el caso de radiación solar y presión baromérica, el coeficiente de correlación es «bajo».

comparando entre las dos variables solamente:

Instrucciones en Python

# Archivos de estación meteorologica
# correlacion entre sensor_EM
import numpy as np
import pandas as pd
import datetime as dt
import os
import matplotlib.pyplot as plt
from matplotlib.dates import DateFormatter, DayLocator

# INGRESO
archivoEM = "2021_10_EstMeteorologica.csv"
sensor_EM = ["TEMP","Humidity","Solar_rad","Bar_press."]
sensor_EM_u = ["(°C)","(%)","(kW)","(hPa)"]

# PROCEDIMIENTO
# tabla estacion meteorologica
tablaEM = pd.read_csv(archivoEM, sep=';',decimal=',')
tablaEM = pd.DataFrame(tablaEM)
tablaEM_n = len(tablaEM)

# Matriz de Correlación
mediciones = tablaEM[sensor_EM].copy()
correlacion_matriz = mediciones.corr()

# SALIDA
print('Matriz de correlación:')
print(correlacion_matriz)

Para añadir las gráficas, se añaden instrucciones para seleccionar las variables en intervalos de fechas.

# tabla estacion meteorologica
fechaformatoEM = "%d/%m/%Y %H:%M:%S"
# fechas concatenando columnas de texto
tablaEM['fecha'] = tablaEM['Date']+' '+tablaEM['Time']
tablaEM['fecha'] = pd.to_datetime(tablaEM['fecha'],
                                  format=fechaformatoEM)
# grafica
sensorA = "TEMP"
sensorB = "Humidity"
fechaA = "2021/10/01"
fechaB = "2021/10/07"

ti = tablaEM["fecha"].loc[tablaEM["fecha"].between(fechaA,fechaB)]
xi = tablaEM[sensorA].loc[tablaEM["fecha"].between(fechaA,fechaB)]
yi = tablaEM[sensorB].loc[tablaEM["fecha"].between(fechaA,fechaB)]

# grafica compara sensores en tiempo
fig_Pt, graf1 = plt.subplots()
graf1.scatter(ti,xi, marker = '.',
              label=sensorA)

eje12 = graf1.twinx()  # segundo eje del grafico
eje12.scatter(ti,yi, marker = '.',
              color='orange',label=sensorB)
graf1.set_ylabel(sensorA, color='blue')
eje12.set_ylabel(sensorB, color='orange')
plt.show()

# grafica compara sensores
plt.scatter(xi,yi)
plt.xlabel(sensorA, color='blue')
plt.ylabel(sensorB, color='orange')
plt.show()

 

Correlación bits con ruido

Dada una secuencia de bits, al transmitirlos se distorsionó con ruido aditivo. En el receptor será dificil discriminar los valores de los bits tomados de la señal recibida.

La correlación permite «limpiar» un poco la señal en el receptor y una mejor estimación de los bits recibidos.

Para el ejemplo se usa una secuencia de bits 01101001 como señal a transmitir, el ruido aditivo es de tipo normal con media m=0 y varianza σ2=1

Para referencia del proceso, se usa un punto rojo en el centro de la ventana de tiempo de cada bit.

La correlación se realiza con un bit y la señal con ruido.

Observe que luego de la correlación es más sencillo discriminar si lo recibido fué un bit 0 o un bit 1.

señal de bits sin ruido
[0 1 1 0 1 0 0 1]
estimado de bits en el receptor
[ 0.002  0.985  1.07   0.054  0.947 -0.079 -0.006  1.149]

Instrucciones en Python

# secuencia de bits con Ruido
# usar correlación para mejorar la senal
# estimar el bit recibido
# Tarea: convertir estimado en bits

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as st

# INGRESO

senalbit = np.array([0, 1, 1, 0, 1, 0, 0, 1])

# PROCEDIMIENTO

# señal moduladora son escalones por cada bit
# o forma rectangular de duración anchobit
anchobit = 128
mitadbit = anchobit//2
senal = np.repeat(senalbit, anchobit)
muestras=len(senal)

# reloj para observar la mitad del bit
reloj = np.arange(mitadbit, muestras, anchobit)

# Ruido normal o Gausiano
media = 0
varianza = 1
ruido = st.norm.rvs(media,varianza,muestras)

# Añade el ruido a la señal
senalruido = senal + ruido

# Referencia de un bit para correlación
modelo = np.ones(anchobit)

# correlación modelo de bits
correlacion = np.correlate(senalruido, modelo, mode='same')
# normaliza el resultado [0,1]
correlacion = correlacion / anchobit

# estimado de bits en el receptor
estimado=correlacion[reloj]

# SALIDA
print('señal de bits sin ruido')
print(senalbit)
print('estimado de bits en el receptor')
np.set_printoptions(precision=3)
print(estimado)

#GRAFICA
plt.figure(1)       # define la grafica
plt.suptitle('señal binaria con ruido gaussiano')

plt.subplot(311)
plt.plot(senal)
plt.plot(reloj, senal[reloj], 'ro')
plt.ylabel('señal original')
plt.margins(0,0.05)

plt.subplot(312)
plt.plot(senalruido)
plt.plot(reloj, senalruido[reloj], 'ro')
plt.ylabel('senal con ruido')
plt.margins(0,0.05)

plt.subplot(313)
plt.plot(correlacion)
plt.ylabel('correlación')
plt.plot(reloj, correlacion[reloj], 'ro')
plt.margins(0,0.05)

plt.show()

pmf – Modulación QPSK

QPSK (Quadrature Phase-Shift Keying)

Este esquema de modulación es conocido también como Quaternary PSK (PSK Cuaternaria), Quadriphase PSK (PSK Cuadrafásica).

Esta modulación digital es representada en el diagrama de constelación por cuatro puntos equidistantes del origen de coordenadas.

Con cuatro fases, QPSK puede codificar dos bits por cada símbolo.

Respecto a un ancho de banda predeterminado, la ventaja de QPSK sobre BPSK está que con el primero se transmite el doble de la velocidad de datos en un ancho de banda determinado en comparación con BPSK, usando la misma tasa de error.

En el caso de la canción procesada en BPSK, se cargan una cantidad de datos que se pueden procesar con la mitad de símbolos, enviando la información por pares.

La pmf del proceso QPSK será:


Resultados del algoritmo:

datos cargados:  8595119
símbolos procesados:  4297559
[[ 539083       0 1580668]
 [      1       0       0]
 [1638724       0  539083]]
pmf[x,y]
[[  1.25439348e-01   0.00000000e+00   3.67806003e-01]
 [  2.32690232e-07   0.00000000e+00   0.00000000e+00]
 [  3.81315067e-01   0.00000000e+00   1.25439348e-01]]
>>> 

Instrucciones en Python

# Modulacion digital QPSK - pmf
# propuesta:edelros@espol.edu.ec
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
# archivo = input('archivo de delta-sigma:' )
narchivo = 'elaguacate_deltasigma_datos.txt'
senal = np.loadtxt(narchivo,dtype=int)

# PROCEDIMIENTO
n = len(senal)

simbolos = [-1,0,1]
m = len(simbolos)

# Codificar de 2 en dos
agrupar = 2
cuenta  = np.zeros(shape=(m,m), dtype=int)
nmax = (n//agrupar)*agrupar
for i in range(0,nmax,agrupar):
    a = senal[i]
    b = senal[i+1]
    f = simbolos.index(a)
    c = simbolos.index(b)
    cuenta[f,c] = cuenta[f,c]+1

k = np.sum(cuenta)
pxy = cuenta/k

# SALIDA
print('datos cargados: ', n)
print('símbolos procesados: ', k)
print(cuenta)
print('pmf[x,y]')
print(pxy)


# Gráfica:
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

xpos, ypos = np.meshgrid(simbolos,simbolos)
xpos = xpos.flatten('F')

ypos = ypos.flatten('F')
zpos = np.zeros_like(xpos)
dx = 0.8 * np.ones_like(zpos)
dy = dx.copy()
dz = pxy.flatten()
ax.bar3d(xpos, ypos, zpos, dx, dy, dz, color='b', zsort='average')
plt.show()

Tarea: Obtener las pmf marginales del ejercicio.

pmf – Modulación BPSK

La modulación por desplazamiento de fase o PSK (Phase Shift Keying) es una forma de modulación angular en que se modifica la fase de la portadora acorde a valores discretos.

BPSK (PSK Binario)

La modulación consiste en el desplazamiento de fase para 2 símbolos.
También conocida como 2-PSK o PRK (Phase Reversal Keying).

Es la modulación más sencilla por  emplear solo 2 símbolos, con 1 bit de información cada uno.

Los símbolos suelen tener un valor de salto de fase de 0º para el 1 y 180º para el 0 (-1), como se muestra en un diagrama de constelación.

En cambio, su velocidad de transmisión es la más baja de las modulaciones de fase.

BPSK – pmf

La pmf de BPSK muestra el uso de cada símbolo durante una transmisión. Por ejemplo: de un ejercicio previo se codificó a Sigma-Delta una canción teniendo como resultado:

elaguacate_deltasigma_datos.txt

elaguacate_deltasigma_parametros.txt


resultado del algoritmo:

cantidad de símbolos:  8595119
cuenta de símbolos: [4297559       1 4297559]
pmf de símbolos:  [  4.99999942e-01   1.16345102e-07   4.99999942e-01]
>>> 

Instrucciones en Python

# PMF de una señal Sigma-Delta
# propuesta:edelros@espol.edu.ec
import numpy as np
import matplotlib.pyplot as plt

# INGRESO 
# archivo = input('archivo delta-sigma:' )
archivo = 'elaguacate_deltasigma_datos.txt'
senal = np.loadtxt(archivo, dtype=int)

# PROCEDIMIENTO
n = len(senal)
simbolos = [-1,0,1]
m = len(simbolos)
cuenta = np.zeros(m, dtype=int)
for i in range(0,n,1):
    bit = senal[i]
    cual = simbolos.index(bit)
    cuenta[cual] = cuenta[cual]+1
pmf = cuenta/n

# SALIDA
print('cantidad de símbolos: ', n)
print('cuenta de símbolos:', cuenta)
print('pmf de símbolos: ', pmf)

# Gráfica
plt.stem(simbolos,pmf)
plt.title('pmf sigma-delta')
plt.xlabel('símbolos')
plt.ylabel('frecuencia relativa')
plt.show()

Suma y Retraso de procesos estocásticos

Referencia: León García ejemplos 10.4 y 10.5 p58, Gubner Ejemplo 10.18 p396

Suma de dos procesos

Encuentre la densidad espectral de potencia de Z(t) = X(t) + Y(t), donde X(t) y Y(t) son procesos conjuntamente estacionarios en el sentido amplio WSS.

Solución

Autocorelación de Z(t) es

R_Z(\tau) = E[Z(t + \tau )Z(t)] = = E[(X(t+ \tau)+ Y(t + \tau ))(X(t) + Y(t))] R_Z(\tau) = R_X(\tau) + R_{XY}(\tau) +R_{YX}(\tau) + R_Y(\tau)

La densidad espectral de potencia se calcula como:

S_Z(f) = Fourier\{ R_X(\tau) + R_{XY}(\tau) +R_{YX}(\tau) + R_Y(\tau) \} S_Z(f) = S_X(f) + S_{XY}(f) +S_{YX}(f) + S_Y(f)

Referencia: León García Ejemplo 10.5 p.582, Gubner Ejemplo 10.18 p396

Retraso

Sea Y(t) = X(t- d), donde d es una constante de retraso y donde X(t) es estacionario en el sentido amplio WSS.
Encuentre RYX(τ), SYX(f), RY(τ) y SY(f)

Solución

Usando las definiciones se tiene que:

R_{YX}(\tau) = E[Y(t + \tau )X(t)] = = E[X(t + \tau - d)X(t)] = R_{YX}(\tau) = R_X (\tau - d)

que usando la propiedad de desplazamiento en tiempo de la transformada de Fourier:

S_{YX}(f) = Fourier\{ R_X(\tau - d) \} = = S_X(f) e^{-j2 \pi fd} = = S_X(f) \cos (2\pi fd) - jS_x(f) sin (2 \pi fd)

Finalmente.

R_Y(\tau) = E[Y(t + \tau)Y(t)] = = E[X(t + \tau - d)X(t - d)] = R_X(\tau) S_Y(f) = Fourier\{ R_Y(\tau) \} = Fourier\{ R_X(\tau) \} = S_Y(f) = S_X(f)

Note que la densidad espectral de potencia crusada es compleka, y que SX(f) = SY (f) sin importar el hecho que X(t) ≠ Y(t). Entonces SX(f) = SY(f) lo que no implica que X(t) = Y(t)

Densidad espectral de Potencia de PM

Referencia: León-García Ejemplo 10.2 p581, Gubner Ejemplo 10.21 p399

Ejercicio

Sea X(t) = a cos(ω t + Θ), donde Θ es uniforme en el intervalo (0,2π) Encontrar la autocovarianza de X(t).
Encuentre la densidad espectral de potencia de SX(f):

Solución

S_X(f) = Fourier\{ R_X (\tau) \}

Autocorrelación

R_X (t_1,t_2) = E[(a \cos(\omega t_1 + \Theta)) (a \cos(\omega t_2 +\Theta))]

Recordando que:

E[g(x)] = \int_{-\infty}^{\infty} g(x) f(x) dx cos(x) cos(y) = \frac{cos(x-y) + cos(x+y) }{2}

se tiene que:

= \int_{-\pi}^{\pi} [a\cos(\omega t_1 + \Theta) a\ cos(\omega t_2 +\Theta)] \frac{1}{2\pi} d\Theta = \int_{-\pi}^{\pi} a^2 \frac{\cos(\omega (t_1 - t_2))+\cos(\omega (t_1 + t_2)+ 2\Theta)}{2} \frac{1}{2\pi} d\Theta = a^2 \int_{-\pi}^{\pi} \frac{cos(\omega (t_1 - t_2))}{2} \frac{1}{2\pi} d\Theta + a^2\int_{-\pi}^{\pi} \frac{cos(\omega (t_1 + t_2 )+ 2\Theta)}{2} \frac{1}{2\pi} d\Theta

El primer integral, el coseno no depende de Θ, mientras que el segundo integral es semejante al intergral de la media y cuyo resultado es cero.

= a^2 \left. \frac{cos(\omega (t_1 - t_2))}{2} \frac{\Theta}{2\pi} \right|_{-\pi}^{\pi} + 0 R_X (t_1,t_2) = \frac{a^2}{2} cos(\omega (t_1 - t_2)) R_X (t_1,t_2) = \frac{a^2}{2} cos(\omega \tau)

La densidad espectral de potencia entonces es:

S_X(f) = Fourier\{ R_X (\tau) \} = Fourier\{ \frac{a^2}{2} cos(\omega \tau) \} = \frac{a^2}{2} Fourier\{ cos(\omega \tau) \} = \frac{a^2}{2} Fourier\{ cos(2\pi f_0 \tau) \}

usando las tablas de transformadas de Fourier:

= \frac{a^2}{2} (\frac{1}{2}\delta (f-f_0) + \frac{1}{2} \delta (f+f_0)) = \frac{a^2}{4} \delta (f-f_0) + \frac{a^2}{4} \delta (f+f_0))

El promedio de potencia de la señal es RX(0) = a2/2.
De toda esta potencia, se concentra en las frecuencias en f0 positiva y negativa, por lo que la densidad espectral de potencia en esta frecuencias es infinita.

Densidad Espectral de Potencia Concepto

Referencia: León García 10.1.1 p578

Sea X(t) un proceso aleatorio contínuo en el tiempo estacionario en el sentido amplio WSS, con media mX y función de autocorrelación RX(t).

Si cambiamos el dominio de la función desde tiempo a frecuencia, lo anterior deberá también cambiarse de dominio. Si la función densidad de probabiliad se cambia de dominio, se conocería como periodograma estimador, y se llega a determinar la densidad espectral de potencia de X(t) definida como:

S_X (f) = lim_{T \rightarrow \infty} \frac{1}{T} E \left[ \left| \tilde{x} (f) \right| ^2\right]

La densidad espectral de potencia de X(t) está dada por la transformada de Fourier de la función de autocorrelación:

S_X (f) = Fourier\{R_x(\tau) \} = \int_{-\infty}^{\infty} R_x(\tau) e^{-j2\pi f \tau} d\tau

La potencia promedio de X(t) se expresa también como:

E[ X^2(t)] = R_X(0) = \int_{-\infty}^{\infty} S_X (f) df

La densiddad espectral de potencia también se relaciona con la autocorrelación y autocovarianza por medio de la transformada de Fourier:

S_X(f) = Fourier\{ C_x (\tau) + m^2_x\}

si consideramos que m_x es un componente constante o «DC»

S_X(f) = Fourier\{ C_x (\tau)\} + m^2_x \delta(f)

ampliando el concepto a densidad espectral de potencia cruzada SX,Y (f) se define como:

S_{X,Y}(f) = Fourier\{ R_{X,Y} (\tau) \}

Autocorrelación Propiedades

Referencia León-García 9.6 p522

Potencia Promedio

La función autocorrelación a τ=0 entrega la potencia promedio (segundo momento) del proceso:

R_{X}(0) = E[X(t)^2]

para todo valor de t

Función par respecto a τ

R_{X}(\tau) = E[X(t+\tau)X(t)] = E[X(t)X(t+\tau)] = R_{X}(-\tau)

Mide la tasa de cambio

La función de autocorrelación es una medida de la tasa de cambio de un proceso aleatorio

P[|X(t+\tau) - X(t)| > \epsilon] = P[(X(t+\tau) - X(t))^2 > \epsilon ^2] \leq \frac{E[(X(t+\tau) - X(t))^2]}{\epsilon ^2} =\frac{2\{R_X(0) - R_X(\tau) \}}{\epsilon ^2}

como resultado de usar la inequidad de Markov

tiene maximo en τ=0

si se usa la inequidad de Cauchy-Schwarz:

E[XY]^2 \leq E[X^2]E[Y^2] R_X(\tau )^2 = E[X(t+ \tau) X(t)]^2 \leq E[X^2(t+ \tau)] E[X^2(t)] = R_X(0) |R_X(\tau)| \leq R_X(0)

Periódica en media cuadrática

si R_X(0) = R_X(d) , entonces R_X(\tau) es periódica con periodo d y X(t) es «periódica en media cuadrática».

R_X(\tau + d)| = R_X(\tau)

se aproxima a el cuadrado de la media cuando τ tiende a infinito

Proceso aleatorio estacionario

Referencia: León-García 9.6, 9.6.1 p518,521; Gubner 10.3, 10.2 p395, p392; Ross 10.7 p654, p656

Estacionario

Un proceso aleatorio discreto o contínuo en el tiempo 
X(t) es estacionario si la distribución conjunta de 
cualquier grupo de muestras No depende de la ubicación 
del tiempo de origen.
F_{X(t_1),..., X(t_k)} (x_1, ... x_k) = F_{X(t_1+\tau),..., X(t_k)} (x_1, ... x_k+\tau)

para cualquier τ, k, o muestras t1, … , tk.

Estacionario con cdf de 1er Orden

Un proceso aleatorio estacionario con cdf de primer orden debe ser independiente del tiempo

F_{X(t)} (x) = F_{X(t+\tau)} (x) = F_X(x)

para todo t, τ, dicho de otra forma, sus resultados son constantes

m_{X}(t) = E[X(t)]=m VAR[X(t)] = E[(X(t)-m)^2] = \sigma ^2

Estacionario con cdf de 2do Orden

Un proceso aleatorio estacionario con cdf de segundo orden debe depender solo de la diferencia de tiempo entre muestras y NO en particular del tiempo de las muestras

F_{X(t_1),X(t_2)} (x_1,x_2) = F_{X(0),X(t_2-t_1)} (x_1,x_2)

para todo t1, t2

R_{X}(t_1,t_2) = R_{X} (t_2-t_1) C_{X}(t_1,t_2) = C_{X} (t_2-t_1)

Estacionario en el Sentido Amplio (WSS o débil)

 Un proceso aleatorio discreto o contínuo en el tiempo 
X(t) es estacionario en el sentido amplio (WSS) 
si satiface que: 
m_X(t) = m C_{X}(t_1,t_2) = C_{X} (t_2-t_1)

De forma semejante, los procesos X(t) y Y(t) son estacionarios conjuntos si ambos son estacionarios en el sentido amplio y si covarianza cruzada depende solo de t1-t2.

R_{X}(t_1,t_2) = R_{X} (\tau) C_{X}(t_1,t_2) = C_{X} (\tau)

Autocovarianza PM

Referencia: León García Ejemplo 9.10 p495, Gubner Ejemplo 10.8 p389, Gubner Ejemplo 10.17 p396

Sea X(t) = cos(ω t + Φ), donde Φ es uniforme en el intervalo (-π,π) Encontrar la autocovarianza de X(t).

la variable aleatoria Φ tiene distribución uniforme en el intervalo, por lo que la función fΦ(φ) es constante = 1/[π – (-π)] = 1/2π.

Recordamos que:

E[g(x)] = \int_{-\infty}^{\infty} g(x) f(x) dx

Media (León-García 4.15 p158):

m_X(t) = E[cos(\omega t + \Phi)] = = \int_{-\pi}^{\pi} cos(\omega t + \Phi) \frac{1}{2\pi} d\Phi = \left. \frac{-1}{2\pi} (sin (\omega t + \Phi)) \right|_{-\pi}^{\pi} = \frac{-1}{2\pi} [sin (\omega t + (-\pi)) - sin (\omega t + \pi)] = 0

Autocovarianza

dado que el valor esperado es cero, la autocovarianza es igual a la autocorrelación

C_{X} (t_1,t_2) = R_X (t_1,t_2) - E[X(t_1,t_2)] = R_X (t_1,t_2) = E[cos(\omega t_1 + \Phi) cos(\omega t_2 +\Phi)]

Recordando que:

E[g(x)] = \int_{-\infty}^{\infty} g(x) f(x) dx cos(x) cos(y) = \frac{cos(x-y) + cos(x+y) }{2}

se tiene que:

= \int_{-\pi}^{\pi} [cos(\omega t_1 + \Phi) cos(\omega t_2 +\Phi)] \frac{1}{2\pi} d\Phi = \int_{-\pi}^{\pi} \frac{cos(\omega (t_1 - t_2))+cos(\omega (t_1 + t_2)+ 2\Phi)}{2} \frac{1}{2\pi} d\Phi = \int_{-\pi}^{\pi} \frac{cos(\omega (t_1 - t_2))}{2} \frac{1}{2\pi} d\Phi + \int_{-\pi}^{\pi} \frac{cos(\omega (t_1 + t_2 )+ 2\Phi)}{2} \frac{1}{2\pi} d\Phi

El primer integral, el coseno no depende de Φ, mientras que el segundo integral es semejante al intergral de la media y cuyo resultado es cero.

= \left. \frac{cos(\omega (t_1 - t_2))}{2} \frac{\Phi}{2\pi} \right|_{-\pi}^{\pi} + 0
C_{X} (t_1,t_2) = R_X (t_1,t_2) =
= \frac{1}{2} cos(\omega (t_1 - t_2))