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()

 

3Eva_IIT2017_T4 Sx(f) para QAM

3ra Evaluación II Término 2017-2018.  Febrero 20, 2018

Tema 4. (25 puntos) La señal de entrada X(t) de un sistema tipo “QAM” son procesos aleatorios A(t) y B(t) independientes con densidades espectrales de potencia mostradas.

X(t) = A(t) cos(2πfct + θ) + B(t) sin(2πfct + θ)

a)     Encuentre la densidad espectral de potencia de la señal QAM, SX(f)

b)     Grafique su respuesta


 

Rúbrica: literal a (15 puntos), literal b (10 puntos)

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()

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))

Autocovarianza en AM

Referencia: León García Ejemplo 9.9 p495, Gubner Ej10.35 p496

Sea X(t) = A cos(2πt), donde A es una variable aleatoria, con un comportamiento semejante a la figura.

Encontrar el valor esperado , la autocorrelación y autocovarianza de de X(t).


El valor esperado se calcula a continuación, note que la media varia respecto a t y que el valor es cero para valores de t donde cos(2πt) =0.

E[X(t)] = E[A cos(2\pi t)] E[X(t)] = E[A] cos(2\pi t)

La autocorrelación es:

R_X(t_1,t_2) = E[A cos(2\pi t_1) A cos(2\pi t_2)] = E[A^{2} cos(2\pi t_1) cos(2\pi t_2)] = E[A^{2}] cos(2\pi t_1) cos(2\pi t_2)

usando:

2 cos(x)cos(y) = cos(x-y) + cos(x+y) cos(x)cos(y) = \frac{ cos(x-y) + cos(x + y)}{2}

se reemplaza:

= E[A^{2}] \frac{1}{2}[cos(2\pi t_1 - 2\pi t_2) + cos(2\pi t_1 + 2\pi t_2)] R_X(t_1,t_2) = E[A^{2}] \frac{[cos(2\pi (t_1 - t_2)) + cos(2\pi (t_1 + t_2))]}{2}

se observa que el valor de autocorrelación depende de las diferencias de tiempo t1 y t2.


La autocovarianza es:

Cov_X(t_1,t_2) = R_X(t_1,t_2) - E[X(t_1)]E[X(t_2)] = E[A^{2}] cos(2\pi t_1) cos(2\pi t_2) - E[A] cos(2\pi t_1)E[A] cos(2\pi t_2) = E[A^{2}] cos(2\pi t_1) cos(2\pi t_2) - E[A]^2 cos(2\pi t_1)cos(2\pi t_2) = (E[A^{2}] - E[A]^2) cos(2\pi t_1)cos(2\pi t_2) = Var[A] cos(2\pi t_1)cos(2\pi t_2)

con el mismo procedimiento de cos(x)cos(y):

Cov_X(t_1,t_2) = Var[A] \frac{[cos(2\pi (t_1 - t_2)) + cos(2\pi (t_1 + t_2))]}{2}

Autocorrelación, Autocovarianza con variable tiempo

Referencia: León-García 9.2.2 p493

Se puede usar los momentos de muestras en el tiempo para parcialmente especificar un proceso aleatorio al resumir la información contenida en las cdf conjuntas.


Procesos aleatorios Contínuos en tiempo

Media:

m_X(t) = E[X(t)] = \int_{-\infty}^{\infty} x f_{X(t)}(x) dx

Varianza:

VAR[X(t)] = \int_{-\infty}^{\infty} ( x - m_X(t))^2 f_{X(t)}(x) dx

donde f_{X(t)}(x) es la pdf de X(t). Note que ambas son funciones determinísticas de tiempo.

Autocorrelación

R_X(t_1,t_2) = E[X(t_1,t_2)] = = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} xy f_{X(t_1),X(t_2)}(x,y) dx dy

Autocovarianza:

C_X(t_1,t_2) = E[\{X(t_1) - m_X(t_1) \} \{ X(t_2) - m_X(t_2)\} C_X(t_1,t_2) = R_X(t_1,t_2) - m_X(t_1) m_X(t_2)

coeficiente de correlación

\rho (t_1, t_2) = \frac{C_x(t_1,t_2)}{\sqrt{C_x(t_1,t_1)}\sqrt{C_x(t_2,t_2)}}

El coeficiente de correlación es la medida en la cual una variable aleatoria puede predecirse como una función lineal de otra.


Procesos aleatorios Discretos en tiempo

Media:

m_X(n) = E[X(n)]

Varianza:

VAR[X(n)] = E[(X(n) - m_X(n))^2]

Autocorrelación

R_X(n_1,n_2) = E[X(n_1,n_2)]

Autocovarianza:

C_X(n_1,n_2) = E[\{X(n_1) - m_X(n_1) \} \{ X(n_2) - m_X(n_2)\} C_X(n_1,n_2) = R_X(n_1,n_2) - m_X(n_1) m_X(n_2)