5.2 Transformada de Fourier - Señales aperiódicas continuas con Sympy-Python

Referencia: Oppenheim 4.1. p285, Lathi 7.2 p680

El planteamiento de Fourier es que una señal aperiódica (no periódica) puede observarse como una señal periódica con periodo infinito.

1. Ejercicio. exponencial decreciente con μ(t)

Referencia: Oppenheim Ejercicio 4.1 p290, Lathi ejemplo 7.1 p685, Hsu Ejemplo 5.2 p218

Considere la señal contínua exponencial decreciente, desarrolle la transformada de Fourier

x(t) =e^{-at} \mu (t) \text{ ; } a \gt 0

el integral a desarrollar,

X(\omega) =\int_{-\infty}^{\infty} e^{-at} \mu (t) e^{-j \omega t} \delta t

X(\omega) =\int_{0}^{\infty} e^{-(a+j\omega)t} \delta t =-\frac{1}{a+j\omega} e^{-(a+j\omega)t} \Big|_{0}^{\infty}

X(\omega) =-\frac{1}{a+j\omega} e^{-(a+j\omega)(\infty)} +\frac{1}{a+j\omega} e^{-(a+j\omega)(0)}

X(\omega) = \frac{1}{a+j\omega}

para a>0

Se obtiene una parte real y otra imaginaria al evaluar F(w) en un intervalo que incluya -a y a.

Transformada Fourier Aperiódica Ej01 Fw Re Imag

También es posible observar la magnitud y fase de F(w) en una gráfica.

Transformada Fourier Aperiódica Ej01 Fw Mg fase

Se observa por ejemplo que el valor de la magnitud |F(0)| = 1/a. También que la fase se encuentra acotada en π/2 y que fase(F(-a)) = -π/4.

Usando el algoritmo con Python se obtiene el siguiente resultado:

 expresion f(t):
 -a*t             
e    *Heaviside(t)

 expresion F(w):
   1   
-------
a + I*w

 |F(w)|:
     1      
------------
   _________
  /  2    2 
\/  a  + w  

 F(w)_fase:
     /w\
-atan|-|
     \a/
>>>

1.1 Algoritmo en Python

# Transformada de Fourier de señales aperiodicas
# blog.espol.edu.ec/telg1001/
import sympy as sym

# INGRESO
t = sym.Symbol('t', real=True,)
w = sym.Symbol('w', real=True)
a = sym.Symbol('a', real=True,positive=True)
T = sym.Symbol('T', real=True,positive=True)
j = sym.I
u = sym.Heaviside(t)
d = sym.DiracDelta(t)

ft = sym.exp(-a*t)*u
#ft = sym.exp(-a*sym.Abs(t))
#ft = sym.Heaviside(t+T)-sym.Heaviside(t-T)
#ft = d

# PROCEDIMIENTO
unilateral = False
ftw = ft*sym.exp(-j*w*t) # f(t,w) para integrar
ftw = sym.expand(ftw)  # expresion de sumas
ftw = sym.powsimp(ftw) # simplifica exponentes

lim_a = 0 # unilateral
if not(unilateral):
    lim_a = -sym.oo
    
# integral de Fourier
Fw_F = sym.integrate(ftw,(t,lim_a,sym.oo))
if Fw_F.is_Piecewise:
    Fw = Fw_F.args[0] # primera ecuacion e intervalo
    Fw = Fw[0] # solo expresion
else:
    Fw = Fw_F
Fw = Fw.simplify()

# Magnitud y Fase
Fw_magn = sym.Abs(Fw)
Fw_fase = sym.atan(sym.im(Fw)/sym.re(Fw))

# SALIDA
print('\n expresion f(t):')
sym.pprint(ft)
print('\n expresion F(w):')
sym.pprint(Fw)
print('\n |F(w)|:')
sym.pprint(Fw_magn)
print('\n F(w)_fase:')
sym.pprint(Fw_fase)

1.2 Algoritmo con Sympy

La librería Sympy incorpora una función para determinar en forma simbólica la transformada de Fourier, se usa la instrucción para simplificar el algoritmo anterior.

# Transformada de Fourier de señales aperiodicas
# blog.espol.edu.ec/telg1001/
import sympy as sym

# INGRESO
t = sym.Symbol('t', real=True,)
w = sym.Symbol('w', real=True)
a = sym.Symbol('a', real=True,positive=True)
T = sym.Symbol('T', real=True,positive=True)
j = sym.I
u = sym.Heaviside(t)
d = sym.DiracDelta(t)

ft = sym.exp(-a*t)*u
#ft = sym.exp(-a*sym.Abs(t))
#ft = u.subs(t,t+T)-u.subs(t,t-T)
#ft = d

# PROCEDIMIENTO
ft = sym.expand(ft)  # expresion de sumas
Fw = sym.fourier_transform(ft,t,w/(2*sym.pi))

# Magnitud y Fase
Fw_magn = sym.Abs(Fw)
Fw_fase = sym.atan(sym.im(Fw)/sym.re(Fw))

# SALIDA
print('\n expresion f(t):')
sym.pprint(ft)
print('\n expresion F(w):')
sym.pprint(Fw)
print('\n |F(w)|:')
sym.pprint(Fw_magn)
print('\n F(w)_fase:')
sym.pprint(Fw_fase)

1.3 Gráfica en Python

Para las gráficas, se dan valores a las constantes a, T, se define el intervalo de la gráfica y el número de muestras a usar.

Luego se sustituye en las ecuaciones resultantes los valores de ‘a‘ con ‘a_k‘ obteniendo valores para las gráficas. La gráfica para f(t) se realiza de forma semejante en la unidad 1. Para la gráfica F(w) se usan una gráfica con parte Real e Imaginaria, otra gráfica para magnitud y fase.

# GRAFICA ----------------
import numpy as np
import matplotlib.pyplot as plt
import telg1001 as fcnm
equivalentes = [{'DiracDelta': lambda x: 1*(x==0)},
                {'Heaviside': lambda x,y: np.heaviside(x, 1)},
                'numpy',]
a_k = 2 # valor de 'a' constante
# Grafica, intervalo tiempo [t_a,t_b]
T_k = 1    # valor T o periodo
t_a = -2*T_k ; t_b = 2*T_k
muestras = 52     # 51 resolucion grafica

ft = ft.subs({a:a_k,T:T_k}) # a tiene valor a_k
Fw = Fw.subs({a:a_k,T:T_k})
Fw_magn = Fw_magn.subs({a:a_k,T:T_k})
Fw_fase = Fw_fase.subs({a:a_k,T:T_k})

figura_ft = fcnm.graficar_ft(ft,t_a,t_b,muestras)

# F(w) real e imaginaria
w_a = -a_k*4 ; w_b = a_k*4
wi  = np.linspace(w_a,w_b,muestras)

# convierte a sympy una constante
Fw = sym.sympify(Fw,w) 
if Fw.has(w): # no es constante
    F_w = sym.lambdify(w,Fw,
                       modules=equivalentes)
else:
    F_w = lambda w: Fw + 0*w
Fwi = F_w(wi) # evalua wi

# F(w) magnitud y fase
# convierte a sympy una constante
Fw_magn = sym.sympify(Fw_magn,w)
if Fw_magn.has(w): # no es constante
    F_w_magn = sym.lambdify(w,Fw_magn,
                            modules=equivalentes)
else:
    F_w_magn = lambda w: Fw_magn + 0*w

# convierte a sympy una constante
Fw_fase = sym.sympify(Fw_fase,w) 
if Fw_fase.has(w): # no es constante
    F_w_fase = sym.lambdify(w,Fw_fase,
                            modules=equivalentes)
else:
    F_w_fase = lambda w: Fw_fase + 0*w

Fwi_magn = F_w_magn(wi) # evalua wi
Fwi_fase = F_w_fase(wi) 

# F(w) real e imaginaria
fig_Fw, graf_Fwi = plt.subplots(2,1)

graf_Fwi[0].plot(wi,np.real(Fwi),label='Re(F(w))',
                 color='orange')
graf_Fwi[0].legend()
graf_Fwi[0].set_ylabel('Re (F(w)) ')
graf_Fwi[0].grid()

graf_Fwi[1].plot(wi,np.imag(Fwi),label='Imag(F(w))',
                 color='brown')
graf_Fwi[1].legend()
graf_Fwi[1].set_xlabel('w')
graf_Fwi[1].set_ylabel('imag(F(w))')
graf_Fwi[1].grid()

plt.suptitle(r'F(w) = $'+str(sym.latex(Fw))+'$')
# plt.show()

# F(w) magnitud y fase
fig_Fwmg, graf_Fw = plt.subplots(2,1)

graf_Fw[0].plot(wi,Fwi_magn,label='F(w)_magn',
                color='orange')
Fwa0 = F_w_magn(0)
Fwa1 = F_w_magn(-a_k)
Fwa2 = F_w_magn(a_k)
if Fw_magn.has(w): # no es constante
    graf_Fw[0].stem([-a_k,a_k],[Fwa1,Fwa2],
                basefmt=' ',linefmt ='--')
    etiqueta1 = '('+str(a_k)+','+str(round(Fwa2,4))+')'
    graf_Fw[0].annotate(etiqueta1, xy=(a_k,Fwa2))
etiqueta0 = '('+str(0)+','+str(round(Fwa0,4))+')'
graf_Fw[0].scatter(0,Fwa0)
graf_Fw[0].annotate(etiqueta0, xy=(0,Fwa0))
graf_Fw[0].legend()
graf_Fw[0].set_ylabel('F(w) magnitud ')
graf_Fw[0].grid()

graf_Fw[1].plot(wi,Fwi_fase,label='F(w)_fase',
                 color='brown')
graf_Fw[1].axhline(np.pi/2,linestyle ='--')
graf_Fw[1].axhline(-np.pi/2,linestyle ='--')
if Fw_magn.has(w): # no es constante
    Fwa1f = F_w_fase(-a_k)
    Fwa2f = F_w_fase(a_k)
    graf_Fw[1].stem([-a_k,a_k],[Fwa1f,Fwa2f],
                basefmt=' ',linefmt ='--')
    etiqueta3 = '('+str(a_k)+','+str(round(Fwa2f,4))+')'
    graf_Fw[1].annotate(etiqueta3, xy=(a_k,Fwa2f))
graf_Fw[1].legend()
graf_Fw[1].set_xlabel('w')
graf_Fw[1].set_ylabel('F(w) fase')
graf_Fw[1].grid()

plt.suptitle(r'F(w) = $'+str(sym.latex(Fw))+'$')

plt.show()

2. Ejercicio. Exponencial decreciente con |t|, función par

Referencia: Oppenheim Ejercicio 4.2 p291, Lathi ejemplo 7. p685, Hsu 5.21 p248

Considere la señal contínua exponencial decreciente, desarrolle la transformada de Fourier

x(t) =e^{-a|t|} \text{ ; } a \gt 0

X(\omega) = \int_{-\infty}^{\infty}e^{-a|t|} e^{-j \omega t} \delta t

X(\omega) = \int_{-\infty}^{0}e^{at} e^{-j \omega t} \delta t + \int_{0}^{\infty}e^{-at} e^{-j \omega t} \delta t

= \int_{-\infty}^{0}e^{at-j \omega t} \delta t + \int_{0}^{\infty}e^{-at-j \omega t} \delta t

= \frac{1}{at-j \omega} e^{(a-j \omega) t}\Big|_{-\infty}^{0} - \frac{1}{a+j\omega}e^{-(at+j \omega) t} \Big| _{0}^{\infty}

= \Big[ \frac{1}{at-j \omega} e^{(a-j \omega) (0)} - \frac{1}{at-j \omega t} e^{(a-j \omega) (-\infty)} \Big] +

- \Big[ \frac{1}{a+j\omega}e^{-(at+j \omega) (\infty)} - \frac{1}{a+j\omega}e^{-(at+j \omega)(0)}\Big]

= \frac{1}{a-j\omega} +\frac{1}{a+j\omega} = \frac{a+j\omega+ a-j\omega}{(a-j\omega)(a+j\omega)}

X(\omega) = \frac{2a}{(a^2+\omega^2)}

Para desarrollar el ejercicio con el algoritmo, la entrada de señal se expresaría en el algoritmo como:

ft = sym.exp(-a*sym.Abs(t)) # Ej 4.2 Oppenheim

el resultado a comparar con el desarrollo analítico es:

 expresion f(t):
 -a*|t|
e      

 expresion F(w):
  2*a  
-------
 2    2
a  + w 

 |F(w)|:
  2*a  
-------
 2    2
a  + w 

 F(w)_fase:
0

Gráfica de F(w) magnitud y fase

Transformada Fourier Aperiódica Ej02 Fw Re Imag

El algoritmo es el mismo que el ejercicio 1, modificando el bloque de ingreso para el problema

3. Ejercicio. Rectangular centrada en origen

Referencia: Oppenheim Ejercicio 4.4 p293, Lathi ejemplo 7.2 p689, Hsu 5.19 p247

Considere la señal pulso rectangular o pulso compuerta (gate), desarrolle la transformada de Fourier

x(t) =\begin{cases}1 && |t|<T_1, \\ 0 && |t|>T_1\end{cases}

X(\omega) = \int_{-T_1}^{T_1} e^{-j\omega t} \delta t

= -\frac{1}{j \omega} e^{-j\omega t}\Big|_{-T_1}^{T_1} = -\frac{1}{j \omega} \Big[ e^{-j\omega (T_1)} - e^{-j\omega (-T_1)}\Big]

= -\frac{1}{j \omega} \Big[ e^{-T_1 j\omega } - e^{jT_1\omega } \Big]

= \frac{1}{j \omega} e^{T_1 j\omega} - \frac{1}{j \omega} e^{-T_1 j\omega }

en este punto es conveniente usar la forma trigonométrica de un exponencial con exponente complejo

= \frac{1}{j \omega} (\cos (T_1\omega)+j \sin(T_1\omega)) - \frac{1}{j \omega} (cos(T_1 \omega) -jsin(T_1\omega))

= \frac{1}{j \omega}\cos (T_1\omega)+j\frac{1}{j \omega} \sin(T_1\omega) - \frac{1}{j \omega} cos(T_1 \omega) +j\frac{1}{j \omega} sin(T_1\omega))

X(\omega) = 2\frac{\sin(T_1\omega)}{\omega}

Para desarrollar el ejercicio con el algoritmo, la señal se expresaría como:

ft = sym.Heaviside(t+T)-sym.Heaviside(t-T)
# Ej 4.4 Oppenheim

obteniendo el siguiente resultado

 expresion f(t):
-Heaviside(-T + t) + Heaviside(T + t)

 expresion F(w):
2*sin(T*w)
----------
    w     

 |F(w)|:
  |sin(T*w)|
2*|--------|
  |   w    |

 F(w)_fase:
0

con la siguiente gráfica f(T)

gráfica de F(w) parte real e imaginaria

Transformada Fourier Aperiódica Ej03 Fw Re Imag

El algoritmo es el mismo que el ejercicio 1, modificando el bloque de ingreso para el problema

4. Ejercicio. Pulso unitario

Referencia: Oppenheim Ejercicio 4.3 p292, Lathi ejemplo 7.3 p693, Hsu 5.1 p218

Considere la señal pulso unitario, desarrolle la transformada de Fourier

x(t) = \delta (t)

X(j\omega) = \int_{-\infty}^{\infty} \delta (t) e^{-j \omega t} \delta t = 1

Es decir la transformada de Fourier tiene componentes de todas las frecuencias.

Transformada Fourier Aperiódica Ej04 Fw Re Imag

El algoritmo entrega el siguiente resultado:

 expresion f(t):
DiracDelta(t)

 expresion F(w):
1

 |F(w)|:
1

 F(w)_fase:
0

El pulso unitario se define en Sympy como:

ft = sym.DiracDelta(t) # Ej 4.3 Oppenheim

En la parte gráfica, el pulsos unitarios se grafican con plt.stem(0,1), no requiriendo el resto de las gráficas para observar el resultado.

El algoritmo es el mismo que el ejercicio 1, modificando el bloque de ingreso para el problema

Tarea: Realizar otros ejercicios con exponenciales para comprobar la operación del algoritmo.

5.2 Transformada de Fourier - Señales aperiódicas continuas con Sympy-Python

1. Ejercicio. exponencial decreciente con μ(t)

1.1 Algoritmo en Python

1.2 Algoritmo con Sympy

1.3 Gráfica en Python

2. Ejercicio. Exponencial decreciente con |t|, función par

3. Ejercicio. Rectangular centrada en origen

4. Ejercicio. Pulso unitario

Señales y Sistemas

Signal & Systems