Series de Fourier

Referencia: Lathi 6.1 pdf406, Chapra 5ed Ejemplo 19.2 p548/pdf572

La serie de Fourier aproxima una señal o función contínua mediante una serie infinita de sinusoides.

f(t) = a_0 + \sum_{k=1}^{\infty}[a_k \cos( k \omega_0 t) + b_k \sin(k \omega_0t)]

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

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.

f(t) = \begin{cases} -1 && -T/2<t<-T/4 \\ 1 && -T/4<t<T/4 \\ -1 && T/4<t<T/2 \end{cases}


Coeficientes ak y bk

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:

Fórmulas simbólicas – Sympy

El desarrollo a papel y lápiz se deja como tarea.

Usando Python se obtiene los siguientes resultados:

 expresión ak:
/1.27323*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/pdf572
import 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 evaluan los n coeficientes ak y bk.

 coeficientes ak: 
[0, 1.27323954473516, 1.55926873300775e-16, 
  -0.424413181578388]

 coeficientes bk: 
[0, 0, 0, 0]

Instrucciones Python

las instrucciones son adicionales al bloque anterior. La evaluación se mantiene usando las expresiones simbólicas usando la instruccion .subs()

# evalua los coeficientes
a0 = ak.subs(k,0)
b0 = bk.subs(k,0)
aki = [a0]
bki = [b0]
i = 1
while not(i>=n):
    avalor = ak.subs(k,i)
    bvalor = bk.subs(k,i)
    aki.append(avalor)
    bki.append(bvalor)
    i = i+1
print('\n coeficientes ak: ')
print(aki)
print('\n coeficientes bk: ')
print(bki)

Expresión de la Serie de Fourier

Encontrados los coeficientes ak y bk, se los usa en la expresión principal

f(t) = a_0 + \sum_{k=1}^{\infty}[a_k \cos( k \omega_0 t) + b_k \sin(k \omega_0t)]

obteniendo la siguiente expresión para la serie de Fourier

 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

# serie de Fourier
serieF = a0 + 0*t 
i = 1
while not(i>=n):
    serieF= serie + aki[i]*sym.cos(i*w0*t)
    serie = serie + bki[i]*sym.sin(i*w0*t)
    i = i+1
# serie = sym.simplify(serie)

print('\n serie de Fourier fs(t): ')
sym.pprint(serie)

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 instruccion 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 la instrucción sym.lambdify()

Intrucciones en Python

# 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()

LTI en dominio frecuencia

Referencia: Oppenheim Ejemplo 4.15 (p.345 pdf)

Ejemplo: Considere un sistema continuo LTI con respuesta al impulso:

h(t) = δ(t-t0)

La respuesta en frecuencia es la transformada de Fourier :

H(jω) = e-jωt0

A cualquier entrada x(t) con transformada X(jω) , la salida en dominio de la frecuencia es:

Y(jω) = H(jω)X(jω) = e-jωt0 X(jω)

que usando la propiedad de desplazamiento:

y(t)=x(t-t0)

Referencia: Oppenheim Ejemplo 4.18 (p.346 pdf)

El Filtrado selectivo en frecuencia se puede representar en un filtro PASA-BAJO ideal:

H(jω) = { 1    |ω|< ωc
        { 0    |ω|> ωc

cuya respuesta en el dominio de la frecuencia es:

h(t) = sen(ωc t)/πt

Como ejemplo: para una frecuencia fs=20 y usando la transformada para una función par se obtiene:

Observaciones: h(t) no es cero para t<0, en consecuencia el filtro pasabajos ideal no es causal. NO es fácil aproximarse mucho a un filtro ideal con componentes reales, por lo que se contruye más facilmente un filtro no ideal.

FFT de una señal coseno

Señal de entrada:

x(t) = \cos (\omega _0 t)

ω0 = (2π fs) en radianes, para el ejemplo fs=20 Hz

x(ω) =π[δ(ω-ω0)+δ(ω+ω0)] con eje en radianes.
La gráfica se muestra en Hz. fs0/2π
Se muestran las partes:
– real e imaginaria de la transformada y
– magnitud y fase

Script de python para usar la transformada rapida de Fourier.
En la seccion de función de entrada, se muestra para cuatro funciones: coseno, constante, impulso y sinc. Para usar una de ellas se elimina el simbolo de comentario (##)

# Transformadas Rapida de Fourier
# entrada funcion par
# Propuesta: edelros@espol.edu.ec
import numpy as np
import scipy.fftpack as fourier
import matplotlib.pyplot as plt

# Definir la funcion de Entrada par
def entradax(t,fs):
    # función matemática CAMBIAR AQUI
    x=np.cos(2*np.pi*fs*t)

##    x=1  #constante
    
##    x=0   # impulso
##    if abs(t)<(1/fs):
##        x=1

##    if t==0:  #para sinc aumentar k=8
##        t=(1/fs)/100000
##    x=(np.sin(2*np.pi*fs*t)/(np.pi*t))

    # función matemática CAMBIAR AQUI
    return(x)

# PROGRAMA
# INGRESO
fs=20       #float(input('frecuencia Señal Hz:'))
ancho=1     #ancho del intervalo
t0=0        #intervalo de t
tn=ancho*1/fs
resolucion=1    #resolucion de transformada

# PROCEDIMIENTO
n= 1024*resolucion 
dt=(tn-t0)/n   # intervalo de muestreo
# Analógica Referencia para mostrar que es par
t=np.arange(-tn,tn,dt)  # eje tiempo analógica
m=len(t)
xanalog=np.zeros(m, dtype=float)
for i in range(0,m):
    xanalog[i]=entradax(t[i],fs)

# FFT: Transformada Rapida de Fourier
# Analiza la parte t>=0 de xnalog muestras[n:2*n]
# Transformada de Fourier
xf=fourier.fft(xanalog[n:2*n])
xf=fourier.fftshift(xf)
# Rango de frecuencia para eje
frq=fourier.fftfreq(n, dt)
frq=fourier.fftshift(frq)

# x[w] real
xfreal=(1/n)*np.real(xf)
# x[w] angulo
xfimag=(1/n)*np.imag(xf)
# x[w] magnitud
xfabs=(1/n)*np.abs(xf)
# x[w] angulo en radianes
xfangle=np.angle(xf)# np.unwrap(np.angle(xf))


#SALIDA
plt.figure(1)   # define la grafica
plt.suptitle('Señal de Entrada')
plt.ylabel('xanalog[t]')
plt.xlabel('tiempo')
plt.plot(t,xanalog)
plt.margins(0,0.05)
plt.grid()

plt.figure(2)   # define la grafica
plt.suptitle('Transformada Rápida Fourier FFT')
ventana=ancho*0.01 # ventana de frecuencia alrededor f=0
ra=int((n/2)*(1-ventana)) #rango de indices
rb=int((n/2)*(1+ventana))+1
a=frq[ra] #rango de ejes uniformes entre graficos
b=frq[rb]
c=np.min([-0.05,1.1*np.min(xfreal)])
d=1.1*np.max(xfabs)
ejes=[a,b,c,d]

plt.subplot(221)    # grafica de 2x2, subgrafica 1
plt.ylabel('x[f] real')
plt.xlabel(' frecuencia (Hz)')
plt.axis(ejes)
plt.stem(frq[ra:rb],xfreal[ra:rb])
plt.grid()

plt.subplot(223)    # grafica de 2x2, subgrafica 2
plt.ylabel('x[f] imag')
plt.xlabel(' frecuencia (Hz)')
c=np.min([-0.05,1.1*np.min(xfimag)])
ejes=[a,b,c,d]
plt.axis(ejes)
plt.stem(frq[ra:rb],xfimag[ra:rb])
plt.grid()

plt.subplot(222)    # grafica de 2x2, subgrafica 3
plt.ylabel('x[f] magnitud')
plt.xlabel(' frecuencia (Hz)')
c=np.min([-0.05,1.1*np.min(xfabs)])
ejes=[a,b,c,d]
plt.axis(ejes)
plt.stem(frq[ra:rb],xfabs[ra:rb])

plt.grid()

plt.subplot(224)    # grafica de 2x2, subgrafica 4
plt.ylabel('x[f] fase radianes')
plt.xlabel(' frecuencia (Hz)')
c=1.1*np.min(xfangle)
d=1.1*np.max(xfangle)
ejes=[a,b,c,d]
plt.axis(ejes)
plt.stem(frq[ra:rb],xfangle[ra:rb])
plt.grid()

plt.show()

En el caso del seno(ω0 t), las graficas resultantes son:


con transformada:

continuar el ejercicio con las otras funciones.

Periodicidad señales analógicas y digitales

Considere la exponencial compleja discreta con frecuencia (ω0 + 2π)

e^{j(\omega _0 + 2 \pi)n} =e^{j 2 \pi n} e^{j \omega _0 n} = e^{j \omega _0 n}

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 e0t 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.ec
import numpy as np
import math
import matplotlib.pyplot as plt

# Definir la funcion para el ejemplo
def ejemplo_funcion(f,t):
    # función matemática CAMBIAR AQUI
    y=np.cos(f*t)
    # función matemática CAMBIAR AQUI
    return(y)

# Programa para graficar la función
# INGRESO
rango=float(input('rangos en periodos de analogica:'))
fa=float(input('frecuencia de analógica:'))
fd=float(input('frecuencia digital:'))

# PROCEDIMIENTO
n=500  # puntos en eje t para el gráfico analogica
t0=-rango*2*np.pi*fa/2
tn=rango*2*np.pi*fa/2
delta=(tn-t0)/n
deltaD=(2*np.pi*fa)/(2*fd)
nd=int((tn-t0)//deltaD +1)   # número de muestras para gráfica
# Dimensiona los arreglos
t=np.zeros(n, dtype=float)  # para la analogica
yanalog=np.zeros(n, dtype=float)
td=np.zeros(nd, dtype=float)    # Para la digital
ydigital=np.zeros(nd, dtype=float)
# valores para analogica
for i in range(0,n):
    t[i]=t0+delta*i     # para el eje x
    yanalog[i]=ejemplo_funcion(fa,t[i])
#valores para digital
for i in range(0,nd):
    td[i]=t0+deltaD*i   #para el eje x
    ydigital[i]=ejemplo_funcion(fa,td[i])

# SALIDA
#Escala y para el grafico
ymax=np.max(yanalog)+0.1*np.max(yanalog)# rango en el eje y
ymin=np.min(yanalog)-0.1*np.max(yanalog)
plt.figure(1)       # define la grafica
plt.suptitle('Señal Analogica vs digital')
plt.subplot(211)    # grafica de 2x1 y subgrafica 1
plt.ylabel('y(t)')
plt.axis((t0,tn,ymin,ymax))
plt.plot(t,yanalog)
plt.subplot(212)    # grafica de 2x1 y subgrafica 2
plt.ylabel('Digital: y[t]')
plt.axis((t0,tn,ymin,ymax))
plt.plot(td,ydigital,'ro')
plt.show()