Transformada Laplace

Para una señal X(t), Su transformada de Laplace esta definida como:

X(s) = \int_{-\infty}^{\infty} x(t) e^{-st} \delta t

la señal x(t) es la inversa de la transformada X(s), que se obtiene como:

x(t) = \frac{1}{2πj} \int_{c-j\infty}^{c+j\infty} x(t) e^{st} \delta t

donde c es una constante seleccionada para asegurar la convergencia de la integral.

Los pares de ecuaciones conocidos como «Pares de la transformada de Laplace» se escriben simbólicamente como:

X(s) = L[x(t)] x(t) = L^{-1}[X(s)]

y propiedades como

L[a_1 x_1(t) + a_2 x_2 (t)] = a_1 X_1(s) + a_2 X_2 (s)

que muestra la linealidad de la transformada

Región de convergencia (ROC)

También conocida como región de existencia, es el grupo de valores de s, en la region del plano complejo, para los cuales la integral de la transformada converge.

Ejemplo

referencia: Lathi ejemplo 4.1. pdf231

Para una señal x(t) = e-atu(t), encuentre la transformada X(s) y su región de convergencia (ROC)

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

pero tomando en cuenta que u(t)=0 para t<0 y u(t)=1 para t≥0:

X(s) = \int_{0}^{\infty} e^{-at} \mu(t) e^{-st} \delta t X(s) = \int_{0}^{\infty} e^{-at} e^{-st} \delta t = \int_{0}^{\infty} e^{-(s+a)t} \delta t = \frac{1}{s+a} \Big|_{0}^{\infty}

Si s es un número complejo y t →∞, el término exponencial no necesariamente desaparece.

El número complejo se puede expresar como:

z = α + jβ
e-zt = e-(α + jβ)t = e-αte-jβt.

El valor |e-jβt| = 1 sin importar el valor de βt.
Entonces si t→∞, e-zt →0 solo si α>0
y e-zt →∞ en el caso que α<0.

La región de convergencia (ROC) de X(s) es que la parte real s>-a.

Tabla de transformadas de Laplace (wikipedia)

Tabla de transformadas y propiedades de Laplace (Universidad de Córdova)

Propiedad de diferenciación

\frac{\delta x}{\delta t} \Leftrightarrow sX(s) - x(0^{-}) \frac{d^2x}{dt^2} \Leftrightarrow s^2X(s) - sx(0^{-}) - x'(0^{-}) \frac{\delta^3x}{\delta t^3} \Leftrightarrow s^3 X(s) - s^2 x(0^{-}) - sx'(0^{-}) - x''(0^{-})

Transformada Laplace Ejercicio01

Ejemplo 1:

Referencia: Oppenheim ejemplo 9.37 p718/pdf746, semejante al ejercicio sistema-modelo LTI Lathi 1.8-1 pdf/p.80. Oppenheim problema 2.61c pdf/p.191

Suponga un sistema LTI causal que se describe mediante la ecuación diferencial:

\frac{d^2y(t)}{dt} + 3 \frac{dy(t)}{dt} + 2y(t) = x(t)

junto con la condición de reposo inicial:

H(s) = \frac{1}{s^2 + 3s +2}

Suponga que x(t) = α μ(t) , entonces la entrada del sistema es x(s) = α

Siempre que las señales de entrada son cero para t<0, la convolución es la multiplicación de los términos. (no aplica para valores t<0)

y(s) = H(s) X(s) = \Big[\frac{1}{s^2 + 3s +2} \Big]\Big[\frac{\alpha}{s}\Big] y(s) = \frac{\alpha}{2}\frac{1}{s} - \alpha\frac{1}{s+1} +\frac{\alpha}{2}\frac{1}{s+2}

por lo que aplicando la anti transformada:

y(t) = \alpha \Big [ \frac{1}{2} - e^{-t} + \frac{1}{2} e^{- 2t} \Big] u(t)

Usando Python y Sympy

Para crear la expresión polinómica simple, no usamos el factor α, que se añade luego de obtener el resultado, se obtiene:

    a         a      a 
--------- - ----- + ---
2*(s + 2)   s + 1   2*s

resultado al que hay que multiplicar por la constante α todos sus términos.

las instrucciones son:

# Oppenheim ejemplo 9.37 p718/pdf746
import sympy as sym

s = sym.Symbol('s')
a = sym.Symbol('a')
Hs = (1/(s**2+3*s+2))*(a/s)
separa = Hs.apart(s)

sym.pprint(separa)

Ejemplo 2

Referencia: Oppenheim ejemplo 9.38 p719/747

Considere el sistema caracterizado por la ecuación diferencial con condiciones iniciales.

\frac{d^2y(t)}{dt^2} + 3 \frac{dy(t)}{dt} + 2y(t) = x(t) y(0^{-}) = \beta \text{, } y'(0^{-}) = \gamma

sea x(t) = αμ(t)

Aplicando la transformada de laplace en ambos lados de la ecuación con condiciones iniciales y aplicando la propiedad de diferenciación:

\frac{d^2y(t)}{dt^2} = s^2Y(s) - y(0^{-})s - y'(0^{-}) = s^2Y(s) - \beta s - \gamma \frac{dy(t)}{dt} = sY(s) - \beta

resultados que se insertan en la ecuacion original:

[s^{2}Y(s) - \beta s - \gamma] + 3[sY(s) - \beta] + 2Y(s) = \frac{\alpha}{s}

para simplificar como:

[s^{2}+3s +2] Y(s) - \beta s - (\gamma + 3 \beta) = \frac{\alpha}{s} [s^{2}+3s +2] Y(s) = \frac{\alpha}{s}+ \beta s + (\gamma + 3 \beta) (s+1)(s+2)Y(s) = \frac{\alpha}{s}+ \beta s + (\gamma + 3 \beta) Y(s) = \frac{\alpha}{s(s+1)(s+2)}+ \frac{\beta s}{(s+1)(s+2)} + \frac{(\gamma + 3 \beta)}{(s+1)(s+2)}

se encuentra que:

Y(s) = \frac{\beta (s+3)}{(s+1)(s+2)} + \frac{\gamma}{(s+2)(s+2)} + \frac{\alpha}{s(s+2)(s+2)}

donde Y(s) es la transformada unilateral de Laplace de y(t)

el último término, representa la respuesta del sistema LTI causal y la condición de reposo inicial, conocido también como la respuesta en estado cero.

Una interpretación análoga se aplica a los primeros dos términos del miembro derecho de la ecuación, que representa la respuesta del sistema cuando la entrada es cero (α=0) , conocida también como respuesta a entrada cero.

Observe que la respuesta a entrada cero es una función lineal de los valores de las condiciones iniciales. Demostrando que la respuesta es la superposición del estado cero y respuesta a entrada cero.

si α=2, β = 3 y γ=-5, al realizar la expansión en fracciones parciales encontramos que:

Y(s) = \frac{1}{s} - \frac{1}{s+1} + \frac{3}{s+2}

con la aplicación de las tablas de transformadas se obtiene:

y(t) = \Big[ 1 - e^{-t} + 3 e^{-2t} \Big] \mu (t) \text{, para } t\gt 0

usando Python y Sympy

se puede trabajar la ecuación para obtener:

ecuacion de entrada
    2                                   a
Ys*s  + 3*Ys*s + 2*Ys - b*s - 3*b - g = -
                                        s
ecuacion Ys:  
        2               
 a + b*s  + 3*b*s + g*s 
[----------------------]
      / 2          \    
    s*\s  + 3*s + 2/    
fracciones parciales de Ys: 
 a    a - 2*b - 2*g   a - 2*b - g
--- + ------------- - -----------
2*s     2*(s + 2)        s + 1   
sustituye a=2, b= 3 y g=-5
  3       1     1
----- - ----- + -
s + 2   s + 1   s
>>>

con las siguientes instrucciones, presentadas por cada sección de la parte de desarrollo analítico, para ir comparando los resultados.

# Oppenheim ejemplo 9.38 p719/747
import sympy as sym

s,a,b,g  = sym.symbols('s a b g')
Ys = sym.Symbol('Ys')

d2y = (s**2)*Ys - b*s - g
dy  = s*Ys - b

ecuacion = sym.Eq(d2y + 3*dy + 2*Ys, a/s)
ecuacionYs = sym.solve(ecuacion,Ys)
parciales = ecuacionYs[0].apart(s)

solucion = parciales.subs({a:2,b:3,g:-5})

print('ecuacion de entrada')
sym.pprint(ecuacion)
print('ecuacion Ys:  ')
sym.pprint(ecuacionYs)
print('fracciones parciales de Ys: ')
sym.pprint(parciales)
print('sustituye a=2, b= 3 y g=-5')
sym.pprint(solucion)

La definición de la ecuacion se realiza usando las variables como símbolos, incluyendo Y(s) como Ys. En éste ejercicio no existe X(s), por lo que no se la incluyó.

Luego se busca la expresión solución para Ys, por simplicidad se la cambia a fracciones parciales, asi es más sencillo usar la tabla de transformadas de Laplace.

Para una solución con los valores de a,b y g, se usa un diccionario {} con las parejas de variables y sus valores. Compare los resultados con lo expuesto en la parte teórica.

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.

FFT python scipy

La librería scipy de python incluye algoritmos para la transformada rápida de Fourier FFT.

Como ejemplo, el resultado para

x(t)=\cos(2\pi fs t)

con fs=20 se presenta en un gráfico con:
– la señal en el tiempo,
– el componente real de FFT
– el componente imaginario de FFT
– la magnitud de FFT- la fase de FFT.

Observe las escalas:la frecuencia se muestra en Hz.

Se recomienda probar el script de python con otras funciones.

# 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):
    x=np.cos(2*np.pi*fs*t)
    return(x)

# PROGRAMA
# INGRESO
t0=0
tn=0.5 #float(input('rango segundos [0,tn]:'))
n= 1024 #int(input(' n muestras en rango:'))

# PROCEDIMIENTO
dt=(tn-t0)/n    # intervalo de muestreo
# Analógica Referencia para mostrar que es par
fs=20
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]
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] imaginario
xfimag=(1/n)*np.imag(xf)
# x[w] magnitud
xfabs=(1/n)*np.abs(xf)
# x[w] angulo
xfangle=(1/n)*np.unwrap(np.angle(xf))

#SALIDA
plt.figure(1)       # define la grafica
plt.suptitle('Transformada Rápida Fourier FFT')

plt.subplot(321)    # grafica de 3x2, subgrafica 1
plt.ylabel('xanalog[t]')
plt.xlabel('tiempo')
plt.plot(t,xanalog)
plt.margins(0,0.05)
plt.grid()

ventana=0.2 # ventana de frecuencia a observar alrededor f=0
ra=int(len(frq)*(0.5-ventana))
rb=int(len(frq)*(0.5+ventana))
plt.subplot(323)    # grafica de 3x2, subgrafica 3
plt.ylabel('x[f] real')
plt.xlabel(' frecuencia (Hz)')
plt.plot(frq[ra:rb],xfreal[ra:rb])
plt.margins(0,0.05)
plt.grid()

plt.subplot(325)    # grafica de 3x2, subgrafica 5
plt.ylabel('x[f] imag')
plt.xlabel(' frecuencia (Hz)')
plt.plot(frq[ra:rb],xfimag[ra:rb])
plt.margins(0,0.05)
plt.grid()

plt.subplot(324)    # grafica de 3x2, subgrafica 4
plt.ylabel('x[f] magnitud')
plt.xlabel(' frecuencia (Hz)')
plt.plot(frq[ra:rb],xfabs[ra:rb])
plt.margins(0,0.05)
plt.grid()

plt.subplot(326)    # grafica de 3x2, subgrafica 6
plt.ylabel('x[f] fase')
plt.xlabel(' frecuencia (Hz)')
plt.plot(frq[ra:rb],xfangle[ra:rb])
plt.margins(0,0.05)
plt.grid()

plt.show()

Para el cálculo de la transformada se han usado los valores de t≥0

la transformada calcula los valores positivos y luego los negativos, por lo que se usa fftshift(xf) antes de graficar para poner el orden apropiado.

Las k muestras debe ser un valor de 2n tal como 512, 1024, 2048, …

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

 

Fracciones Parciales – Separar con Python

Referencia: Oppenheim 9.4 p671/pdf699. Lathi B5 pdf21, Schaum/HSu 3.8.C p120

Se pueden separar en fraciones parciales del siguiente ejemplo:

H(s) = \frac{2s^2 +5}{s^2 +3s+2}

se requiere obtener los residuos, polos y términos directos al separar en fracciones parciales entre dos polinomios P(s)/Q(s).

\frac{P(s)}{Q(s)} = \frac{ceros(1)}{s - polos(1)}+\frac{ceros(2)}{s - polos(2)}+\text{ ... } \text{ ... }+\frac{ceros(n)}{s - polos(n)}+ganancia(s)

Se puede realizar de dos formas usando Python, en ambas se requiere definir los vectores P y Q.

Sympy – fórmulas simbólicas

La primera usando Sympy para crear las fórmulas simbólicas y luego formar una  fracción, separar las partes Sympy.apart()

Para los datos proporcionados en el ejemplo se tiene como resultado:

P = [2, 0, 5]
Q = [1, 3, 2]

Fracciones parciales:
2 - 13/(s + 2) + 7/(s + 1)

con las instrucciones sympy:

# Expansión en fracciones parciales
# P es numerador, Q es denominador

P = [2, 0, 5]
Q = [1, 3, 2]

import sympy as sym
s = sym.Symbol('s')
n = len(P)
m = len(Q)
Ps = 0
for i in range(0,n,1):
    Ps = Ps + P[i]*s**(n-i)
Qs = 0
for i in range(0,m,1):
    Qs = Qs + Q[i]*s**(m-i)
fraccion = Ps/Qs
parciales = sym.apart(fraccion)

# Salida - polinomio
print('Fracciones parciales:')
print(parciales)

Las fórmulas simbólicas permiten agrupar y realizar operaciones entre varios sistemas además de expandir, separar, y evaluar las expresiones. Util principalmente para la sección de transformadas de Laplace.

Un resumen de sympy se muestra en:

http://blog.espol.edu.ec/matg1013/formulas-simbolicas-sympy/


Scipy – Procesar señales

La segunda forma es usando libreria scipy.signal.residue con los valores de P y Q:

ceros(n)   :  [-13.   7.]
polos(n)   :  [-2. -1.]
ganancia(s):  [ 2.]
>>>

Con los resultados el equivalente en Laplace es:

transformado al dominio del tiempo,
H(t) = (−13e−2t + 7e−t)u(t) + 2δ(t)

las instrucciones adicionales a la sección anterior son:

# Usando scipy.signal.residue

import scipy.signal as senal
ceros,polos,ganancia = senal.residue(P,Q)

# SALIDA - coeficientes
print()
print('ceros(n)   : ', ceros)
print('polos(n)   : ', polos)
print('ganancia(s): ', ganancia)

Diagramas de bloques con Z

Considere el sistema LTI causal descrito mediante:

y[n] - (1/4) y[n-1] = x[n]

Ejemplo Openheim 10.28. con función del sistema:

H(z)  = 1/ [1- (1/4) z-1]

Aqui z-1 esla función del sistema de un retardo unitario. el diagrama de bloques en la figura contiene un lazo de retroalimentación.

y[z] - (1/4) z-1 y[z] = x[z]

y[z]  = x[z] + (1/4) z-1 y[z]

Ejemplo Openheim 10.29. Considere un sistema LTI causal con función del sistema:

H(z) = [ 1 - 2z-1 ] / [ 1 - (1/4)z-1 ] 
     = {1/ [ 1 - (1/4)z-1 ]} [ 1 - 2z-1 ]

Diagramas de Bloques con s

Representaciones en diagramas de bloques para los sistemas LTI causales descritos por ecuaciones diferenciales y funciones racionales del sistema. (Oppenheim 9.8.2)


Ejemplo 1

Referencia: Oppenheim ejemplo 9.28 p708/pdf736

Considere un sistema LTI causal

\frac{\delta y(t)}{\delta t} + 3y(t) = x(t) Dy(t) + 3y(t) = x(t) Dy(t) = x(t) -3y(t) y(t) = \frac{1}{D}[x(t) - 3y(t)]

usando s en lugar de D

y(s) = \frac{1}{s} [x(s) - 3 y(s)]

usando la transformada de Laplace, cambiando s por d/dt desde la primera ecuación, la función del sistema es:

( s+3) y(s) = x(s)

por lo que y(s)/x(s) es:

H(s) = \frac{1}{s+3}

al cambiar el signo en el coeficiente de la multiplicación de la trayectoria de retroalimentación, obtenemos una forma de diagrama de bloques similar a la anterior, pero con la que se verifica que:

H(s) = \frac{\frac{1}{s}}{1+\frac{3}{s}} = \frac{1}{s+3}

que confirma que:

\frac{Y(s)}{X(s)} = H(s) = \frac{H_1(s)}{1+H_1(s)H_2(s)}

siendo H1 el bloque superior, y H2 el bloque inferior en el diagrama.


Ejemplo 2

Referencia: Oppenheim ejemplo 9.30 p711/pdf739

Considere el sistema causal de segundo orden con la funcion del sistema:

H(s) = \frac{1}{(s+1)(s+2)} = \frac{1}{s^2 + 3s +2}

la entrada x(t) y la salida y(t) para este sistema satisfacen la ecuación diferencial:

\frac{\delta^2(t)}{\delta t^2} + 3 \frac{\delta y(t)}{dt} + 2 y(t) = x(t) (s^2 + 3s + 2)y(t) = x(t)

pero también se puede armar como:

H(s) = \frac{1}{(s+1)(s+2)} = \frac{k_1}{s+1} + \frac{k_2}{s+2}

siendo el numerador :

k_1s + 2k_1 + k_2s + k_2 = (k_1 + k_2)s + (2k_1 +k_2)

que se observa que debe dar:
(k_1 + k_2)s + (2k_1 +k_2) = 0s + 1

entonces

k_1 + k_2 = 0 2k_1 + k_2 = 1

con lo que

k1 = 1 y
K2 = -1
H(s) = \frac{1}{s+1} -\frac{1}{s+2}

Si

H(s) = \frac{1}{(s+1)(s+2)}

se escribe como

H(s) = \Big[\frac{1}{s+1}\Big]\Big[\frac{1}{s+2}\Big]

que se representa como la serie o cascada de dos sistemas de primer orden.

que por cierto, también es la convolución de dos sistemas en serie.

También es posible realizar las fracciones parciales con Sympy y se obtiene:

    1       1  
- ----- + -----
  s + 2   s + 1
>>>

usando las instrucciones:

# Oppenheim 9.30 p711/pdf739
import sympy as sym

s = sym.Symbol('s')
Hs = 1/((s+1)*(s+2))
separa = Hs.apart()

sym.pprint(separa)

Ejemplo 3

referencia: Oppenheim ejemplo 9.31.p712/pdf739

Considere la función del sistema

H(s) = \frac{2s^2 + 4s - 6}{s^2 + 3s + 2}

se puede reescribir como:

H(s) = \Big[\frac{1}{s^2 + 3s + 2} \Big] \Big[ 2s^2 + 4S - 6 \Big]

y si las derivadas se toman del sistema anterior, se tiene una forma mas simplificada:

también, con un poco de trabajo, reordenando las ecuaciones, se podría convertir en:

H(s) = \Big[ 2\frac{s-1}{s+2} \Big] \Big[\frac{s+3}{s+1} \Big]

o también en :

H(s) = 2 + \frac{6}{s+2} - \frac{8}{s+1}

que requiere un sistema en forma paralelo.
QUEDA COMO TAREA HACER EL DIAGRAMA.

Se presenta el ejercicio usando Sympy:

      6       8  
2 + ----- - -----
    s + 2   s + 1
>>>

que se obtiene usando las instrucciones:

# Oppenheim 9.31.p712/pdf739
import sympy as sym

s = sym.Symbol('s')
Hs = (2*s**2+4*s-6)/(s**2+3*s+2)
separa = Hs.apart()

sym.pprint(separa)

Modulación Delta Sigma

La modulación Delta-Sigma (∑Δ ) codifica una señal analógica a digital generando una secuencia de +1 y -1 (impulsos) que representan la diferencia entre la señal analógica muestreada y la señal digital acumulada. Tiene aplicaciones en sintetizadores de frecuencia, fuentes de poder conmutadas y controladores de motor.

El sistema en diagrama de bloques en simulink se muestra en la figura:

La señal de entrada (Sine Wave) se cuantifica en (Zero-Order Hold).
Luego se obtiene el signo (sign) de la diferencia entre la señal de entrada y la señal digital acumulada (∑) que genera la secuencia de +1 y -1 en la salida del codificador.
El valor de la ganancia Δ=0.3 representa el valor del escalón de subida o bajada de la señal digital acumulada, representada en color magenta.

La secuencia de salida +1 y -1, de nuevo acumulada en el decodificador, permite obtener una versión digital cercana a la señal analógica inicial.

Simulación usando Python: CODIFICADOR

La señal de entrada para el ejemplo tiene frecuencia fs=1 Hz, Δ=deltaY=0.3, rango [0,tn] hasta 2 segundos con divisiones k=40 en el rango de tiempo.

>>> 
rango [0,tn]:2
Secciones en el rango k:40
[ 0  1  1  1  1 -1  1 -1 -1 -1 -1 -1 -1 -1 -1  1 -1  1  1  1  1  1  1  1  1 -1  1 -1 -1 -1 -1 -1 -1 -1 -1  1 -1  1  1  1]
[  0.05   0.3   40.  ]

El script genera la gráfica y los archivos para datos y parámetros en formato texto.

La señal de entrada y los valores pueden ser modificados en el script adjunto:

# Modulacion Delta - Codificador
# entrada x(t), salida: y[n]
# propuesta:edelros@espol.edu.ec
import numpy as np
import matplotlib.pyplot as plt

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

# PROGRAMA para la modulación
# INGRESO
t0=0
tn=float(input('rango [0,tn]:'))
k=int(input('Secciones en el rango k:'))

# PROCEDIMIENTO
# Analógica Referencia
fs=1
m=500  # puntos en eje t para gráfico analógica
dm=(tn-t0)/m
t=np.arange(0,tn,dm)  # eje tiempo analógica
xanalog=np.zeros(m, dtype=float)
for i in range(0,m):
    xanalog[i]=entradax(t[i],fs)

# Codificar Sigma-Delta
deltaY=0.3          # Tamaño delta en eje Y
deltaT=(tn-t0)/k    # Tamaño delta en eje tiempo
td=np.arange(0,tn,deltaT)    # tiempo muestreo
muestra=np.zeros(k, dtype=float)    # analógica para comparar
xdigital=np.zeros(k, dtype=float)   # digital
ysalida=np.zeros(k, dtype=int)      # Salida de +1|-1

td[0]=t0
muestra[0]=entradax(td[0],fs)
ysalida[0]=0
for i in range(1,k):
    muestra[i]=entradax(td[i],fs) # referencia analógica
    diferencia=muestra[i]-xdigital[i-1]
    if (diferencia>0):
        bit=1
    else:
        bit=-1
    xdigital[i]=xdigital[i-1]+bit*deltaY
    ysalida[i]=bit
parametros=np.array([deltaT,deltaY,k])

# SALIDA
print(ysalida)
print(parametros)
np.savetxt('sigmadelta_datos.txt',ysalida,fmt='%i')
np.savetxt('sigmadelta_parametros.txt',parametros)
# Graficar
plt.figure(1)       # define la grafica
plt.suptitle('Codificador Sigma-Delta')
plt.subplot(211)    # grafica de 2x1 y subgrafica 1
plt.ylabel('x(t), x[n]')
xmax=np.max(xanalog)+0.1*np.max(xanalog) # rango en el eje y
xmin=np.min(xanalog)-0.1*np.max(xanalog)
plt.axis((t0,tn,xmin,xmax)) # Ajusta ejes
plt.plot(t,xanalog, 'g')
plt.axis((t0,tn,xmin,xmax))
plt.plot(td,xdigital,'bo')  # Puntos x[n]
plt.step(td,xdigital, where='post',color='m')
plt.subplot(212)    # grafica de 2x1 y subgrafica 2
plt.ylabel('y[n]')
plt.axis((0,k,-1.1,1.1))
plt.plot(ysalida, 'bo')     # Puntos y[n]
puntos=np.arange(0,k,1)     #posicion eje x para escalon
plt.step(puntos,ysalida, where='post')
plt.show()

DECODIFICADOR

Para observar el resultado, se decodifican los datos y se obtiene el siguiente resultado:

>>> 
entrada: [ 0  1  1  1  1 -1  1 -1 -1 -1 -1 -1 -1 -1 -1  1 -1  1  1  1  1  1  1  1  1 -1  1 -1 -1 -1 -1 -1 -1 -1 -1  1 -1  1  1  1]
salida: 
[  0.00000000e+00   3.00000000e-01   6.00000000e-01   9.00000000e-01
   1.20000000e+00   9.00000000e-01   1.20000000e+00   9.00000000e-01
   6.00000000e-01   3.00000000e-01  -1.11022302e-16  -3.00000000e-01
  -6.00000000e-01  -9.00000000e-01  -1.20000000e+00  -9.00000000e-01
  -1.20000000e+00  -9.00000000e-01  -6.00000000e-01  -3.00000000e-01
  -1.11022302e-16   3.00000000e-01   6.00000000e-01   9.00000000e-01
   1.20000000e+00   9.00000000e-01   1.20000000e+00   9.00000000e-01
   6.00000000e-01   3.00000000e-01  -1.11022302e-16  -3.00000000e-01
  -6.00000000e-01  -9.00000000e-01  -1.20000000e+00  -9.00000000e-01
  -1.20000000e+00  -9.00000000e-01  -6.00000000e-01  -3.00000000e-01]

realizadousando el siguiente script de python:

# Modulacion Sigma-Delta Decodificador
#  entrada y[n], salida: x[n]
# propuesta:edelros@espol.edu.ec
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
# Lectura de datos desde archivos
yentrada=np.loadtxt('sigmadelta_datos.txt',dtype=int)
datos=np.loadtxt('sigmadelta_parametros.txt',dtype=float)

# PROCEDIMIENTO
deltaT=datos[0] # Tamaño delta en eje tiempo
deltaY=datos[1] # Tamaño delta en eje Y
k=len(yentrada) # número de muestras
xdigital=np.zeros(k, dtype=float)
punto=np.zeros(k, dtype=int)      # número de muestra
td=td=np.zeros(k, dtype=float)    # tiempo muestreo
# DECOdifica Sigma-Delta
xdigital[0]=yentrada[0]
punto[0]=0
td[0]=0
for i in range(1,k):
    punto[i]=i
    td[i]=deltaT*i
    xdigital[i]=xdigital[i-1]+yentrada[i]*deltaY

# SALIDA
print('entrada:', yentrada)
print('salida:',xdigital,)
# Graficar
plt.figure(1)       # define la grafica
plt.suptitle('Decodifica Sigma_Delta')
plt.subplot(211)    # grafica de 2x1 y subgrafica 1
plt.ylabel('entrada: y[n]')
plt.axis((0,k-1,-1.1,1.1)) # Ajusta ejes
plt.plot(punto,yentrada,'bo') # Puntos y[n]
plt.step(punto,yentrada, where='post')
plt.subplot(212)    # grafica de 2x1 y subgrafica 2
plt.ylabel('salida: x[t]')
xmax=np.max(xdigital)+0.1*np.max(xdigital)# rango en el eje y
xmin=np.min(xdigital)-0.1*np.max(xdigital)
plt.axis((0,td[k-1],xmin,xmax)) # Ajusta ejes
plt.plot(td,xdigital,'bo')
plt.step(td,xdigital, where='post', color='m')
plt.show()

Tarea: Realice observaciones cambiando:
a) la frecuencia de la señal de entrada,
b) el valor para deltaY
c) el rango tnd) el número de secciones en el tiempo
e) la forma de la señal de entrada a triangular, diente de sierra, exponencial periódica, etc.

Delta Sigma con audio

Usando los conceptos descritos para abrir un archivo de audio .wav, realice en python lo siguiente:

a) Codificador: Realice la codificación en ∑Δ para el audio del ejemplo, y crear el archivo correspondiente al resultado.

b) Decodificador: Reconstruya la señal de audio a partir del archivo del literal anterior, de ∑Δ a .wav, y pruebe su resultado ejecutando el archivo con un programa como windows media player.

c) Realice las observaciones correspondientes al proceso realizado.

Trabajo extra, opcional

d) Cree un archivo de audio en formato .wav dictando una frase como la mostrada, y repita el proceso anterior.

e) Realice las observaciones correspondientes al proceso realizado.

Anexo: Codificador
Se adjunta un algoritmo en python de referencia al que hay que añadir la parte de manejo del audio.

– En el bloque de ingreso complete las instrucciones para leer el archivo. Observe los nombres de las variables.

# Modulacion Delta - Codificador de audio
# entrada x(t), salida: ysalida[n]
# propuesta: edelros@espol.edu.ec
import numpy as np
import matplotlib.pyplot as plt
import scipy.io.wavfile as waves
import scipy.signal as senal

# INGRESO
# archivo=......
# fsonido, sonido = ......

# PROCEDIMIENTO
# Analógica de referencia
deltaT=1/(2*np.pi*fsonido)
xanalog=sonido[:,0] # Canal izquierdo
xmax=np.max(xanalog)
nsonido=len(xanalog)
t0=0
tn=nsonido*deltaT
t=np.arange(0,tn,deltaT)

# Codificacion Sigma-Delta
salto=1 #usa todas las muestras sin saltar
deltaY=xmax//10
nd=nsonido//salto
td=np.arange(0,nd,1) # eje tiempo[n] digital
xdigital=np.zeros(nd,dtype=float)
ysalida=np.zeros(nd,dtype=int)
muestra=xanalog[0]
i=1
while not(i>=nd):
    muestra=xanalog[i*salto] # referencia analógica
    diferencia=muestra-xdigital[i-1]
    if (diferencia>0):
        bit=1
    else:
        bit=-1
    xdigital[i]=xdigital[i-1]+bit*deltaY
    ysalida[i]=bit
    i=i+1
fsalida=fsonido/salto
parametros=np.array([deltaT,deltaY,fsalida])

# SALIDA
print(ysalida)
print(parametros)
# Graba archivo
np.savetxt('sigmadelta_datos.txt',ysalida,fmt='%i')
np.savetxt('sigmadelta_parametros.txt',parametros)
# Grafico
#Escala y para el grafico
a=int(nsonido*0.35)       # RANGO ejex del gráfico
b=int(nsonido*0.355)
plt.figure(1)       # define la grafica
plt.suptitle('Codificador Sigma-Delta Audio.wav')
plt.subplot(311)    # grafica de 3x1 y subgrafica 1
plt.ylabel('x(t)')
xmax=np.max(xanalog)+0.1*np.max(xanalog)# rango en el eje y
xmin=np.min(xanalog)-0.1*np.max(xanalog)
plt.axis((a*deltaT, b*deltaT, xmin, xmax))
plt.plot(t[a:b],xanalog[a:b], 'g')
plt.ylabel('xanalogica[n]')
plt.subplot(312) # grafica de 3x1 y subgrafica 2
plt.ylabel('xdigital[n]')
mxa=np.max(xdigital[a:b])  # rango en el eje
mna=np.min(xdigital[a:b])
plt.axis((a,b,mxa,mna))
#plt.plot(td,xdigital,'bo')
plt.step(td[a:b],xdigital[a:b], where='post',color='m')
plt.subplot(313) # grafica de 3x1 y subgrafica 3
plt.ylabel('ysalida[n]')
plt.axis((a,b,-1.1,1.1))
plt.step(td[a:b],ysalida[a:b], where='post',color='b')
plt.show()

Anexo Decodificador

El resultado gráfico al decodificar el archivo de ∑Δ se muestra, se debe añadir la instrucción para generar el archivo de audio y poder escucharlo.

modifique la seccion de #Salida para crear el arhivo de audio con nombre: ‘sigmadeltaaudio.wav’.

# Modulacion Sigma-Delta Decodificador
# entrada yentrada[n], salida: x[t]
# propuesta:edelros@espol.edu.ec
import numpy as np
import matplotlib.pyplot as plt
import scipy.io.wavfile as waves

# INGRESO
yentrada=np.loadtxt('sigmadelta_datos.txt',dtype=int)
datos=np.loadtxt('sigmadelta_parametros.txt',dtype=float)
deltaD=datos[0]
deltaY=datos[1]
fsonido=datos[2]

# PROCEDIMIENTO
k=len(yentrada) # número de muestras para gráfica
# Arreglos para datos
xdigital=np.zeros(k, dtype='int16')
punto=np.zeros(k, dtype=int)
td=td=np.zeros(k, dtype=float) # eje n digital
# calcula los valores para digital
xdigital[0]=yentrada[0]
punto[0]=0
td[0]=0
for i in range(1,k):
 punto[i]=i
 td[i]=deltaD*i
 xdigital[i]=xdigital[i-1]+yentrada[i]*deltaY

# SALIDA
print(xdigital)
#Archivo de audio
# CREAR EL ARCHIVO DE AUDIO ***********************

#Gráfico
a=int(k*0.25) # RANGO ejex del gráfico
b=int(k*0.252)
plt.figure(1) # define la grafica
plt.suptitle('Decodifica Sigma_Delta')
plt.subplot(211) # grafica de 2x1 y subgrafica 1
plt.ylabel('entrada: y[n]')
plt.axis((a,b,-1.1,1.1))
#plt.plot(punto[a:b],yentrada[a:b],'b')
plt.step(punto[a:b],yentrada[a:b], where='post', color='b')
plt.subplot(212) # grafica de 2x1 y subgrafica 2
plt.ylabel('salida: x[t]')
xmax=np.max(xdigital[a:b]) # rango en el eje
xmin=np.min(xdigital[a:b])
plt.axis((a*deltaD,b*deltaD,xmin,xmax))
#plt.plot(td[a:b],xdigital[a:b],'m')
plt.step(td[a:b],xdigital[a:b], where='post', color='m')
plt.show()

La gráfica muestra solo el intervalo de tiempo entre a=0.25 y b=0.252 del total de tiempo del audio, el propósito es la observacion visual complementará a la de escuchar el audio resultante.