Lo que para señales y sistemas representa un componente importante en el análisis de sistemas LTI, al poder representar los sitemas por bloques en serie y paralelo semejante a lo realizado con la transformada de Laplace:
La convergencia de la transformada de Fourier se garantiza bajo ciertas condiciones. El análisis de Fourier para el estudio de sistemas LTI se restringe a los que sus respuestas al impulso tienen transformada de Fourier. Los detalles se mostrarán en cada ejercicio presentado.
Ejemplo 1. h(t) como respuesta a un impulso desplazado
Se puede ver que la exponencial con frecuencia ω0 + 2π es la misma que aquella con la frecuencia ω0.
Existe una situación diferente en el caso contínuo, en el cual las señales ejω0t son todas distintas para distintos valores de ω0.
En el caso discreto, éstas señales no son distintas, ya que la señal con frecuencia ω0 es idéntica a las señales con frecuencias ω0 ± 2π , ω0 ± 4π, … y las que le siguen. Por lo tanto, al considerar las exponenciales complejas, necesitamos tomar en cuenta solamente un intervalo de frecuencia de longitud 2π dentro del cual se escoge ω0.
ω0N = 2π m
o de otra forma
ω0/2π = m/N
Para el caso de m=1 y N=2
otra prueba con m=7 y N=4
con m=2 y N=1
El código en python para realizar las gráficas es:
# Señales discretas, y(t) a y[t]# propuesta:edelros@espol.edu.ecimport numpy as np
import math
import matplotlib.pyplot as plt
# Definir la funcion para el ejemplodefanalogica(f,t):
# función matemática CAMBIAR AQUI
y=np.cos(f*t)
return(y)
# Programa para graficar la función# INGRESO
rango = 4 # float(input('rangos en periodos de analogica:'))
m = 1 # float(input('frecuencia de analógica:'))
N = 2 # float(input('frecuencia digital:'))# PROCEDIMIENTO# grafica analogica
puntoscontinuos = 500
t0 = -rango*2*np.pi*m/2
tn = rango*2*np.pi*m/2
t = np.linspace(t0,tn,puntoscontinuos+1)
yanalog = analogica(m,t)
# grafica digital
deltaD = (2*np.pi*m)/(2*N)
muestreo = int((tn-t0)//deltaD +1)
td = np.linspace(t0,tn,muestreo)
ydigital = analogica(m,td)
# SALIDA - GRAFICA#Escala y para el grafico
margen = 0.1*np.max(yanalog)
ymax = np.max(yanalog)+margen
ymin = np.min(yanalog)-margen
plt.figure(1)
plt.suptitle('Señal Analogica vs digital')
plt.subplot(211) # grafica de 2x1 arriba
plt.plot(t,yanalog)
plt.axis((t0,tn,ymin,ymax))
plt.ylabel('y(t)')
plt.subplot(212) # grafica de 2x1 abajo
plt.plot(td,ydigital,'ro')
plt.axis((t0,tn,ymin,ymax))
plt.ylabel('Digital: y[t]')
plt.show()
donde los coeficientes de la ecuación se calculan como:
a_k = \frac{2}{T} \int_0^T f(t) \cos(k \omega_0 t) \delta t b_k = \frac{2}{T} \int_0^T f(t) \sin(k \omega_0 t) \delta t
Ejemplo
Referencia: Lathi 6.5 p617, Chapra 5ed Ejemplo 19.2 p548/pdf572
Utilice la serie de Fourier continua para aproximar la función de onda cuadrada o rectangular. Para el cálculo del ejemplo se usará hasta 4 términos de la serie.
Para determinar las expresiones de los coeficientes, en Python se usa la libreria Sympy que nos facilita el desarrollo las fórmulas simbólicas ak y bk.
Si requiere revisar el concepto se adjunta el enlace:
Usando Python se obtiene los siguientes resultados:
expresión ak:
/1.2732*sin(1.5707*k) - 0.6366*sin(3.1415*k) for And(k > -oo, k < oo, k != 0)
|--------------------------------------------
< k
|
\ 0 otherwise
expresión bk:
0
Nota: se han truncado los decimales a 4 para facilitar la lectura en la pantalla.
usando formato latex para la expresión, generado por Sympy se obtiene:
a_k = \begin{cases} \frac{1.2732\sin (1.5707 k) - 0.6366 \sin(3.1415 k )}{k} & \text{for}\: k > -\infty \wedge k < \infty \wedge k \neq 0 \\0 & \text{otherwise} \end{cases} b_k = 0
Instrucciones en Python
# Serie de Fourier, con n coeficientes# Ref.Chapra 5ed Ejemplo 19.2 p548/pdf572import numpy as np
import sympy as sym
import matplotlib.pyplot as plt
# INGRESO
T = 2*np.pi
t = sym.Symbol('t')
ft = sym.Piecewise((-1,t <-T/2),
(-1,t <-T/4),
( 1,t < T/4),
(-1,t < T/2),
(-1,True),)
# intervalo
a = -T/2
b = T/2
# número de coeficientes
n = 4
# PROCEDIMIENTO
k = sym.Symbol('k')
w0 = 2*np.pi/T
# Términos ak para coseno
enintegral = ft*sym.cos(k*w0*t)
yaintegrado = sym.integrate(enintegral,(t,a,b))
ak = (2/T)*yaintegrado
ak = sym.simplify(ak)
# Términos bk para seno
enintegral = ft*sym.sin(k*w0*t)
yaintegrado = sym.integrate(enintegral,(t,a,b))
bk = (2/T)*yaintegrado
bk = sym.simplify(bk)
print(' expresión ak:')
sym.pprint(ak)
print('\n ak formato latex')
print(sym.latex(ak))
print('\n expresión bk:')
sym.pprint(bk)
print('\n bk formato latex')
print(sym.latex(bk))
Evaluación de coeficientes
Con las expresiones obtenidas en el bloque anterior, se evalúan los n coeficientes ak y bk, substituyendo en las expresiones los valores en cada índice i y se obtiene como resultado:
obteniendo la siguiente expresión para la serie de Fourier como fs(t)
serie de Fourier fs(t):
1.27323954473516*cos(1.0*t)
+ 1.55926873300775e-16*cos(2.0*t)
- 0.424413181578388*cos(3.0*t)
Instrucciones en Python
Las instrucciones se añaden a continuación de los bloques anteriores,
# serie de Fourier
serieF = a0 + 0*t
i = 1
whilenot(i>=n):
serieF = serieF + aki[i]*sym.cos(i*w0*t)
serieF = serieF + bki[i]*sym.sin(i*w0*t)
i = i+1
# serie = sym.simplify(serie)print('\n serie de Fourier fs(t): ')
sym.pprint(serieF)
Gráficas de f(t) y Serie de Fourier
Para comparar la función f(t) con la aproximación en series de Fourier, se usa el intervalo [a,b] con una cantidad de muestras usando la instrucción np.linspace() y guardando el resultado en ti.
Para evaluar los puntos ti en cada expresión, por simplicidad se convierte la expresión de su forma simbólica a numérica lambda, usando sym.lambdify()
Instrucciones en Python
Las instrucciones para graficar las series de Fourier van a continuación de las anteriores,
# Para evaluación numérica
fsn = sym.lambdify(t,serieF)
ftn = sym.lambdify(t,ft)
# Evaluación para gráfica
muestras = 41
ti = np.linspace(a,b,muestras)
fi = ftn(ti)
fsi = fsn(ti)
# SALIDA# Grafica
plt.plot(ti,fi,label = 'f(t)')
etiqueta = 'coeficientes = '+ str(n)
plt.plot(ti,fsi,label = etiqueta)
plt.xlabel('ti')
plt.legend()
plt.title('Serie de Fourier')
plt.show()
Tarea: Realizar el ejercicio, aumentando el número de términos a 8
Espectro de Fourier de amplitud y fase
Referencia: Lathi 6.1 p599, Chapra 19.3 p551
Los espectros ofrecen información que no aparece en el dominio del tiempo. La gráfica de frecuencias ofrece una representación rápida de la estructura de armónicas, que son como las huellas dactilares que ayudan a caracterizar y entender la forma de una señal complicada.
Para mostrar este espectro de frecuencias se incrementó el número de terminos de la serie a n=11, para observar mejor la forma de la gráfica.
El algoritmo final como una integración de las partes presentadas se usa en la página siguiente con algunos ejemplos tradicionales de la transformada de Fourier.
Instrucciones Python integradas
Se incorporan mejoras en la presentación de resultados, al considerar valores casicero, dígitos a mostrar del los coeficientes, presentación de tablas de coeficientes.
# Serie de Fourier, con n coeficientes# Ref.Chapra 5ed Ejemplo 19.2 p548/pdf572import numpy as np
import sympy as sym
import matplotlib.pyplot as plt
# INGRESO
T = 2*np.pi
t = sym.Symbol('t')
ft = sym.Piecewise((-1,t <-T/2),
(-1,t <-T/4),
( 1,t < T/4),
(-1,t < T/2),
(-1,True),)
# intervalo
a = -T/2
b = T/2
# número de coeficientes
n = 11
# Evaluación para gráfica
muestras = 101
digitos = 4 # numpy set_print_options
casicero = 1e-15
# PROCEDIMIENTO
k = sym.Symbol('k')
w0 = 2*np.pi/T
# Términos ak para coseno
enintegral = ft*sym.cos(k*w0*t)
yaintegrado = sym.integrate(enintegral,(t,a,b))
ak = (2/T)*yaintegrado
ak = sym.simplify(ak)
# Términos bk para seno
enintegral = ft*sym.sin(k*w0*t)
yaintegrado = sym.integrate(enintegral,(t,a,b))
bk = (2/T)*yaintegrado
bk = sym.simplify(bk)
print(' expresión ak:')
sym.pprint(ak)
#print('\n ak formato latex')#print(sym.latex(ak))print('\n expresión bk:')
sym.pprint(bk)
#print('\n bk formato latex')#print(sym.latex(bk))# evalua los coeficientes
a0 = ak.subs(k,0)/2
b0 = bk.subs(k,0)
aki = [a0]
bki = [b0]
n_i = [0]
i = 1
whilenot(i>=n):
avalor = ak.subs(k,i)
bvalor = bk.subs(k,i)
aki.append(avalor)
bki.append(bvalor)
n_i.append(i)
i = i+1
##print('\n coeficientes ak: ')##print(aki)##print('\n coeficientes bk: ')##print(bki)# tabla de valores ak, Bk
tabla = np.concatenate([[np.array(n_i,dtype = int)],
[np.array(aki,dtype = float)],
[np.array(bki,dtype = float)]],
axis=0)
tabla = np.transpose(tabla)
np.set_printoptions(precision=digitos)
print('\n coeficientes: \n [ k, \t\t ak, \t bk ] ')
print(tabla)
# serie de Fourier expresión
serieF = a0 + 0*t
i = 1
whilenot(i>=n):
serieF = serieF + aki[i]*sym.cos(i*w0*t)
serieF = serieF + bki[i]*sym.sin(i*w0*t)
i = i+1
# serie = sym.simplify(serie)print('\n serie de Fourier fs(t): ')
# sym.pprint(serieF)print(a0)
for i inrange(1,n,1):
termino = '+ 'if aki[i]<0:
termino = '- '
aki_str = str(abs(aki[i]))
ifabs(aki[i])<casicero:
aki_str = '0'
termino = termino + aki_str+' cos('+str(i*w0)+' t)'if bki[i]>=0:
termino = termino + '+ 'else:
termino = termino + '- '
bki_str = str(abs(bki[i]))
ifabs(bki[i])<casicero:
bki_str = '0'
termino = termino + bki_str+' sin('+str(i*w0)+' t)'print(termino)
# Grafica - Para evaluacion numerica
fsn = sym.lambdify(t,serieF)
ftn = sym.lambdify(t,ft)
# Evaluacion para grafica
ti = np.linspace(a,b,muestras)
fi = ftn(ti)
fsi = fsn(ti)
# SALIDA# Grafica f(t) y Fourier en t
figura, graf_ftT0 = plt.subplots()
graf_ftT0.plot(ti,fi,label = 'f(t)')
etiqueta = 'coeficientes = '+ str(n)
graf_ftT0.plot(ti,fsi,label = etiqueta)
graf_ftT0.set_xlabel('t')
graf_ftT0.set_ylabel('f(t)')
graf_ftT0.legend()
graf_ftT0.grid()
graf_ftT0.set_title('Serie de Fourier de f(t)')
#plt.show()# espectro de frecuencia Amplitud y fase
k_i = np.arange(0,n,1,dtype=float)
ak_i = np.array(aki,dtype=float)
bk_i = np.array(bki,dtype=float)
ck_i = np.sqrt(ak_i**2 + bk_i**2)
cfs_i = np.zeros(len(ak_i))
condicion = (abs(ak_i)>=casicero)
pendiente = -bk_i[condicion]/ak_i[condicion]
cfs_i[condicion] = np.arctan(pendiente)
tabla = np.concatenate([tabla,
np.transpose([ck_i]),
np.transpose([cfs_i])
],
axis=1)
np.set_printoptions(precision=4)
print('\n coeficientes magnitud y fase: ')
print(' [k,\t\t ak,\t bk,\t\t Ck,\t fase_k ]')
print(tabla)
# grafica de espectro de frecuencia
figura, graf_Fw = plt.subplots(2,1)
graf_Fw[0].stem(k_i,ck_i,label='Ck_magnitud')
graf_Fw[0].set_ylabel('Ck')
graf_Fw[0].set_title('Espectro de frecuencia')
graf_Fw[0].legend()
graf_Fw[0].grid()
graf_Fw[1].stem(k_i,cfs_i,label='Ck_fase')
graf_Fw[1].legend()
graf_Fw[1].set_ylabel('Ck_fase')
graf_Fw[1].set_xlabel('k')
graf_Fw[1].grid()
plt.show()