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

 

Referencia: Problema León García 10.10 p.636

Sean X(t) y Y(t) procesos independientes estacionarios en el sentido amplio.

Defina Z(t) = X(t) Y(t)

a) Muestre que Z(t) es estacionario en el sentido amplio (WSS).

b) Encuentre R_Z(τ) y S_Z(f)

Solución propuesta:

E[Z(t)] = E[X(t) Y(t)]

que por ser independientes,

= E[ X(t) ] E[ Y(t) ] = m_X m_Y

R_Z(τ) = E[ Z(t) Z(t+τ) ]
= E[ X(t)Y(t) X(t+τ)Y(t+τ) ]
= E[ X(t)X(t+τ) Y(t)Y(t+τ) ]
= E[ X(t)X(t+τ)] E[ Y(t)Y(t+τ) ]
= R_X(τ) R_Y(τ)

que tambien dependen solo de τ, por lo que Z(t) es WSS

S_Z(f) = F[R_Z(τ)]
= F[R_X(τ) R_Y(τ)]
= S_X(f) * S_Y(f) (convolución)

Referencia: Convolutions | Why X+Y in probability is a beautiful mess. 3Blue1Brown. 27-Junio-2023.

Referencia: Problema León García 10.6 p.636

Sea Z(t) = X(t) + Y(t)

¿Bajo qué condiciones S_Z(f) = S_X(f) + S_Y(f)?

Solución Propuesta:

E[ Z(t)Z(t+τ) ] = E[ [X(t) + Y(t)][X(t+τ) + Y(t+τ)] ]

= E[ X(t)X(t+τ) + Y(t) X(t+τ) + X(t)Y(t+τ) + Y(t)Y(t+τ) ]

= E[  X(t)X(t+τ) ] +E[ Y(t) X(t+τ) ] + E[ X(t)Y(t+τ) ] + E[ Y(t)Y(t+τ) ]

R_Z(τ) = R_X(τ) + R_YX(τ)  + R_XY(τ)  + R_Y(τ)

S_Z(f) = S_X(f) + S_YX(f)  + S_XY(f)  + S_Y(f)

Si X(t) y Y(t) son ortogonales, entonces:

R_XY(τ)  = R_YX(τ) = 0

S_Z(f) = S_X(f) + 0  + 0  + S_Y(f)

S_Z(f) = S_X(f) +  S_Y(f)

 

Referencia: Problema León García 10.4 p.635

a) Encuentre la función de autocorrelación correspondiente a la densidad espectral de potencia de la figura.

b) Encuentre el total de la potencia promedio

c) Grafique la potencia en el rango de |f| > f₀ como función de f₀>0.

Solución propuesta:

a)
S_X(f) = A \prod \Big(\frac{f}{2f_2} \Big) + (B-A) \bigwedge \Big( \frac{f}{f_1} \Big)

R_X(\tau) = F^{-1}\Big[ A \prod \Big(\frac{f}{2f_2} \Big) \Big] + F^{-1} \Big[ (B-A) \bigwedge \Big( \frac{f}{f_1} \Big) \Big] = 2A f_2 [Sa (2\pi f_2 \tau)] + (B-A) f_1 [Sa (\pi f_1 \tau) ]^2

b)
P = \int_{-\infty}^{\infty} S_X(f) \delta f

= A(2f_2) + (2f_1)\frac{(B-A)}{2} = 2A f_2 + (B-A)f_1

c)

la potencia en función de la frecuencia es par, por lo que se integra entre 0 y f₀ y se duplica para el rango entre [-f₀, f₀]
2\int_{0}^{f_0} S_X(f) \delta f =

primera sección: 0 < f₀ < f₁

2\int_{0}^{f_0} \Big[ \Big( -\frac{B-A}{f_1} \Big)f +B\Big] \delta f = 2 \Big[ \Big(-\frac{B-A}{f_1}\Big) \frac{ f^2}{2} +Bf \Big] \Big|_{0}^{f_0} = \Big(-\frac{B-A}{f_1}\Big) f_0^2 +2Bf_0 = 2Bf_0 -\frac{B-A}{f_1} f_0^2

segunda sección: f₁ < f₀ < f₂
= 2[ \frac{B-A}{2}f_1 + A(f_0 - f_1) \Big]

# leon- garcia 10.4
# literal c
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
A = 1
B = 2
f1 = 1
f2 = 2

n = 50
final = 4
# PROCEDIMIENTO
f = np.linspace(0,final,n)
P = np.zeros(n,dtype=float)
for i in range(0,n,1):
    if f[i]= f1 and f[i]f2:
        P[i] = 2*(((B+A)/2)*f1 + A*(f2-f1))

# SALIDA Grafica
plt.plot(f,P)
plt.vlines(f1,0,2.5*B, color='m', linestyle='dashed')
plt.vlines(f2,0,2.5*B, color='m', linestyle='dashed')
plt.show()

Referencia: Problema León García 10.3 p.635

a) Encuentre la densidad espectral de potencia S_Y(f) de un proceso aleatorio con función de autocorrelación R_X(τ) cos(2π f₀ τ), donde R_X(τ) es también una función de autocorrelación.

b) Grafique S_Y(f) si R_X(τ) como en el problema 10.1a.

Solución propuesta:

R_Y(\tau) = R_X(\tau) \cos(2\pi f_0 \tau) S_Y(f) = F\Big[ R_X(\tau) \cos(2\pi f_0 \tau) \Big] = F\Big[ R_X(\tau) \frac {e^{j2\pi f_0 \tau} + e^{-j2\pi f_0 \tau}}{2} \Big] = \frac{1}{2} F\Big[ R_X(\tau) e^{j2\pi f_0 \tau} \Big] + \frac{1}{2} F\Big[ R_X(\tau) e^{-j2\pi f_0 \tau} \Big] = \frac{1}{2} S_X(f-f_0) + \frac{1}{2} S_X(f+f_0)

donde S_X (f) = F[R_X(τ)]

1. PSD Rectangular

Referencia: Problema León García 10.2 p.635

Sea p(x) una función rectangular. ¿ R_X(τ) = p(τ/T) es una función de autocorrelación?

Solución propuesta:

T=2, es función rectangular es el ancho de la base.

S_Y(f) = F\Big[ \prod \Big(\frac{\tau}{T} \Big) \Big] = 2AT \frac {Sen(2\pi f \tau)}{2\pi fT} = AT Sa (\pi f\tau)

La función S_X(f) es negativa para algunos rangos de f. Dado que la densidad espectral de potencia es no negativa, la función rectangular en el tiempo no es una función de autocorrelacón válida.

---

2. PSD Triangular

 Referencia: Problema León García 10.1 p.635

Sea g(x) una función triangular.

amplitud = A
T = 1 , en triangulares T es la mitad de la base del triángulo.

a) Encuentre la densidad espectral de potencia correspondiente a R_X(τ) = g(τ/T)

b) Encuentre la autocorrelación correspondiente a la densidad espectral de potencia S_X(f) = g(f/W)

Solución propuesta:
a)

S_x(f) = F\Big[ g \Big(\frac{\tau}{T} \Big) \Big] = AT \Big(\frac{sin\frac{\omega T}{2}}{\frac{\omega T}{2}} \Big)^2 = AT \big[Sa(\pi f T) \big]^2

b)

S_x(f) = F\Big[ g \Big(\frac{f}{W} \Big) \Big] R_X (\tau) = AW \Big(\frac{sin\frac{W \tau}{2}}{\frac{W \tau}{2}} \Big)^2 = AW \big[Sa(\pi f \tau) \big]^2

 

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.

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

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 R_YX(τ), S_YX(f), R_Y(τ) y S_Y(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 S_X(f) = S_Y (f) sin importar el hecho que X(t) ≠ Y(t). Entonces S_X(f) = S_Y(f) lo que no implica que X(t) = Y(t)

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 S_X(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 R_X(0) = a²/2.
De toda esta potencia, se concentra en las frecuencias en f₀ positiva y negativa, por lo que la densidad espectral de potencia en esta frecuencias es infinita.