-
Estudiante: Juan Daniel Madrigal Solís
-
Carné: B64043
-
Estudiante: David Mairena Castro
-
Carné: B43969
-
Estudiante: Sofia Villalobos Brenes
-
Carné: B98464
-
Grupo: 4.12
La modulación digital es una de las aplicaciones del análisis de procesos estocásticos, y es parte de los sistemas digitales de comunicación. Este proyecto presenta una introdución a tópicos fundamentales de la ingeniería de comunicaciones para simular un sistema de transmisión de imágenes de baja resolución.
Las bibliotecas de Python de interés para este proyecto son:
# Para manipular imágenes (Python Imaging Library)
from PIL import Image
# Para manipular 'arrays' de pixeles y bits, señales y operaciones
import numpy as np
# Para visualizar imágenes y señales
import matplotlib.pyplot as plt
# Para medir el tiempo de simulación
import time
from PIL import Image
import numpy as np
def fuente_info(imagen):
'''Una función que simula una fuente de
información al importar una imagen y
retornar un vector de NumPy con las
dimensiones de la imagen, incluidos los
canales RGB: alto x largo x 3 canales
:param imagen: Una imagen en formato JPG
:return: un vector de pixeles
'''
img = Image.open(imagen)
return np.array(img)
import numpy as np
def rgb_a_bit(array_imagen):
'''Convierte los pixeles de base
decimal (de 0 a 255) a binaria
(de 00000000 a 11111111).
:param imagen: array de una imagen
:return: Un vector de (1 x k) bits 'int'
'''
# Obtener las dimensiones de la imagen
x, y, z = array_imagen.shape
# Número total de elementos (pixeles x canales)
n_elementos = x * y * z
# Convertir la imagen a un vector unidimensional de n_elementos
pixeles = np.reshape(array_imagen, n_elementos)
# Convertir los canales a base 2
bits = [format(pixel, '08b') for pixel in pixeles]
bits_Rx = np.array(list(''.join(bits)))
return bits_Rx.astype(int)
import numpy as np
def canal_ruidoso(senal_Tx, Pm, SNR):
'''Un bloque que simula un medio de trans-
misión no ideal (ruidoso) empleando ruido
AWGN. Pide por parámetro un vector con la
señal provieniente de un modulador y un
valor en decibelios para la relación señal
a ruido.
:param senal_Tx: El vector del modulador
:param Pm: Potencia de la señal modulada
:param SNR: Relación señal-a-ruido en dB
:return: La señal modulada al dejar el canal
'''
# Potencia del ruido generado por el canal
Pn = Pm / pow(10, SNR/10)
# Generando ruido auditivo blanco gaussiano (potencia = varianza)
ruido = np.random.normal(0, np.sqrt(Pn), senal_Tx.shape)
# Señal distorsionada por el canal ruidoso
senal_Rx = senal_Tx + ruido
return senal_Rx
import numpy as np
def bits_a_rgb(bits_Rx, dimensiones):
'''Un blque que decodifica el los bits
recuperados en el proceso de demodulación
:param: Un vector de bits 1 x k
:param dimensiones: Tupla con dimensiones de la img.
:return: Un array con los pixeles reconstruidos
'''
# Cantidad de bits
N = len(bits_Rx)
# Se reconstruyen los canales RGB
bits = np.split(bits_Rx, N / 8)
# Se decofican los canales:
canales = [int(''.join(map(str, canal)), 2) for canal in bits]
pixeles = np.reshape(canales, dimensiones)
return pixeles.astype(np.uint8)
La modulación 8-PSK tiene ocho símbolos posibles con tres bits b1 b2 b3 por símbolo.
El "diagrama de constelación" de la modulación es:
Ahí es posible verificar la correspondencia de los bits b1 b2 b3 (000, 001, 010, 011, 100, 101, 110, 111) con los puntos de la constelación.
- (50%) Realice una simulación del sistema de comunicaciones como en la sección 3.2., pero utilizando una modulación 8-PSK en lugar de una modulación BPSK. Deben mostrarse las imágenes enviadas y recuperadas y las formas de onda.
4.1.1 Se establecen las funciones necesarias
import numpy as np
def modulador_8_PSK(bits, fc, mpp):
'''Un método que simula el esquema de modulación digital 8-PSK.
:param bits: Vector unidimensional de bits
:param fc: Frecuencia de la portadora en Hz
:param mpp: Cantidad de muestras por periodo de onda portadora
:return: Un vector con la señal modulada
:return: Un valor con la potencia promedio [W]
:return: La onda portadora1 c1(t) = cos(2πfct)
:return: La onda portadora2 c2(t) = sin(2πfct)
'''
# 1. Parámetros de la 'señal' de información (bits)
N = len(bits) # Cantidad de bits
# 2. Construyendo un periodo de la señal portadora c(t)
Tc = 1 / fc # periodo [s]
t_periodo = np.linspace(0, Tc, mpp) # mpp: muestras por período
portadora1 = np.cos(2*np.pi*fc*t_periodo) # cos(2πfct)
portadora2 = np.sin(2*np.pi*fc*t_periodo) # sin(2πfct)
# 3. Inicializar la señal modulada s(t)
t_simulacion = np.linspace(0, N*Tc, N*mpp)
senal_Tx = np.zeros(t_simulacion.shape)
# 4. Asignar las formas de onda según los bits (8-PSK)
h = np.sqrt(2)/2
for i in range(0, N, 3):
if bits[i] == 1 and bits[i+1] == 1 and bits[i+2] == 1:
senal_Tx[i*mpp: (i+1)*mpp] = portadora1 * 1 + portadora2 * 0
elif bits[i] == 1 and bits[i+1] == 1 and bits[i+2] == 0:
senal_Tx[i*mpp: (i+1)*mpp] = portadora1 * h + portadora2 * h
elif bits[i] == 0 and bits[i+1] == 1 and bits[i+2] == 0:
senal_Tx[i*mpp: (i+1)*mpp] = portadora1 * 0 + portadora2 * 1
elif bits[i] == 0 and bits[i+1] == 1 and bits[i+2] == 1:
senal_Tx[i*mpp: (i+1)*mpp] = portadora1 * -h + portadora2 * h
elif bits[i] == 0 and bits[i+1] == 0 and bits[i+2] == 1:
senal_Tx[i*mpp: (i+1)*mpp] = portadora1 * -1 + portadora2 * 0
elif bits[i] == 0 and bits[i+1] == 0 and bits[i+2] == 0:
senal_Tx[i*mpp: (i+1)*mpp] = portadora1 * -h + portadora2 * -h
elif bits[i] == 1 and bits[i+1] == 0 and bits[i+2] == 0:
senal_Tx[i*mpp: (i+1)*mpp] = portadora1 * 0 + portadora2 * -1
else:
senal_Tx[i*mpp: (i+1)*mpp] = portadora1 * h + portadora2 * -h
# 5. Calcular la potencia promedio de la señal modulada
Pm = (1 / (N*Tc)) * np.trapz(pow(senal_Tx, 2), t_simulacion)
return senal_Tx, Pm, portadora1, portadora2
def demodulador_8_PSK(senal_Rx, portadora1, portadora2, mpp):
'''Un método que simula un bloque demodulador
de señales, bajo un esquema BPSK. El criterio
de demodulación se basa en decodificación por
detección de energía.
:param senal_Rx: La señal recibida del canal
:param portadora: La onda portadora c(t)
:param mpp: Número de muestras por periodo
:return: Los bits de la señal demodulada
'''
# Cantidad de muestras en senal_Rx
M = len(senal_Rx)
# Cantidad de bits (símbolos) en transmisión
N = int(M / mpp)
# Vector para bits obtenidos por la demodulación
bits_Rx = np.zeros(N)
# Vector para la señal demodulada
senal_demodulada = np.zeros(senal_Rx.shape)
# Pseudo-energía de un período de la portadora1
Es1 = np.sum(portadora1 * portadora1)
# Pseudo-energía de un período de la portadora2
Es2 = np.sum(portadora2 * portadora2)
h = np.sqrt(2)/2
# Demodulación
for i in range(N):
# Producto interno de dos funciones
producto1 = senal_Rx[i*mpp: (i+1)*mpp] * portadora1
Ep1 = np.sum(producto1)
producto2 = senal_Rx[i*mpp: (i+1)*mpp] * portadora2
Ep2 = np.sum(producto2)
senal_demodulada[i*mpp: (i+1)*mpp] = producto1 + producto2
# Criterio de decisión por detección de energía
if i % 3 == 0:
if (Ep1 >= (1+h)/2*Es1 and -h/2*Es2 <= Ep2 <= h/2*Es2):
bits_Rx[i] = 1
bits_Rx[i+1] = 1
bits_Rx[i+2] = 1
elif (h/2*Es1 <= Ep1 <= (1+h)/2*Es1 and
h/2*Es2 <= Ep2 <= (1+h)/2*Es2):
bits_Rx[i] = 1
bits_Rx[i+1] = 1
bits_Rx[i+2] = 0
elif (-h/2*Es1 <= Ep1 <= h/2*Es1 and Ep2 >= (1+h)/2*Es2):
bits_Rx[i] = 0
bits_Rx[i+1] = 1
bits_Rx[i+2] = 0
elif (-(1+h)/2*Es1 <= Ep1 <= -h/2*Es1 and
h/2*Es2 <= Ep2 <= (1+h)/2*Es2):
bits_Rx[i] = 0
bits_Rx[i+1] = 1
bits_Rx[i+2] = 1
elif (Ep1 <= -(1+h)/2*Es1 and -h/2*Es1 <= Ep2 <= h/2*Es2):
bits_Rx[i] = 0
bits_Rx[i+1] = 0
bits_Rx[i+2] = 1
elif (-(1+h)/2*Es1 <= Ep1 <= -h/2*Es1 and
-(1+h)/2*Es2 <= Ep2 <= -h/2*Es2):
bits_Rx[i] = 0
bits_Rx[i+1] = 0
bits_Rx[i+2] = 0
elif (-h/2*Es1 <= Ep1 <= h/2*Es1 and Ep2 <= -(1+h)/2*Es2):
bits_Rx[i] = 1
bits_Rx[i+1] = 0
bits_Rx[i+2] = 0
else:
bits_Rx[i] = 1
bits_Rx[i+1] = 0
bits_Rx[i+2] = 1
return bits_Rx.astype(int), senal_demodulada
import numpy as np
import matplotlib.pyplot as plt
import time
def Modulacion_8PSK(SNR):
# Parámetros
fc = 5000 # frecuencia de la portadora
mpp = 20 # muestras por periodo de la portadora
# Iniciar medición del tiempo de simulación
inicio = time.time()
# 1. Importar y convertir la imagen a trasmitir
imagen_Tx = fuente_info(requests.get(
'https://github.com/DavidMairena/Tema4/blob/main/arenal.jpg?raw=true',
stream=True).raw)
dimensiones = imagen_Tx.shape
# 2. Codificar los pixeles de la imagen
bits_Tx = rgb_a_bit(imagen_Tx)
# 3. Modular la cadena de bits usando el esquema BPSK
senal_Tx, Pm, portadora1, portadora2 = modulador_8_PSK(bits_Tx, fc, mpp)
# 4. Se transmite la señal modulada, por un canal ruidoso
senal_Rx = canal_ruidoso(senal_Tx, Pm, SNR)
# 5. Se desmodula la señal recibida del canal
bits_Rx, senal_demodulada = demodulador_8_PSK(senal_Rx, portadora1,
portadora2, mpp)
# 6. Se visualiza la imagen recibida
imagen_Rx = bits_a_rgb(bits_Rx, dimensiones)
Fig = plt.figure(figsize=(10, 6))
# Cálculo del tiempo de simulación
print('Duración de la simulación: ', time.time() - inicio)
# 7. Calcular número de errores
errores = sum(abs(bits_Tx - bits_Rx))
BER = errores/len(bits_Tx)
print('{} errores, para un BER de {:0.4f}.'.format(errores, BER))
# Mostrar imagen transmitida
ax = Fig.add_subplot(1, 2, 1)
imgplot = plt.imshow(imagen_Tx)
ax.set_title('Transmitido')
# Mostrar imagen recuperada
ax = Fig.add_subplot(1, 2, 2)
imgplot = plt.imshow(imagen_Rx)
ax.set_title('Recuperado')
Fig.tight_layout()
plt.imshow(imagen_Rx)
# Visualizar el cambio entre las señales
fig, (ax1, ax2, ax3) = plt.subplots(nrows=3, sharex=True, figsize=(14, 7))
# La señal modulada por (8-PSK)
ax1.plot(senal_Tx[0:600], color='g', lw=2)
ax1.set_ylabel('$s(t)$')
# La señal modulada al dejar el canal
ax2.plot(senal_Rx[0:600], color='b', lw=2)
ax2.set_ylabel('$s(t) + n(t)$')
# La señal demodulada
ax3.plot(senal_demodulada[0:600], color='m', lw=2)
ax3.set_ylabel('$b^{\prime}(t)$')
ax3.set_xlabel('$t$ / milisegundos')
fig.tight_layout()
plt.show()
return senal_Tx
4.1.2 Simulación para la modulación 8-PSK:
# Se realizan la prueba con diferentes relaciones señal-a-ruido del canal (SNR)
SNR1 = -5 # Mayor ruido que la señal, va a tener muchos errores
senal_Tx1 = Modulacion_8PSK(SNR1)
SNR2 = 5 # Una proporción entre la señal y el ruido
senal_Tx2 = Modulacion_8PSK(SNR2)
SNR3 = 15 # La señal es mucho clara que el ruido
# En este caso si la demodulación está bien hecha no debería haber ningún error
senal_Tx = Modulacion_8PSK(SNR3)
- (30%) Realice pruebas de estacionaridad y ergodicidad a la señal modulada
senal_Txy obtenga conclusiones sobre estas.
import numpy as np
import matplotlib.pyplot as plt
# Tiempo de muestreo
t_x = np.linspace(0, 0.01, 100)
R = [1, -1]
# Formas de la onda:
X_t = np.empty((4, len(t_x)))
# Nueva figura
plt.figure()
# Matriz con los posibles valores de cada función
for i in R:
x1 = i * np.cos(2 * (np.pi) * fc * t_x) + i * np.sin(2 * (np.pi) * fc * t_x)
x2 = -i * np.cos(2 * (np.pi) * fc * t_x) + i * np.sin(2 * (np.pi) * fc * t_x)
X_t[i, :] = x1
X_t[i+1, :] = x2
plt.plot(t_x, x1, lw=4, color='g')
plt.plot(t_x, x2, lw=4, color='c')
# Se determina un promedio de las 4 realizaciones para cada instante:
PR = [np.mean(X_t[:, i]) for i in range(len(t_x))]
plt.plot(t_x, PR, lw=6, color='k', label='Promedio de realizaciones')
# Se grafica el resultado teórico del valor esperado:
S = senal_Tx.astype('float')
RT = np.mean(S) * t_x
plt.plot(t_x, RT, '-.', lw=3, color='y', label='Valor teórico esperado')
# Se muestran las realizaciones, junto al promedio calculado y teórico:
plt.title('Realizaciones del proceso aleatorio $X(t)$')
plt.xlabel('$t$')
plt.ylabel('$x_i(t)$')
plt.legend()
plt.show()
Estacionaridad
La estacionaridad está relacionada con la repetición de los patrones de la onda a lo largo del tiempo. Es decir, su variancia, media y posición no tendrán cambios abruptos a lo largo de un tiempo indefinido. A paritr de la gráfica generada, notamos que tanto su promedio como las ondas generadas serán invariables y continuarán con el patrón observado.
Ergodicidad
Cuando hablamos de ergodicidad, estamos hablando de un efecto que ocurre en algunos sistemas mecánicos en donde el promedio de las realizaciones a cada instante es igual al promedio del sistema a largo plazo. Esto se demuestra en la gráfica generada en donde ambos se mantienen constante a traves del tiempo.
- (20%) Determine y grafique la densidad espectral de potencia para la señal modulada
senal_Tx.
from scipy import fft
# Transformada de Fourier
senal_f = fft(senal_Tx)
# Muestras de la señal
Nm = len(senal_Tx)
# Número de símbolos (198 x 89 x 8 x 3)
Ns = Nm // mpp
# Tiempo del símbolo = periodo de la onda portadora
Tc = 1 / fc
# Tiempo entre muestras (período de muestreo)
Tm = Tc / mpp
# Tiempo de la simulación
T = Ns * Tc
# Espacio de frecuencias
f = np.linspace(0.0, 1.0/(2.0*Tm), Nm//2)
# Gráfica
plt.figure(figsize=(18,10))
plt.plot(f, 2.0/Nm * np.power(np.abs(senal_f[0:Nm//2]), 2))
plt.xlim(0, 20000)
plt.grid()
plt.show()
