Meu primeiro programa Python - computação científica com Python

Depois de uma resistência inicial muito grande, as evidências mostraram que eu devo migrar de Scilab/Matlab para Python. Não vai ser muito trivial, mas vamos tentar. Nesse sentido, o código abaixo é nossa primeira experiência de fazer algo útil em Python. Como sempre, código com poucos comentários. Após rodar o código, o resultado final é o gráfico acima.

---

import numpy as np
import matplotlib.pyplot as plt
import math as mt

N = 100000;
#### Gerando simbolos +1, -1 e o ruido gaussiano normalizado: 

x = np.random.randn(N,1);
a = np.sign(x);
Pa = sum(a*a)/N; Pa = Pa[0];
n = np.random.randn(N,1);
n = n - np.mean(n);
Pn = sum(n*n)/N; Pn = Pn[0]; sn = np.sqrt(Pn);
n = n/sn;
Pn = sum(n*n)/N; Pn = Pn[0];
print('Conferindo a potência do ruído:', Pn);

#### SNR e taxa de erro:

k = 1;
tx =  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0];
txt = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0];
vdb = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11];
while k<13:
    A = np.sqrt(1/(10**(vdb[k-1]/10)));  # amplitude do ruido     

    nn = A*n;

    ### sinal + ruido:     

    y = a + nn;

    #### Erros:     

    ad = np.sign(y);
    erros = ad - a;
    n_erros = 0.5*np.sum(np.abs(erros));
    tx_erros = n_erros/N;
    z = A*np.sqrt(2);
    Pett = 0.5*mt.erfc(1/z);

    txt[k - 1] = Pett;
    tx[k - 1] = tx_erros;  # np.log10(tx_erros);     

    print(vdb[k - 1], tx[k - 1]);
    k = k + 1;

### Gráfico:

plt.semilogy(vdb, tx,'-o');
plt.semilogy(vdb,txt,'m');
plt.title('Taxa de erro de bit - sinal BPSK');
plt.xlabel('$SNR - dB$');
plt.ylabel('$BER$');
plt.legend(('Simulação','Analítico'));
plt.grid(True);
plt.show();