quarta-feira, 28 de dezembro de 2016

Meu segundo programa Python - números complexos e equalização



Neste segundo programa Python, apresentamos um código para equalização de um sinal QAM após passar por um canal complexo. Não conseguimos fazer o comando "numpy.convolve" funcionar adequadamente e criamos um código (pouco eficiente) para realizar a convolução. Os gráficos podem ser vistos acima. Sobre equalização - ver aqui no blog. Segue o código:


import numpy as np
import matplotlib.pyplot as plt

plt.close('all');

N = 1000;  # número de símbolos
#### Gerando simbolos +1, -1 e o ruido guassiano normalizado:
x1 = np.random.randn(N,1);
x2 = np.random.randn(N,1);
a = np.sign(x1) + complex(0,1)*np.sign(x2);
n = np.random.randn(N,1) + complex(0,1)*np.random.randn(N,1);
n = 0.01*n;  # ruído
plt.figure();
plt.plot(np.real(a+n/10),np.imag(a+n/10),'r.'); plt.grid();
plt.title('Sinal 4QAM');

# canal complexo:
h = [complex(0.85,-0.30), complex(0.5,0.30), complex(0.15,0.10)]; h0 = h[0]; h1 = h[1]; h2 = h[2];


# convolução:
ac = complex(0,0)*np.ones((N,1));
ac[0] = a[0]*h0; 
ac[1] = a[1]*h0 + a[0]*h1; 
ac[2] = a[2]*h0 + a[1]*h1 + a[0]*h2; 

k=3;
while k    a0 = a[k]; a1 = a[k-1]; a2 = a[k-2];
    ac[k] = a0*h0 + a1*h1 + a2*h2;
    k = k + 1;

# ruido após o canal:
an = ac + n;
plt.figure();
plt.plot(np.real(an),np.imag(an),'r.'); plt.grid();
plt.title('Sinal 4QAM após o canal');

########## Equalizador:
M = 7;
w =  complex(0,0)*np.ones((M,1));
x =  complex(0,0)*np.ones((M,1));
y =  complex(0,0)*np.ones((N,1));
e =  complex(0,0)*np.ones((N,1));
mi = 0.25;
k = 0;
while k    
    sp = 0;
    n = 0;
    while n < M        
            if (k-n)>-1:
            sp = sp + an[k-n]*w[n];
            x[n] = an[k-n];
        n = n + 1;
    y[k] = sp;
    e[k] = a[k] - y[k];
    vx2 = np.sum(np.abs(x)*np.abs(x)) + 0.01;
    w = w + mi*(e[k])*np.conj(x)/vx2;
    k = k + 1;

# Gráficos:
plt.figure();
plt.plot(np.real(y[N-500:N]),np.imag(y[N-500:N]),'r.'); plt.title('Sinal 4QAM equalizado'); plt.grid();

Pe = np.sum(np.abs(e[N-500:N])*np.abs(e[N-500:N]))/(501);
plt.figure();
plt.plot(np.abs(e)); plt.grid();
plt.title('Curva de Aprendizado - Erro. Pe = '+str(Pe)); plt.show();

Nenhum comentário:

Postar um comentário