De la web www.matpic.com: agosto 2015

Errata del libro Señales y Sistemas de Oppenheim.

Ecuación 1.39 es exp(-j0.5t) (página 20).
Ecuación 1.123 es x[n].
Ecuación 2.14 (página 86).
Ecuación 2.21 (página 89).
Ecuación 2.30 el h lleva sombrero (página 96).
Ecuación 3.61 falta 't' en el exponente.
Ecuación 3.124. El argumentos de H es jkwo. No con la exponencial.
Problema 3.8. La señal x(t) es real e impar.
Página 305. Debe ser mayúscula la x.
Página 310, ejemplo 4.13: la última ecuación es -1, no -2.
Ecuación 4.45 (página 313) fata 1 sobre 2Pi.
Ec 4.49 exponente w0 falta.
Luego de la ecuación 4.83 el exponente de la segunda exponencial
Ejemplo 5.2. El menos de la seg exponencial
Ejemplo 5.7 Entre paréntesis va: y otros múltiplos impares de Pi.
Ecuación 5.47 (página 380) el límite inferior es -Inf.
Ejemplo 5.9. Para 0, “si n es impar”.
Figura 5.17. No es omega 0, es omega c.
Ecuación 9.90 es aR.
Figura 9.24 es aR (ver etiqueta también), no R/a.
Página 687, luego de 9.94 “...poles at s=-1+-j3”.
Ecuación 7.18 es Xc mayúscula.
Página 916, ecuación A.53 no es A12, es A21.

Modulación AM.

#AM modulation
# Funcion SINC
import numpy as np
import pylab as plt
from scipy.fftpack import fft

N = 500 # number of sample points
dt = 1. / 1000 # sample spacing

npp=9
nN=N*npp
t = np.linspace(-N*dt, N*dt, nN)

y=np.sinc(150*t)*np.cos(2*np.pi*400*t)
#y=np.abs(np.sinc(t))
plt.figure(1)
plt.plot(t,y, linewidth=2)
plt.grid(True), plt.xlabel('Tiempo')
plt.title(u'x(t)=sinc(150t)*cos(2 $\pi$ (400)t)')
plt.xlim([-6./200,6./200])

# FFT
yf = fft(y)
tf = npp*np.linspace(-1./(4.*dt), 1./(4.*dt), nN)
spectrum = 1./nN * np.concatenate([np.abs(yf[nN/2:nN]), np.abs(yf[0:nN/2])])

#figure1 = plt.figure(4, (10, 5))
plt.figure(2)
plt.plot(tf, spectrum, linewidth=2)
plt.grid()
plt.xlim([-700,700])
plt.title(u'Espectro de magnitud |X(j$\omega$)|')
plt.xlabel('Frecuencia [Hz]')
plt.ylabel('Magnitud |X(j$\omega$)|')

Transformada de Fourier de la función Sinc.

import numpy as np
import pylab as plt
from scipy.fftpack import fft

N = 500 # number of sample points
dt = 1. / 1000 # sample spacing

npp=3
nN=N*npp
t = np.linspace(-N*dt, N*dt, nN)

y=np.sinc(200*t)
#y=np.abs(np.sinc(t))
plt.figure(1)
plt.plot(t,y)
plt.grid(True), plt.xlabel('Tiempo')
plt.title(u'Sinc(200t)')
plt.xlim([-6./200,6./200])

# FFT
yf = fft(y)
tf = npp*np.linspace(-1./(4.*dt), 1./(4.*dt), nN)
spectrum = 1./nN * np.concatenate([np.abs(yf[nN/2:nN]), np.abs(yf[0:nN/2])])

#figure1 = plt.figure(4, (10, 5))
plt.figure(2)
plt.plot(tf, spectrum, '-')
plt.grid()
plt.xlim([-300,300])
plt.title(u'Espectro de magnitud |X(j$\omega$)|')
plt.xlabel('Frecuencia [Hz]')
plt.ylabel('Magnitud |X(j$\omega$)|')

Espectro unilateral de un suma de senoides en Python.

import numpy as np
#from scipy import signal
from scipy.fftpack import fft
import matplotlib.pyplot as plt

N = 1000 # number of sample points
dt = 1. / 500 # sample spacing

frequency1 = 50.
frequency2 = 150.
nN=N*4
t = np.linspace(0.0, N*dt, nN)
s1 = 0.8*np.sin(2*np.pi * frequency1 * t)
s2 = 0.4* np.sin(2*np.pi * frequency2 * t)
y = s1 + s2

plt.figure(1)
plt.plot(t, s1)
plt.grid()
plt.title('Senoide de '+str(frequency1)+' Hz')
plt.ylabel('Amplitud')
plt.xlabel('Tiempo [s]')
plt.axis([0, 8/frequency2, -1.5, 1.5])

plt.figure(2)
plt.plot(t, s2)
plt.grid()
plt.title('Senoide de '+str(frequency2)+' Hz')
plt.ylabel('Amplitud')
plt.xlabel('Tiempo [s]')
plt.axis([0, 8/frequency2, -1.5, 1.5])


plt.figure(3)
plt.plot(t, s1+s2)
plt.grid()
plt.title('Suma de senoides de '+str(frequency1)+' Hz'+' y '+str(frequency2)+' Hz')
plt.ylabel('Amplitud')
plt.xlabel('Tiempo [s]')
plt.axis([0, 8/frequency2, -1.5, 1.5])
# FFT
yf = fft(y)
tf = np.linspace(.0, 1./(2.*dt), N/2)
spectrum = 2./nN * np.abs(yf[0:N/2])

#figure1 = plt.figure(4, (10, 5))
plt.figure(4)
plt.plot(tf, spectrum, '-')
plt.grid()
plt.title(u'Espectro de magnitud |X(j$\omega$)|')
plt.xlabel('Frequencia [Hz]')
plt.ylabel('Magnitud |X(j$\omega$)|')

Espectro bilateral de un suma de senoides en Python.

import numpy as np
#from scipy import signal
from scipy.fftpack import fft,fftshift
import matplotlib.pyplot as plt

N = 1000 # number of sample points
Fs=500.
dt = 1. / Fs # sample spacing

frequency1 = 50.
frequency2 = 150.
npp=4
nN=N*npp
t = np.linspace(0.0, N*dt, nN)
s1 = 0.8*np.sin(2*np.pi * frequency1 * t)
s2 = 0.4* np.sin(2*np.pi * frequency2 * t)
y = s1 + s2

plt.figure(1)
plt.plot(t, s1)
plt.grid()
plt.title('Senoide de '+str(frequency1)+' Hz')
plt.ylabel('Amplitud')
plt.xlabel('Tiempo [s]')
plt.axis([0, 8/frequency2, -1.5, 1.5])

plt.figure(2)
plt.plot(t, s2)
plt.grid()
plt.title('Senoide de '+str(frequency2)+' Hz')
plt.ylabel('Amplitud')
plt.xlabel('Tiempo [s]')
plt.axis([0, 8/frequency2, -1.5, 1.5])


plt.figure(3)
plt.plot(t, s1+s2)
plt.grid()
plt.title('Suma de senoides de '+str(frequency1)+' Hz'+' y '+str(frequency2)+' Hz')
plt.ylabel('Amplitud')
plt.xlabel('Tiempo [s]')
plt.axis([0, 8/frequency2, -1.5, 1.5])
# FFT
yf = fft(y)
tf = npp*np.linspace( -1./(2.*dt), 1./(2.*dt), nN)
spectrum = 2./nN * np.abs(fftshift(yf))

#figure1 = plt.figure(4, (10, 5))
plt.figure(4)
plt.plot(tf, spectrum, '-')
plt.grid()
plt.title(u'Espectro de magnitud |X(j$\omega$)|')
plt.xlabel('Frequencia [Hz]')
plt.ylabel('Magnitud |X(j$\omega$)|')
#plt.axis([-200, 200, 0., 1.5])

Coeficientes de la serie de Fourier conforme T tiende al infinito.

import numpy as np
import pylab as pl
 
T1=1
T=4*T1
w0=2*np.pi/T
k=np.arange(-40*w0,40*w0,w0)
ak=(2.0*T1)*np.sinc(k)
pl.subplot(4,1,1)
pl.stem(k,ak*T)
pl.hold(True)
ke=np.arange(-40*w0,40*w0,w0/10)
ake=(2.0*T1)*np.sinc(ke)
pl.plot(ke,ake*T,'r--')
pl.hold(False)
pl.grid(True)
pl.xlim([-15,15])
 
T=8*T1
w0=2*np.pi/T
k=np.arange(-40*w0,40*w0,w0)
ak=(2.0*T1)*np.sinc(k)
pl.subplot(4,1,2)
pl.stem(k,ak*T)
pl.hold(True)
ke=np.arange(-40*w0,40*w0,w0/10)
ake=(2.0*T1)*np.sinc(ke)
pl.plot(ke,ake*T,'r--')
pl.hold(False)
pl.grid(True)
pl.xlim([-15,15])
 
T=16*T1
w0=2*np.pi/T
k=np.arange(-40*w0,40*w0,w0)
ak=(2.0*T1)*np.sinc(k)
pl.subplot(4,1,3)
pl.stem(k,ak*T)
pl.hold(True)
ke=np.arange(-40*w0,40*w0,w0/10)
ake=(2.0*T1)*np.sinc(ke)
pl.plot(ke,ake*T,'r--')
pl.hold(False)
pl.grid(True)
pl.xlim([-15,15])
 
T=32*T1
w0=2*np.pi/T
k=np.arange(-80*w0,80*w0,w0)
ak=(2.0*T1)*np.sinc(k)
pl.subplot(4,1,4)
pl.stem(k,ak*T)
pl.hold(True)
ke=np.arange(-40*w0,40*w0,w0/10)
ake=(2.0*T1)*np.sinc(ke)
pl.plot(ke,ake*T,'r--')
pl.hold(False)
pl.grid(True)
pl.xlim([-15,15])

Función sinc en Python.

# Funcion SINC
import numpy as np
import pylab as pl

t=np.linspace(-10,10,1000)
s=np.sinc(t)
s_mod=np.abs(np.sinc(t))
pl.figure(1)
pl.plot(t,s)
pl.grid(True), pl.xlabel('Tiempo')

pl.figure(2)
pl.plot(t,s_mod)
pl.grid(True), pl.xlabel('Tiempo')

Serie de Fourier de una onda triangular.

import pylab
import numpy as np
#Reconstruccion
A=1.
a0=1./2
t=np.linspace(-3,3,10000)
s_per=0
for k in range(-30,30):     
    #ak=((-1)^k-1)*A/(np.pi^2*k^2)
    if k!=0:
        s_per=s_per+((-1)**k-1)*A/(np.pi**2*k**2)*np.exp(1j*k*2*np.pi*t)
s_per=s_per+a0    
pylab.plot(t,s_per), pylab.show()
pylab.grid(True)

Serie de Fourier de una onda cuadrada.

import pylab
import numpy as np
#Con SINC
k=np.arange(-30,30)
ak=(1/8.0)*np.sinc(k/8.0)*np.exp(-1j*k*np.pi/8.0)
m=np.abs(ak)
f=np.angle(ak)
#pylab.figure(2)
#pylab.stem(k,m,markerfmt=" ")
#pylab.show()
#Reconstruccion
t=np.linspace(-2,2,10000)
s_per=0
for k in range(-30,30):     
    #print s_per
    s_per=s_per+(1/8.0)*np.sinc(k/8.0)*np.exp(-1j*k*np.pi/8.0)*np.exp(1j*k*2*np.pi*t)    
pylab.plot(t,s_per), pylab.show()

Convolución continua con Python.

#By Mark Wickert from Signals and Systems For Dummies
def expo_conv(t,A,alpha):
    import numpy as np
    y = np.zeros(len(t))
    for k, tk in enumerate(t):
        if tk >= 0:
            y[k] = A/alpha*(1 - np.exp(-alpha*tk))
    return y

import numpy as np
import pylab as pl
import ssd
    
t = np.arange(-4,14,.01)
xc2 = ssd.step(t)
hc2 = ssd.step(t)*np.exp(-1*t)
yc2_num,tyc2 = ssd.conv_integral(xc2,t,hc2,t,('r','r'))
#pl.subplot(211)
pl.plot(tyc2,yc2_num,linewidth=2)
pl.grid(True)
pl.xlabel('Tiempo (t)')
pl.title('y(t)=x(t)*h(t)')
pl.axis([-8,10,-0.5,1.5])
#pl.subplot(212)
#pl.plot(t,expo_conv(t,1,1))

Convolución discreta con Python.

#By Mark Wickert from Signals and Systems For Dummies
import numpy as np
import pylab as pl
import ssd

n=np.arange(-4,6)
xd1=2*ssd.drect(n,4)
hd1=1.5*ssd.dimpulse(n)-0.5*ssd.drect(n+1,3)
#hd1=-0.5*ssd.dimpulse(n+1)+1*ssd.dimpulse(n)-0.5*ssd.dimpulse(n-1)
yd1_num,nd1=ssd.conv_sum(xd1,n,hd1,n)
#pl.subplot(1,3,1)
pl.figure(1)
pl.stem(n,xd1)
pl.title('x[n]'), pl.axis([-4,5,-1,3]),pl.grid(True)
#pl.subplot(1,3,2)
pl.figure(2)
pl.stem(n,hd1),pl.grid(True)
pl.title('h[n]')
#pl.subplot(1,3,3)
pl.figure(3)
pl.stem(nd1,yd1_num), pl.axis([-4,5,-2,2]),pl.grid(True)
pl.title('y[n]=x[n]*h[n]')

Descomposición de una secuencia discreta en una combinación lineal de impulsos desplazados ponderados.

import numpy as np
import pylab as pl
from impulse_step_ramp import *

#-----------------------------------------------------
def udelta(n):
    y=(n==0)*1 #Convert bool to int [...*1]
    #i=np.where(y==1)
    return y
#-----------------------------------------------------

n=np.arange(-10,10)
s=udelta(n+6)-3*udelta(n+5)-4*udelta(n+4)+1*udelta(n+3)-5*udelta(n+2)+7*udelta(n+1)-1*udelta(n)+\
3*udelta(n-1)+3*udelta(n-2)-5*udelta(n-3)+2*udelta(n-4)-1*udelta(n-5)-1*udelta(n-6)
pl.stem(n,s)
pl.axis([-10,10,-10,10])
pl.grid(True)
pl.xlabel('Tiempo [n]')
for i, txt in enumerate(n):
    pl.annotate(s[i], (n[i],s[i]+0.5))
pl.show

Senoidal, exponencial compleja y señal amortiguada con Python.

Señal coseno.

import numpy as np
import pylab as pl

t=np.linspace(0,1,100)
f0=5
s=np.cos(2*np.pi*f0*t-np.pi/2)
pl.plot(t,s)
pl.grid('on')
pl.xlabel('Tiempo [s]')
pl.ylabel('Amplitud')
pl.title('Senoidal')
pl.show()

Exponencial continua real.

import numpy as np
import pylab as pl
n=np.arange(-15,15)
e=np.exp(-0.1*n)
pl.plot(n,e,'r')
pl.grid(True)
pl.show()

Relación de Euler para seno y coseno.

import numpy as np
import pylab as pl

t=np.linspace(0,0.5,100)
coseno_s=0.5*(np.exp(1j*4*np.pi*t)+np.exp(1j*4*np.pi*t))
seno_s=0.5*((np.exp(1j*4*np.pi*t)+np.exp(1j*4*np.pi*t)))/1j
pl.subplot(1,2,1)
pl.plot(t,coseno_s)
pl.grid('on')
pl.xlabel('Tiempo [s]')
pl.ylabel('Amplitud')
pl.title('Coseno')
pl.subplot(1,2,2)
pl.plot(t,seno_s)
pl.grid('on')
pl.xlabel('Tiempo [s]')
#pl.ylabel('Amplitud')
pl.title('Seno')
pl.show()

Seno amortiguado.

c=np.sin(2*np.pi*n/15)
s_a=e*c
#pl.stem(n,s_a)
pl.plot(n,e,n,-e)
pl.plot(n,s_a,'r')
pl.grid(True)
pl.show()

Simulación de un sistema para transmisión de datos binarios en Python

Uno de los mejores libros que he revisado este año es Signals and Systems For Dummies de Mark Wickert. En una de sus aplicaciones titulada Simulate the System for Transmitting Binary Data in Python añadí una leve corrección al código en la parte de gráficas y llamado a funciones de la librería ssd.py.

import numpy as np
import pylab as pl
import ssd
from scipy import signal

r, b, data0 = ssd.BPSK_tx(100000,10,1.5,0,'src')
# 100000 symbols, Ns=10, Df = 1.5*Rb, 0dB down
r = ssd.cpx_AWGN(r,100,10) # EsN0=100dB, Ns=10
Pr,f = ssd.my_psd(r,2**10,Fs=10)
pl.figure(1)
pl.plot(f,10*np.log10(Pr))
z = signal.lfilter(b,1,r)# b is the SRC filter
#pl.figure(2)
ssd.eye_plot(np.real(z[2000:6000]),20)#20 samp wind
yI,yQ = ssd.scatter(z[2000:6000],10,0)
pl.figure(3)
pl.plot(yI,yQ,'.')
pl.axis('equal')
r, b, data0 = ssd.BPSK_tx(100000,10,1.5,0,'src')
r = ssd.cpx_AWGN(r,20,10) # EsN0=20dB, Ns=10
Pr,f = ssd.my_psd(r,2**10,Fs=10)
pl.figure(4)
pl.plot(f,10*np.log10(Pr))
z = signal.lfilter(b,1,r)
#pl.figure(5)
ssd.eye_plot(np.real(z[2000:6000]),20)#20 samp wind
#ssd.scatter_plot(z[2000:6000]*np.exp(1j*np.pi/5),10,0)
yI,yQ = ssd.scatter(z[2000:6000]*np.exp(1j*np.pi/5),10,0)
pl.figure(6)
pl.plot(yI,yQ,'.')
pl.axis('equal')