1.5 Espectro – Producto de sinusoides, AM y BW con Sympy-Python

[ producto cos ] [ ejercicio ] [ analítico ] [ algoritmo ] [gráfica ] [ AM y BW]
..


1. Producto de Señales

Un modelo de señal de gran utilidad es el producto de dos sinusoides, usada en difusión radial conocida como Modulación de Amplitud (AM).

[ producto cos ] [ ejercicio ] [ analítico ] [ algoritmo ] [gráfica ] [ AM y BW]
..


2. Ejemplo – Espectro de Producto de sinusoides

Referencia: McClellan ejemplo 3.2 p75

Realice el espectro de frecuencias para una señal resultado del producto de dos sinusoides con frecuencias 1/2 Hz y 5 Hz.

x(t) = \cos (\pi t) \sin(10 \pi t)

[ producto cos ] [ ejercicio ] [ analítico ] [ algoritmo ] [gráfica ] [ AM y BW]
..


3. Desarrollo Analítico

Lo primero es convertir x(t) usando la fórmula de Euler:

x(t) = \Big( \frac{e^{j \pi t} +e^{-j \pi t}}{2} \Big)\Big( \frac{e^{j 10 \pi t} - e^{-j 10 \pi t}}{2j} \Big) = \frac{1}{4}e^{-j \pi /2}e^{j 11 \pi t} + \frac{1}{4}e^{-j \pi /2}e^{j 9 \pi t} + \frac{1}{4}e^{j \pi /2}e^{-j 9 \pi t} + \frac{1}{4}e^{j \pi /2}e^{-j 11 \pi t} = \frac{1}{2}\cos (11 \pi t - pi/2) +\frac{1}{2}\cos(9 \pi t - \pi/2)

Se observa que para el espectro de frecuencias se tienen 4 componentes, que las frecuencias originales de 5 Hz y 1/2 Hz ya no se muestran en el espectro

espectro Senales Producto Cos 01

La gráfica de domino en tiempo con los componentes individuales es:

espectro Senales Producto Cos t 02

[ producto cos ] [ ejercicio ] [ analítico ] [ algoritmo ] [gráfica ] [ AM y BW]
..


4. Algoritmo en Python

El algoritmo para el ejercicio reutiliza las funciones descritas en la sección anterior para el espectro de frecuencias. Por lo que éste caso  es una aplicación del concepto.

El resultado del algoritmo es:

x(t): cos(pi*t)*cos(10*pi*t)
x_senales: 
senal:   cos(pi*t)*cos(10*pi*t)
  euler: exp(11*I*pi*t)/4 + exp(9*I*pi*t)/4 + exp(-9*I*pi*t)/4 + exp(-11*I*pi*t)/4
x_simplif: 
  cos(9*pi*t)/2 + cos(11*pi*t)/2
x_espectro:
freq : [-5.5 -4.5  4.5  5.5]
ampl : [1/4 1/4 1/4 1/4]
fase : [0 0 0 0]
etiq : ['1/4' '1/4' '1/4' '1/4']

Instrucciones e Python

# ejemplo 3.2 p75 Espectro de dos lados
# telg1034 DSP fiec-espol edelros@espol.edu.ec
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
import telg1034 as dsp

# variables
from telg1034 import t,A,w,f,p,pi,DosPi,I,equivalentes

# INGRESO
x0 = 0
x1 = sym.cos(pi*t)
x2 = sym.cos(10*pi*t)
x  = x1*x2

# PROCEDIMIENTO
x_senales = dsp.cos_x_components_add(x)
x_conteo = len(x_senales)
Xe_senales = dsp.cos_to_euler_one_term(x_senales)
Xe_espectro = dsp.cos_spectrum_list(x_senales)

# x_senales simplificados
Xe_conteo = len(Xe_senales)
x_senales_simple = []
x_simplif = sym.S.Zero
for i in range(Xe_conteo):
    unaSenal = Xe_senales[i]
    Xe_to_cos = unaSenal.rewrite(sym.cos)
    x_simplif = x_simplif + Xe_to_cos

# SALIDA
print('x(t):',x)
print('x_senales: ')
for i in range(x_conteo):
    print('senal:  ',x_senales[i])
    print('  euler:',Xe_senales[i])
print('x_simplif: ')
print(' ',x_simplif)
print('x_espectro:')
for unparam in Xe_espectro:
    print(unparam,':',Xe_espectro[unparam])

[ producto cos ] [ ejercicio ] [ analítico ] [ algoritmo ] [gráfica ] [ AM y BW]
..


5. Gráfica de espectro de frecuencias

espectro Senales Producto Cos 01

Instrucciones adicionales para la gráfica de espectro son las mismas que en la sección anterior

# GRAFICAS de espectro de frecuencias ---------
freq = Xe_espectro['freq']
magnitud = Xe_espectro['ampl']
etiqueta = Xe_espectro['etiq']
mag_max = float(max(magnitud))
freq_max = float(max(freq))

# grafica 
graf_dx = 0.12
fig_espectro = plt.figure()
graf_fasor = fig_espectro.add_subplot()
# grafica magnitud
graf_fasor.set_xlim([-freq_max*(1+graf_dx),freq_max*(1+graf_dx)])
graf_fasor.set_ylim([0,mag_max*(1+graf_dx*2)])
graf_fasor.axhline(0,color='black')
graf_fasor.axvline(0,linestyle='dotted',color='grey')
graf_fasor.stem(freq,magnitud)
# etiquetas de fasor
for k in range(0,len(freq),1):
    texto = etiqueta[k]
    x_text = freq[k]
    y_text = magnitud[k]
    plt.annotate(texto,xy=(x_text,y_text), xytext=(0,5),
                 textcoords='offset points',ha='center')
graf_fasor.grid()
graf_fasor.set_xlabel('freq Hz')
graf_fasor.set_ylabel('amplitud')
graf_fasor.set_title('Espectro: x_senales')

#plt.show()

La gráfica de tiempo, permite observar el efecto del producto se señales.

espectro Senales Producto Cos t 02

Se reutilizan las instrucciones de la sección 1.3, añadiendo las variantes de estilo de línea para destacar la señal resultante del producto.

En este caso se está usando la lista de frecuencias del espectro de señales, lo que no incluye las frecuencias individuales de cada componente al desaparecer en el resultado las frecuencia originales. por lo que la «envolvente» de menor frecuencia no domina sobre las demás.
Considerar el caso presentado para completar el algoritmo si se requiere mayor precisión.

# GRAFICAR x(t) ---------------------------
# INGRESO
x_senales_graf = [x1,x2,x] # para graficar
x_etiqueta = ['x1(t)','x2(t)','x1(t)x2(t)']
x_tipolinea = ['dashed','dotted','solid']

tramos_T = 20       # tramos por periodo
periodos_graf = 2   # periodos en una grafica

# PROCEDIMIENTO GRAFICA
x_conteo_graf = len(x_senales_graf)
# etiquetas si están vacias
if len(x_etiqueta)==0: 
    if x_conteo_graf > 1:
        for i in range(0,x_conteo_graf-1,1):
            x_etiqueta.append('x'+str(i)+'(t)')
            x_tipolinea.append('dotted')
    x_etiqueta.append('x(t)')
    x_tipolinea.append('solid')
    
[T_min,T_max] = dsp.periodo_minmax(freq)
x_conteo = len(x_senales_graf)

muestras_T = tramos_T + 1
# intervalo de t entre [a,b] en pasos dt
a = 0
b = T_max*periodos_graf
dt = T_min/tramos_T # tamaño de paso,tramo, muestreo
ti = np.arange(a,b+dt,dt)

xi_k = [] # muestras de cada señal x
for unasenal in x_senales_graf:
    # x(t) a forma numérica lambda t:
    if unasenal.has(t):
        xk = sym.lambdify(t,unasenal)
    else: # es constante
        xk = lambda t: unasenal + 0*t
    xki = xk(ti) # evalúa en intervalo
    xi_k.append(np.copy(xki))

# graficar
fig_x = plt.figure()
graf_x = fig_x.add_subplot()
for i in range(0,x_conteo_graf,1):
    color_i = dsp._graf_lineacolor_i(i)
    graf_x.plot(ti,xi_k[i], color=color_i,
                linestyle=x_tipolinea[i],
                label=x_etiqueta[i])
    graf_x.axhline(0,color='black')
graf_x.set_title('Señales xk(t)')
graf_x.set_xlabel('t (segundos)')
graf_x.set_ylabel('x_senales')
graf_x.grid()
graf_x.legend()
plt.show()

[ producto cos ] [ ejercicio ] [ analítico ] [ algoritmo ] [gráfica ] [ AM y BW]
..


6. Modulación de Amplitud y Ancho de Banda

En comunicación por radio difusión se conoce el término Modulación de amplitud AM (amplitude modulation). AM es el proceso de multiplicar una señal sinusoidal de alta frecuencia fc (portadora) con una señal de baja frecuencia como la voz o música (envolvente)

x(t) = envolvente (portadora)

x(t) = v(t) \cos ( 2 \pi f_c t)

El concepto de ancho de banda se usa con la representación del espectro de frecuencias de la señal y destaca una de sus características. Para el caso de AM cuyo espectro tiene como frecuencia mas alta (fc + fΔ ) y referenciada a la frecuencia de la portadora fc , se considera importante el grupo de frecuencias positivas  en el intervalo entre la mas alta y las mas baja conocidas como ancho de banda o BW (Bandwitdth) de la señal.

(f_c + f_{\Delta})-(f_c - f_{\Delta}) = 2f_{\Delta}

6.1 Ejercicio

Referencia: McClellan ejemplo 3.4 p78 y ejercicio 3.3 p80

Realice el espectro de frecuencias para una señal AM con portadora cos(400πt) y envolvente (5+4cos(40 πt)).

La señal envolvente tiene un término DC suficientemente grande que la mantiene en la parte superior (positiva)  para todo ‘t’. El término DC de la envolvente permite simplificar la implementación de un circuito detector de envolvente en el radio-receptor, la envolvente es la voz o música que se quiere escuchar.

Obseve  que la señal de la portadora es de 200Hz y la envolvente es 20Hz, con una relacion 1/10. En el espectro de frecuencia, la mas alta es 220 Hz y la mas baja es 180 Hz, la frecuencia central o portadora es 200 Hz. El ancho de banda es (220-180) = 40 Hz. que es la freciencia de la envolvente.

[ producto cos ] [ ejercicio ] [ analítico ] [ algoritmo ] [gráfica ] [ AM y BW]

6.2 Algoritmo con Python

A partir del algoritmo del producto de señales,  se añade el cálculo de las frecuencia máxima y mínima de x(t).

# ancho de banda
freq = Xe_espectro['freq']
freq_max = float(max(freq))
freq_min = float(min(freq[freq>0])) # freq positivas [freq>0]
freq_centro = (freq_max+freq_min)/2
BW = freq_max-freq_min

con lo que se obtiene como resultado

x(t): (4*cos(40*pi*t) + 5)*cos(400*pi*t)
x_senales: 
senal:   5*cos(400*pi*t)
  euler: 5*exp(400*I*pi*t)/2 + 5*exp(-400*I*pi*t)/2
senal:   4*cos(40*pi*t)*cos(400*pi*t)
  euler: exp(440*I*pi*t) + exp(360*I*pi*t) + exp(-360*I*pi*t) + exp(-440*I*pi*t)
x_simplif: 
  2*cos(360*pi*t) + 5*cos(400*pi*t) + 2*cos(440*pi*t)
   freq_min: 180.0
freq_centro: 200.0
   freq_max: 220.0
ancho de banda BW: 40.0
x_espectro:
freq : [-220. -200. -180.  180.  200.  220.]
ampl : [1 5/2 1 1 5/2 1]
fase : [0 0 0 0 0 0]
etiq : ['1' '5/2' '1' '1' '5/2' '1']

con gráficas de espectro de frecuencia y x(t)

espectro Senales Producto AM 01.

Senales Producto AM 02

Algoritmo en Python con gráficas

# ejemplo 3.4 p78 Modulación de Amplitud (AM) o Amplitud Modulada
# telg1034 DSP fiec-espol edelros@espol.edu.ec
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
import telg1034 as dsp

# variables
from telg1034 import t,A,w,f,p,pi,DosPi,I,equivalentes

# INGRESO
x1 = 5+4*sym.cos(40*pi*t)
x2 = sym.cos(400*pi*t)
x = x1*x2

# PROCEDIMIENTO
x_senales = dsp.cos_x_components_add(x)
x_conteo = len(x_senales)
Xe_senales = dsp.cos_to_euler_one_term(x_senales)
Xe_espectro = dsp.cos_spectrum_list(x_senales)

# x_senales simplificados
Xe_conteo = len(Xe_senales)
x_senales_simple = []
x_simplif = sym.S.Zero
for i in range(Xe_conteo):
    unaSenal = Xe_senales[i]
    Xe_to_cos = unaSenal.rewrite(sym.cos)
    x_simplif = x_simplif + Xe_to_cos

# ancho de banda
freq = Xe_espectro['freq']
freq_max = float(max(freq))
freq_min = float(min(freq[freq>0])) # freq positivas [freq>0]
freq_centro = (freq_max+freq_min)/2
BW = freq_max-freq_min

# SALIDA
print('x(t):',x)
print('x_senales: ')
for i in range(x_conteo):
    print('senal:  ',x_senales[i])
    print('  euler:',Xe_senales[i])
print('x_simplif: ')
print(' ',x_simplif)
print('   freq_min:',freq_min)
print('freq_centro:',freq_centro)
print('   freq_max:',freq_max)
print('ancho de banda BW:',BW)
print('x_espectro:')
for unparam in Xe_espectro:
    print(unparam,':',Xe_espectro[unparam])

# GRAFICAS de espectro de frecuencias ---------
freq = Xe_espectro['freq']
magnitud = Xe_espectro['ampl']
etiqueta = Xe_espectro['etiq']
mag_max = float(max(magnitud))
freq_max = float(max(freq))

# grafica 
graf_dx = 0.12
fig_espectro = plt.figure()
graf_fasor = fig_espectro.add_subplot()
# grafica magnitud
graf_fasor.set_xlim([-freq_max*(1+graf_dx),freq_max*(1+graf_dx)])
graf_fasor.set_ylim([0,mag_max*(1+graf_dx*2)])
graf_fasor.axhline(0,color='black')
graf_fasor.axvline(0,linestyle='dotted',color='grey')
graf_fasor.stem(freq,magnitud)
# etiquetas de fasor
for k in range(0,len(freq),1):
    texto = etiqueta[k]
    x_text = freq[k]
    y_text = magnitud[k]
    plt.annotate(texto,xy=(x_text,y_text), xytext=(0,5),
                 textcoords='offset points',ha='center')
graf_fasor.grid()
graf_fasor.set_xlabel('freq Hz')
graf_fasor.set_ylabel('amplitud')
graf_fasor.set_title('Espectro: x_senales')

#plt.show()

# GRAFICAR x(t) ---------------------------
# INGRESO
x_senales_graf = [x1,x2,x_simplif] # para graficar
x_etiqueta = ['x1(t)','x2(t)','x(t)']
x_tipolinea = ['dashed','dotted','solid']

tramos_T = 20       # tramos por periodo
periodos_graf = 4   # periodos en una grafica

# PROCEDIMIENTO GRAFICA
x_conteo_graf = len(x_senales_graf)
# etiquetas si están vacias
if len(x_etiqueta)==0: 
    if x_conteo_graf > 1:
        for i in range(0,x_conteo_graf-1,1):
            x_etiqueta.append('x'+str(i)+'(t)')
            x_tipolinea.append('dotted')
    x_etiqueta.append('x(t)')
    x_tipolinea.append('solid')
    
[T_min,T_max] = dsp.periodo_minmax(freq)
x_conteo = len(x_senales_graf)

muestras_T = tramos_T + 1
# intervalo de t entre [a,b] en pasos dt
a = 0
b = T_max*periodos_graf
dt = T_min/tramos_T # tamaño de paso,tramo, muestreo
ti = np.arange(a,b+dt,dt)

xi_k = [] # muestras de cada señal x
for unasenal in x_senales_graf:
    # x(t) a forma numérica lambda t:
    if unasenal.has(t):
        xk = sym.lambdify(t,unasenal)
    else: # es constante
        xk = lambda t: unasenal + 0*t
    xki = xk(ti) # evalúa en intervalo
    xi_k.append(np.copy(xki))

# graficar
fig_x = plt.figure()
graf_x = fig_x.add_subplot()
for i in range(0,x_conteo_graf,1):
    color_i = dsp._graf_lineacolor_i(i)
    graf_x.plot(ti,xi_k[i], color=color_i,
                linestyle=x_tipolinea[i],
                label=x_etiqueta[i])
    graf_x.axhline(0,color='black')
graf_x.set_title('Señales xk(t)')
graf_x.set_xlabel('t (segundos)')
graf_x.set_ylabel('x_senales')
graf_x.grid()
graf_x.legend()
plt.show()

[ producto cos ] [ ejercicio ] [ analítico ] [ algoritmo ] [gráfica ] [ AM y BW]

1.4 Espectro – Suma de sinusoides y la fórmula de Euler con Sympy-Python

[ espectro ] [ ejercicio ] [ analítico ] [ algoritmo ] [ gráfica ]
..


1. Espectro de una suma de sinusoides

Referencia: McClellan 3.1 p70

El método mas general para construir nuevas señales a partir de sinusoides es la combinación lineal, donde una señal se compone de la suma de una constante y N sinusoides, cada una con diferente frecuencia, amplitud y fase.

x(t) = A_0 + \sum_{k =1}^{N} A_k \cos ( 2 \pi f_k t + \varphi_k)

que como exponenciales complejos

x(t) = X_0 + \sum_{k =1}^{N} \Re \Big\{ X_k e^{ j 2 \pi f_k t } \Big\} X_k = A_k e^{ \varphi_k }

donde X0 = A0 corresponde al componente de una constante real.

La fórmula inversa de Euler permite escribir la señal x(t) como:

x(t) = X_0 + \sum_{k =1}^{N} \Re \Big\{ \frac{X_k}{2} e^{ j 2 \pi f_k t} + \frac{X_k^*}{2} e^{ -j 2 \pi f_k t} \Big\}

[ espectro ] [ ejercicio ] [ analítico ] [ algoritmo ] [ gráfica ]
..


2. Ejemplo – Espectro de dos lados

Referencia: McClellan ejemplo 3.1 p71

Determine el espectro de la siguiente señal.

x(t) = 10 + 14 \cos(200 \pi t - \pi/3) + 8 cos(500 \pi t + \pi/2)

[ espectro ] [ ejercicio ] [ analítico ] [ algoritmo ] [ gráfica ]
..


3. Desarrollo analítico

La expresión se separa en los componentes de cada término de la suma

x_1(t) = 10 x_2(t) = 14 \cos(200 \pi t - \pi/3) x_3(t) = 8 cos(500 \pi t + \pi/2)

se escribe cada término en su forma exponencial

x_1(t) = 10 e^{j0t} = 10 x_2(t) = 7 e^{- j \pi/3} e^{j2\pi(100) t} + 7 e^{j \pi/3} e^{-j 2 \pi (100) t} x_3(t) = 4 e^{- j \pi/2} e^{j 2\pi(250) t} + 4 e^{j \pi/2} e^{-j 2 \pi (250) t}

Para la gráfica de espectro de señal se requiere frecuencia, amplitud y fasor, que pueden construir con los parámetros de cada señal usando la función  dsp.cos_args_one_term(señal).

El resultado del algoritmo para el ejercicio es:

x_senales: 
senal:   10
  euler: 10
senal:   -8*sin(500*pi*t)
  euler: 4*I*exp(500*I*pi*t) - 4*I*exp(-500*I*pi*t)
senal:   14*sin(200*pi*t + pi/6)
  euler: 7*exp(-I*pi/3)*exp(200*I*pi*t) + 7*exp(I*pi/3)*exp(-200*I*pi*t)
x_espectro:
freq : [-250. -100.    0.  100.  250.]
ampl : [4 7 10 7 4]
fase : [-pi/2 pi/3 0 -pi/3 pi/2]
etiq : ['4$ e^j(- \\frac{\\pi}{2})$' '7$ e^j(\\frac{\\pi}{3})$' 10
 '7$ e^j(- \\frac{\\pi}{3})$' '4$ e^j(\\frac{\\pi}{2})$']
>>>

espectro de señales por frecuencia

[ espectro ] [ ejercicio ] [ analítico ] [ algoritmo ] [ gráfica ]
..


4. Algoritmo con Python

4.1 Señal suma de términos

Para este ejercicio, la señal es una suma de varios componentes: constantes y senoidales.

# INGRESO
x0 = 10
x1 = 14*sym.cos(200*pi*t - pi/3)
x2 = 8*sym.cos(500*pi*t + pi/2)
x = x0 + x1 + x2

Para simplificar el análisis como en los ejercicios anteriores se separan los términos de la suma en una lista x_senales, al revisar si x es una suma x.is_Add:

# Revisa expresión x(t)
x = sym.sympify(x) # expresión a sympy, por si es constante
x = sym.expand(x) # simplific parentesis (A + 3)*(cos(2*pi*f*t +w))
if x.is_Add: # separa términos suma
    x_senales = list(x.args)    
else:
    x_senales = [x]

Una señal de la x_senales se convierte al formato Euler con las instrucciones

Xj = unasenal.rewrite(sym.exp)  # Xj con forma Euler
Xj = sym.expand(Xj)

Por ejemplo, si se usa x1(t), el resultado es la amplitud por la expresión Euler,

>>> x1
14*sin(200*pi*t + pi/6)
>>> x1.rewrite(sym.exp)
-7*I*(-exp(I*(-200*pi*t - pi/6)) + exp(I*(200*pi*t + pi/6)))

para disponer de términos suma, se simplifica la expresión con:

>>> sym.expand(x1.rewrite(sym.exp))
-7*I*exp(I*pi/6)*exp(200*I*pi*t) + 7*I*exp(-I*pi/6)*exp(-200*I*pi*t)
>>>

Si se aplican las instrucciones dadas a cada señal en x_senales, se obtiene X_senales:

x_senales: 
senal:   10
  euler: 10
senal:   -8*sin(500*pi*t)
  euler: 4*I*exp(500*I*pi*t) - 4*I*exp(-500*I*pi*t)
senal:   14*sin(200*pi*t + pi/6)
  euler: 7*exp(-I*pi/3)*exp(200*I*pi*t) + 7*exp(I*pi/3)*exp(-200*I*pi*t)

Sympy simplifica automáticamente e(j π/2) como I y  e(-j π/2) como -I

Las instrucciones en Python quedan como:

# ejemplo 3.1 p71 x_senales a formato_euler
# telg1034 DSP fiec-espol edelros@espol.edu.ec
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
import telg1034 as dsp

# variables
from telg1034 import t,A,w,f,p,pi,DosPi,I,equivalentes

# INGRESO
x0 = 10
x1 = 14*sym.cos(200*pi*t - pi/3)
x2 = 8*sym.cos(500*pi*t + pi/2)
x = x0 + x1 + x2

# PROCEDIMIENTO
# Revisa expresión x(t)
x = sym.sympify(x) # expresión a sympy, por si es constante
x = sym.expand(x) # simplific parentesis (A + 3)*(cos(2*pi*f*t +w))
if x.is_Add: # separa términos suma
    x_senales = list(x.args)    
else:
    x_senales = [x]
# Analiza x_senales
x_conteo = len(x_senales)
X_senales = []
for i in range(0,x_conteo,1):
    unasenal = x_senales[i]
    unasenal = sym.sympify(unasenal) # expresión a sympy, por si es constante
    cond1 = unasenal.has(t)
    cond2 = unasenal.has(sym.cos) or unasenal.has(sym.sin)
    if cond1 and cond2: # tiene variable t y es senoidal
        if unasenal.has(sym.sin): # evita sin()
            unasenal = unasenal.rewrite(sym.cos)
        Xj = unasenal.rewrite(sym.exp)  # Xj con forma Euler
        Xj = sym.expand(Xj)
    else: # es una constante
        Xj = unasenal
    X_senales.append(Xj)

# SALIDA
print('x_senales: ')
for i in range(x_conteo):
    print('senal:  ',x_senales[i])
    print('  euler:',X_senales[i])

4.2. Espectro de frecuencia de dos lados

Para realizar la gráfica de espectro de frecuencias, se requiere obtener las listas de frecuencia, amplitud, fase y etiquetas en un diccionario con las listas de valores correspondientes. La creación de la gráfica se simplifica al usar agrupar las instrucciones en una función denominada cos_spectrum_list(x_senales), para obtener los resultados según lo indicado.

Los argumentos de los términos Euler se interpretan como frecuencia, amplitud y fase, se crea euler_args_one_term(X) a partir de lo realizado para la función coseno. Considerando que para la fase es el factor multiplicador con exp() sin variable t, puede ajustarse si aparece un factor con I, o puede considerar el signo en -I.

También, los procedimientos para convertir las señales a la forma Euler, se convierten en funciones de un solo término, cos_to_euler_one_term(x_senales), verificando que cada uno sea SinSumas. En caso de tener sumas se muestra un mensaje que indica los componentes de x_senales que requieran corregirse.

El resultado del algoritmo para el ejercicio es:

x_senales: 
senal:   10
  euler: 10
senal:   -8*sin(500*pi*t)
  euler: 4*I*exp(500*I*pi*t) - 4*I*exp(-500*I*pi*t)
senal:   14*sin(200*pi*t + pi/6)
  euler: 7*exp(-I*pi/3)*exp(200*I*pi*t) + 7*exp(I*pi/3)*exp(-200*I*pi*t)
x_espectro:
freq : [-250. -100.    0.  100.  250.]
ampl : [4 7 10 7 4]
fase : [-pi/2 pi/3 0 -pi/3 pi/2]
etiq : ['4$ e^j(- \\frac{\\pi}{2})$' '7$ e^j(\\frac{\\pi}{3})$' 10
 '7$ e^j(- \\frac{\\pi}{3})$' '4$ e^j(\\frac{\\pi}{2})$']
>>>

Instrucciones en Python

# ejemplo 3.1 p71 x_senales a formato_euler
# telg1034 DSP fiec-espol edelros@espol.edu.ec
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
import telg1034 as dsp

# variables
from telg1034 import t,A,w,f,p,pi,DosPi,I,equivalentes

# INGRESO
x0 = 10
x1 = 14*sym.cos(200*pi*t - pi/3)
x2 = 8*sym.cos(500*pi*t + pi/2)
x = x0 + x1 + x2

# PROCEDIMIENTO
def cos_x_components_add(x):
    ''' x(t) separa términos como suma de componentes,
    simplifica paréntesis
    '''
    x = sym.sympify(x) # expresión a sympy, por si es constante
    x = sym.expand(x) # simplifica paréntesis (A + 3)*(cos(2*pi*f*t +w))
    if x.is_Add: # separa términos suma
        x_senales = list(x.args)    
    else:
        x_senales = [x]
    return(x_senales)
    
def cos_to_euler_one_term(x_senales):
    ''' convierte lista x_senales a la forma Euler
    entregando la lista X_senales de cada señal en forma Euler
    '''
    x_conteo = len(x_senales)
    X_senales = []
    SinSumas = True # Analizar un solo cos o sin
    SinSumas_cual = []
    for i in range(0,x_conteo,1):
        unasenal = x_senales[i]
        unasenal = sym.sympify(unasenal) # expresión a sympy, por si es constante
        if unasenal.is_Add:
            SinSumas = False
        cond1 = unasenal.has(t)
        cond2 = unasenal.has(sym.cos) or unasenal.has(sym.sin)
        if cond1 and cond2 and SinSumas: # tiene variable t y es senoidal
            if unasenal.has(sym.sin): # evita sin()
                unasenal = unasenal.rewrite(sym.cos)
            Xj = unasenal.rewrite(sym.exp)  # Xj con forma Euler
            Xj = sym.expand(Xj)
        else: # es una constante
            Xj = unasenal
        if not(SinSumas):
            SinSumas_cual.append(i)
        X_senales.append(Xj)
    if SinSumas == False:
        SinSumas_conteo = len(SinSumas_cual)
        msg ='x_senales tiene '+str(SinSumas_conteo)+'suma de terminos,'
        msg = msg + 'usar solo un término de la forma:'
        msg = msg + '\n A*cos(w*t+p) o también A*sin(w*t+p) ... \n'
        print('cos_to_euler_one_term:',x,msg)
        for i in SinSumas_cual:
            print('  revisar: ',x_senales[i])
    return(X_senales)

def cos_spectrum_list(x_senales):
    '''dado una lista de señales de un solo término,
    se obtiene la lista de frecuencias, amplitudes y etiquetas
    para crear la gráfica de espectro de frecuencias
    '''
    # Analiza x_senales
    x_conteo = len(x_senales)
    x_freq_spectr = [] ; x_ampl_spectr = []
    x_etiq_spectr = [] ; x_fase_spectr = []
    for i in range(0,x_conteo,1):
        unasenal = x_senales[i]
        X_senal = unasenal.rewrite(sym.exp)
        X_senal = sym.expand(X_senal)
        if X_senal.is_Add:
            X_senales = X_senal.args
        else:
            X_senales = [X_senal]
        for un_X in X_senales:
            parametro = euler_args_one_term(un_X)
            freq_k = float(parametro['freq_Hz'])
            ampl_k = parametro['amplitud']
            fase_k = parametro['fase']
            x_freq_spectr.append(freq_k)
            x_ampl_spectr.append(ampl_k)
            x_fase_spectr.append(fase_k)
            texto = str(ampl_k)
            if fase_k!=sym.S.Zero:
                texto = texto+f'$ e^j('+sym.latex(fase_k)+')$'
            x_etiq_spectr.append(texto)
    # ordenar por freq_Hz
    x_freq_spectr = np.array(x_freq_spectr)
    x_ampl_spectr = np.array(x_ampl_spectr)
    x_etiq_spectr = np.array(x_etiq_spectr)
    x_fase_spectr = np.array(x_fase_spectr)
    orden = np.argsort(x_freq_spectr)
    x_espectro = {'freq':x_freq_spectr[orden],
                  'ampl':x_ampl_spectr[orden],
                  'fase':x_fase_spectr[orden],
                  'etiq':x_etiq_spectr[orden],}
    return(x_espectro)

def euler_args_one_term(X):
    ''' versión 1. Amplitud, frecuencia en Hz y rad, fase de exp()
    referenciado a un término coseno.
    Entrega diccionario con 'amplitud', 'freq_Hz','freq_rad','fase','periodo'.
    '''
    parametro={} # parámetros en diccionario
    amplitud = sym.S.One ; T = sym.S.Zero ; fase =  sym.S.Zero
    freq_Hz = sym.S.Zero ; freq_rad =  sym.S.Zero

    X = sym.sympify(X) # expresión a sympy, por si es constante
    X = sym.expand(X)  # expresión con al menos una suma
    SinSumas = True # Analizar un solo exp()
    if X.is_Add:
        SinSumas = False

    if SinSumas: # un solo termino: amplitud*exp(Iw*t+I*p)
        #X = sym.powsimp(X)
        partes = X.args
        if len(partes)==0: # amplitud es una constante.
            amplitud = X  # no tiene args
        for unaparte in partes:
            cond1 = unaparte.has(sym.exp)
            cond2 = unaparte.has(t)
            if cond1 and cond2: # exp(I*w*t+I*p)
                argumento_exp = sym.expand(unaparte.args[0]/I)
                [freq_Hz,freq_rad,fase_k] = dsp._cos_arg_freqfase(argumento_exp)
                fase = fase + fase_k
            if cond1 and not(cond2): # exp(p*I)
                argumento_exp = sym.expand(unaparte.args[0]/I)
                fase = fase+argumento_exp
            if not(cond1) and not(cond2): # amplitud
                if not(unaparte.has(I)):
                    amplitud = amplitud*unaparte
                else: # amplitud*I o -amplitud*i
                    amplitud = sym.expand(amplitud*unaparte/I)
                    fase = fase + pi/2
                    if amplitud<0:
                        amplitud = abs(amplitud)
                        fase = fase - pi
            if not(cond1) and cond2: # I*w*t+I*p
                argumento_exp = sym.expand(unaparte/I)
                [freq_Hz,freq_rad,fase] = dsp._cos_arg_freqfase(argumento_exp)
                
        # amplitud, revisión numérica
        if amplitud.is_Number:
            if amplitud<0: # amplitud con signo positivo
                amplitud = np.abs(amplitud)
                fase = fase + sym.pi
            if isinstance(amplitud,sym.Float):
                amplitud = float(amplitud)
            if isinstance(amplitud,sym.Integer):
                amplitud = int(amplitud)
            # mantiene forma racional (a/b) de amplitud
        # periodo T
        if freq_Hz != sym.S.Zero:
            T = 1/freq_Hz
        else: # es una constante
            T = sym.S.NaN
    if not(X.has(t)): # es una constante
        amplitud = X
    # parametro fasor y fasor_complex
    fasor = amplitud*sym.exp(I*fase)
    fasor_complex = sym.expand(fasor,complex=True).evalf()
    parametro['amplitud'] = amplitud
    parametro['freq_Hz']  = freq_Hz
    parametro['freq_rad'] = freq_rad
    parametro['fase'] = fase
    parametro['periodo'] = T
    parametro['fasor'] = fasor
    parametro['fasor_complex'] = fasor_complex
                
    if SinSumas == False:
        msg ='x(t) tiene suma de terminos, usar solo un término de la forma:'
        msg = msg + '\n A*exp(I*(w*t+p)) o también A*exp(-I*(w*t+p)) ... \n'
        print('euler_args_one_term:',x,msg)
    return(parametro)

# operaciones con señal x(t)
x_senales = cos_x_components_add(x)
x_conteo = len(x_senales)
X_senales = cos_to_euler_one_term(x_senales)
X_espectro = cos_spectrum_list(x_senales)

# SALIDA
print('x_senales: ')
for i in range(x_conteo):
    print('senal:  ',x_senales[i])
    print('  euler:',X_senales[i])
print('x_espectro:')
for unparam in X_espectro:
    print(unparam,':',X_espectro[unparam])

[ espectro ] [ ejercicio ] [ analítico ] [ algoritmo ] [ gráfica ]
..


5. Gráfica de espectro de frecuencias

La gráfica de espectro de frecuencias se obtiene con los datos de x_espectro realizado en el bloque anterior. Se incorpora cada fasor por frecuencia como las etiquetas creada por la función cos_spectrum_list(x_senales).

espectro de señales por frecuencia

Instrucciones adicionales para la gráfica

# GRAFICAS de espectro de frecuencias ---------
freq = X_espectro['freq']
magnitud = X_espectro['ampl']
etiqueta = X_espectro['etiq']
mag_max = float(max(magnitud))
freq_max = float(max(freq))

# grafica 
graf_dx = 0.12
fig_espectro = plt.figure()
graf_fasor = fig_espectro.add_subplot()
# grafica magnitud
graf_fasor.set_xlim([-freq_max*(1+graf_dx),freq_max*(1+graf_dx)])
graf_fasor.set_ylim([0,mag_max*(1+graf_dx*2)])
graf_fasor.axhline(0,color='black')
graf_fasor.axvline(0,linestyle='dotted',color='grey')
graf_fasor.stem(freq,magnitud)
# etiquetas de fasor
for k in range(0,len(freq),1):
    texto = etiqueta[k]
    x_text = freq[k]
    y_text = magnitud[k]
    plt.annotate(texto,xy=(x_text,y_text), xytext=(0,5),
                 textcoords='offset points',ha='center')
graf_fasor.grid()
graf_fasor.set_xlabel('freq Hz')
graf_fasor.set_ylabel('amplitud')
graf_fasor.set_title('Espectro: x_senales')

plt.show()

[ espectro ] [ ejercicio ] [ analítico ] [ algoritmo ] [ gráfica ]

1.3 Señal – Suma por Fasores de igual frecuencia con Sympy-Python

[ fasores ] [ ejemplo ] [ analítico ] [ algoritmo ] [ gráfica] [ función ]
..


1. Suma por Fasores

Referencia: McClellan 2.6 p48

Con señales de igual frecuencia en el mismo medio de transmisión las señales se suman. Una forma simplificada de sumar señales con diferente amplitud, fase e igual frecuencia es por medio de fasores , considerando:

x(t) = \sum_{k=1}^{N} A_k \cos(\omega_0 t + \varphi_k) = A \cos(\omega_0 t + \varphi)

Se puede representar los cosenos en forma compleja

A \cos(\omega_0 t + \varphi) = \Re \Big\{ A e^{j(\omega_0 t + \varphi)} \Big\} = \Re \Big \{ A e^{j\varphi} e^{j \omega_0 t} \Big\} \sum_{k=1}^{N} A_k e^{j\varphi_k} = A e^{j \varphi}

[ fasores ] [ ejemplo ] [ analítico ] [ algoritmo ] [ gráfica] [ función ]
..


2. Ejemplo – regla de suma de fasores

Referencia: McClellan 2.6.3 p51

Sumar dos señales sinusoidales de igual frecuencia usando fasores:

x_1 (t) = 1.7\cos(20 \pi t + 70\pi / 180) x_2 (t) = 1.9\cos(20 \pi t + 200\pi / 180)

[ fasores ] [ ejemplo ] [ analítico ] [ algoritmo ] [ gráfica] [ función ]
..


3. Desarrollo analítico

La frecuencia es 10 Hz para ambas señales.
Las suma de x1(t) + x2(t) mediante fasores tiene cuatro pasos:

1. cada señal se escribe usando fasores:

X_1 = 1.7e^{j 70\pi / 180} X_2 = 1.9e^{j 200\pi / 180}

2. cada fasor se escribe en forma compleja (rectangular):

X_1 = 0.5814 + j 1.5975 X_2 = -1.7854 - j 0.6498

3. se suman las partes reales e imaginarias

X_3 = X_1 + X_2 = (0.5814 + j 1.5975) + (-1.7854 - j 0.6498) = -1.204 + j 0.9477

4.Se reescribe a la forma polar y luego a senoidal con la frecuencia de las señales

X_3 = 1.5322 e^{j 141.79\pi/180} x_3(t) = 1.5322\cos(20\pi t +141.79\pi /180)

[ fasores ] [ ejemplo ] [ analítico ] [ algoritmo ] [ gráfica] [ función ]
..


4. Algoritmo con Sympy

Las señales en el bloque de ingreso son:

x1 = 1.7*sym.cos(20*pi*t+70*pi/180)
x2 = 1.9*sym.cos(20*pi*t+200*pi/180)

Desarrollando solo para x1(t) se obtienen los parámetros de un término coseno,

x1_args = dsp.cos_args_one_term(x1)
amplitud = x1_args['amplitud']
fase = x1_args['fase']

El fasor se obtiene escribiendo las expresiones con amplitud y fase,

fasor1 = amplitud*sym.exp(I*fase)

para la forma de números complejos, se usa sym.expand(fasor1,complex=True),

>>> fasor1
1.7*exp(7*I*pi/18)
>>> sym.expand(fasor1,complex=True)
1.7*cos(7*pi/18) + 1.7*I*sin(7*pi/18)
>>> sym.expand(fasor1,complex=True).evalf()
0.581434243653637 + 1.59747745533604*I
>>> 

sin embargo la expresión muestra el cos() para la parte real y sin() para la imaginaria. La expresión las evalúa la funciones senoidales añadiendo .evalf() al final de la instrucción, quedando la instrucción de la siguiente forma:

fasor1_complex = sym.expand(fasor1,complex=True).evalf()

El proceso se repite para la señal x2(t) para luego sumar los fasores en forma compleja, siempre y cuando las frecuencias sean iguales.

Para mantener la amplitud del fasor_suma con signo positivo, evitando que Sympy simplifique la fase con π se usa la instrucción sym.UnevaluatedExpr(fase).

El resultado del algoritmo se muestra a continuación:

fasor1: 1.7*exp(7*I*pi/18)
fasor1_complex: 0.581434243653637 + 1.59747745533604*I
fasor2: 1.9*exp(-8*I*pi/9)
fasor2_complex: -1.78541597949323 - 0.649838272318771*I
fasor_suma_complex: -1.20398173583959 + 0.947639183017274*I
fasor_suma: 1.53218538089389*exp(-0.666817831701547 + pi)
x_suma: 1.53218538089389*cos(10*t*(2*pi) - (0.666817831701547 + pi))
>>>

Las instrucciones en Python son:

# ejemplo 2.6.3 p52 Suma de Fasores y graf_polar
# telg1034 DSP fiec-espol edelros@espol.edu.ec
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
import telg1034 as dsp

# variables continuas
from telg1034 import t,A,w,f,p,pi,DosPi,I,equivalentes

# INGRESO
x1 = 1.7*sym.cos(20*pi*t+70*pi/180)
x2 = 1.9*sym.cos(20*pi*t+200*pi/180)

# PROCEDIMIENTO
# señal x1
x1_args = dsp.cos_args_one_term(x1)
amplitud = x1_args['amplitud']
fase = x1_args['fase']
fasor1 = amplitud*sym.exp(I*fase)
# fasor números complejos
fasor1_complex = sym.expand(fasor1,complex=True).evalf()

# señal x2
x2_args = dsp.cos_args_one_term(x2)
amplitud = x2_args['amplitud']
fase = x2_args['fase']
fasor2 = amplitud*sym.exp(I*fase)
# fasor números complejos
fasor2_complex = sym.expand(fasor2,complex=True).evalf()

# suma de fasores de x1(t) y x2(t) en la misma freq
if x1_args['freq_Hz'] == x2_args['freq_Hz']:
    freq_iguales = True
    fasorsum_complex = fasor1_complex + fasor2_complex
    
    amplitud = sym.Abs(fasorsum_complex)
    fase = sym.arg(fasorsum_complex)
        
    fasor_suma = amplitud*sym.exp(sym.UnevaluatedExpr(fase))
    x_suma = amplitud*sym.cos(DosPi*x1_args['freq_Hz']*t+sym.UnevaluatedExpr(fase))
    x_suma_args = dsp.cos_args_one_term(x_suma)
    
else: # diferentes frecuencias
    freq_iguales = False
    fasorsum_complex = 0+0j
    fasor_suma = fasor1 +fasor2
    x_suma = x1+x2

# SALIDA
print('fasor1:',fasor1)
print('fasor1_complex:',fasor1_complex)
print('fasor2:',fasor2)
print('fasor2_complex:',fasor2_complex)
if freq_iguales:
    print('fasor_suma_complex:',fasorsum_complex)
    print('fasor_suma:',fasor_suma)
else:
    print('  las frecuencias no son iguales en x1(t) y x2(t)')
print('x_suma:',x_suma)

[ fasores ] [ ejemplo ] [ analítico ] [ algoritmo ] [ gráfica] [ función ]
..


5. Gráfica de fasores

La gráfica de fasores es de tipo polar con flechas que parten desde el centro y se dibujan usando la amplitud y fase por cada fasor.

Primero se genera la lista de x_amplitud y x_fase, la suma se añade solamente si las señales son de igual frecuencia.

La forma polar de la gráfica se indica con projection='polar' al definir el sub-gráfico. Como referencia se añade puntos al final de cada fasor y se diferencian al usar un color diferente.

La instrucción de flechas, requiere el punto de partida (0,0) y el punto de llegada para cada fasor (x_fase[i],x_amplitud[i]).

Las etiquetas del gráfico se muestran en radianes
Suma fasores igual frecuencia 01

Instrucciones adicionales para la gráfica polar

# GRAFICA de suma de fasores
x_etiqueta = ['x1(t)','x2(t)']
x_amplitud = [float(x1_args['amplitud']),
              float(x2_args['amplitud'])]
x_fase = [float(x1_args['fase']),
          float(x2_args['fase'])]
if freq_iguales:
        x_amplitud.append(float(x_suma_args['amplitud']))
        x_fase.append(float(x_suma_args['fase']))
        x_etiqueta.append('x_suma')

# graficar
fig_polar  = plt.figure()
graf_polar = fig_polar.add_subplot(projection='polar')
punto = graf_polar.plot(x_fase,x_amplitud,'o',xunits='radians')
for i in range(0,len(x_amplitud),1):
    if (i<=9): # color de la línea
        lineacolor = 'C'+str(i)
    else:
        numcolor = i%10
        lineacolor = 'C'+str(numcolor)
    graf_polar.quiver(0,0,x_fase[i],x_amplitud[i],
                      angles='xy', scale_units='xy', scale=1,
                      label = x_etiqueta[i],
                      color = lineacolor)
# actualiza etiquetas de gráfico polar
etiquetas = ['0', '$\pi$/4', '$\pi$/2','3$\pi$/4',
             '$\pi$','5$\pi$/4','3$\pi$/2','7$\pi$/4']
dtheta = int(360/len(etiquetas))
lines, labels = plt.thetagrids(range(0, 360,dtheta),
                               etiquetas)
graf_polar.set_title('Fasores Xk a '+str(x1_args['freq_Hz'])+ ' Hz')
graf_polar.legend()
plt.show()

[ fasores ] [ ejemplo ] [ analítico ] [ algoritmo ] [ gráfica] [ función ]
..


6. Función

6.1 Incorporar parámetro fasor a dsp.cos_args_one_term()

El fasor puede considerarse como parte de los parámetros de una señal, por lo que el cálculo se incorpora a la instrucción dsp.cos_args_one_term() como parámetro: 'fasor' y 'fasor_complex'.

La disponibilidad del parámetro 'fasor' puede facilitar las operaciones de suma de señales de igual frecuencia cuando se usen mas componentes de la señal. La simplificación del algoritmo se muestra en el siguiente ejercicio:

6.2 Ejercicio

Referencia: McClellan ejercicio 2.8 p54

Considere las dos señales:

x_1 (t) = 5\cos(2 \pi (100) t + \pi /3) x_2 (t) = 4\cos(2 \pi(100) t + \pi / 4)

6.3 Algoritmo en Python

Considerando que se pueden usar mas señales de igual frecuencia para sumar los fasores, se procede a usar x_senales como una lista de señales a ser analizada, junto a x_etiqueta con las etiquetas para las gráficas.

# INGRESO
x1 = 5*sym.cos(DosPi*100*t + pi/3)
x2 = 4*sym.cos(DosPi*100*t - pi/4)

x_senales = [x1,x2] # sumar fasores de igual frecuencia
x_etiqueta = ['x1(t)','x2(t)']

La revisión de las frecuencias de cada señal se realiza con el diccionario x_freq que tendra el 'conteo' de señales con la misma frecuencia y 'fasor_complex' que se acumulará cada vez que se encuentre una señal de igual frecuencia en la lista x_senales.

Luego se realiza una revisión de las frecuencias que tienen conteo>1 para generar la senal x_suma y añadir a la lista x_senales con la etiqueta correspondiente en x_etiqueta. Se debe actualizar el x_conteo que es usado para enumerar las gráficas.

La salida del algoritmo para el ejercicio se presenta como:

freq_Hz: 100
  conteo : 2
  fasor_complex : 5.32842712474619 + 1.501699894176*I
  senal : 5.53599477925144*cos(100*t*(2*pi) + 0.27470298976466)

Las simplificación del algoritmo usando la función actualizada y las señales en una lista se muestran a continuación:

Instrucciones en Python.

# ejercicio 2.8 p54 Suma de Fasores y graf_polar
# telg1034 DSP fiec-espol edelros@espol.edu.ec
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
import telg1034 as dsp

# variables continuas
from telg1034 import t,A,w,f,p,pi,DosPi,I,equivalentes

# INGRESO
x1 = 5*sym.cos(DosPi*100*t + pi/3)
x2 = 4*sym.cos(DosPi*100*t - pi/4)

x_senales = [x1,x2] # sumar fasores de igual frecuencia
x_etiqueta = ['x1(t)','x2(t)']

# PROCEDIMIENTO
x_conteo = len(x_senales)
# sumar fasores por frecuencias iguales
x_freq = {}
for unasenal in x_senales:
    xk_args = dsp.cos_args_one_term(unasenal)    
    freq_k = xk_args['freq_Hz']
    fasor_complex_k = xk_args['fasor_complex']
    if freq_k in x_freq: # sumar fasores en la freq_k
        x_freq[freq_k]['conteo'] = x_freq[freq_k]['conteo'] + 1
        fasor_complex = x_freq[freq_k]['fasor_complex']
        x_freq[freq_k]['fasor_complex'] = fasor_complex + fasor_complex_k
    else: # primer fasor en freq_k
        x_freq[freq_k] = {}
        x_freq[freq_k]['conteo'] = 1
        x_freq[freq_k]['fasor_complex'] = fasor_complex_k

# señal x_suma(t) por frecuencias 
for unafreq in x_freq:
    if x_freq[unafreq]['conteo']>1: # existe x_suma
        fasor_complex = x_freq[unafreq]['fasor_complex']
        amplitud = sym.Abs(fasor_complex)
        fase = sym.arg(fasor_complex)
        fasor_suma = amplitud*sym.exp(sym.UnevaluatedExpr(fase))
        x_suma = amplitud*sym.cos(DosPi*unafreq*t+sym.UnevaluatedExpr(fase))
        x_freq[unafreq]['senal'] = x_suma
        # actualiza lista de señales
        x_senales.append(x_suma)
        x_etiqueta.append('x_'+str(unafreq)+'Hz')
        x_conteo = len(x_senales)

# SALIDA
for unafreq in x_freq:
    print('freq_Hz:',unafreq)
    for unparametro in x_freq[unafreq]:
        print(' ',unparametro,':', x_freq[unafreq][unparametro])

6.4 funciones para gráficas

En el desarrollo de las gráficas una operación recurrente es encontrar el periodo máximo T_max, por lo que periodo_minmax(x_freq) que se incorpora a telg1034.py.

De forma semejante se requiere un color diferente por cada componente de señal en x_senales dado por su índice en la lista con la función auxiliar . De esa forma se facilita la identificación entre la gráfica de fasores y la gráfica de tiempo _graf_lineacolor_i(i).

Para las gráficas de la suma de fasores y señales en tiempo se requiere realizar las listas de parámetros para cada gráfica.

Se destaca la relación entre las gráfica de fasores y en tiempo al rotar graf_polar como se muestra en la imagen.

grafica suma fasores 02

Instrucciones adicionales para las gráficas

# GRAFICA de suma de fasores ---------
# PROCEDIMIENTO de suma de fasores
x_amplitud = [] ; x_fase = [] ; x_ceros = []
for unasenal in x_senales:
    parametro = dsp.cos_args_one_term(unasenal)
    freq_k = parametro['freq_Hz']
    ampl_k = parametro['amplitud']
    fase_k = parametro['fase']
    x_amplitud.append(float(ampl_k))
    x_fase.append(float(fase_k))
    # marca un cero cos(wt+p)=0; wt+p=pi/2
    # t_cero = (pi/2-p)/w = (pi/2-p)/(2*pi*f)
    uncero = sym.nan # si es una constante
    if abs(freq_k)>0:
        uncero = (pi/2 - fase_k)/(2*pi*freq_k)
        uncero = np.around(float(uncero),decimals=15)
    x_ceros.append(uncero)

# graficar
def _graf_lineacolor_i(i):
    ''' función auxiliar, asigna el color en la paleta de matplotlib
    acorde al número i entre 0 y 9
    '''
    if (i<=9): # color de la línea
        lineacolor = 'C'+str(i)
    else:
        numcolor = i%10
        lineacolor = 'C'+str(numcolor)
    return(lineacolor)

fig_polar  = plt.figure()
graf_polar = fig_polar.add_subplot(projection='polar')
punto = graf_polar.plot(x_fase,x_amplitud,'o',xunits='radians')
for i in range(0,len(x_amplitud),1):
    graf_polar.quiver(0,0,x_fase[i],x_amplitud[i],
                      angles='xy', scale_units='xy', scale=1,
                      label = x_etiqueta[i],
                      color = _graf_lineacolor_i(i))
# actualiza etiquetas de gráfico polar
etiquetas = ['0', '$\pi$/4', '$\pi$/2','3$\pi$/4',
             '$\pi$','5$\pi$/4','3$\pi$/2','7$\pi$/4']
dtheta = int(360/len(etiquetas))
lines, labels = plt.thetagrids(range(0, 360,dtheta),
                               etiquetas)
graf_polar.set_title('Fasores Xk')
graf_polar.legend()
#plt.show()


# GRAFICAR x(t) ---------------------------
# INGRESO
tramos_T = 20       # tramos por periodo
periodos_graf = 2   # periodos en una grafica

# PROCEDIMIENTO GRAFICA
# periodo máximo y mínimo de señales
def periodo_minmax(x_freq):
    ''' función auxiliar: obtiene T_max y T_min
    de una lista de frecuencias x_freq
    ingresa: x_freq
    salida:  [T_max,Tmin].
    Si la frecuencia es cero, T_max=1,T_min=1
    '''
    T_max = 1 ; T_min = 1
    freq = np.array(list(x_freq))
    freq_min = float(np.min(np.abs(freq)))
    freq_max = float(np.max(np.abs(freq)))
    d_freq = np.abs(freq_max - freq_min) # ancho de banda
    if freq_max>0:
        T_min = 1/freq_max
    if d_freq>0: # min y max no son iguales
        if freq_min>0:
            T_max = 1/d_freq
        else:  # freq_min==0, selecionar otro T_max
            freq_min = float(np.min(freq[freq>0]))
            T_max = 1/freq_min 
    else: # min y max son iguales
        T_max = 1/freq_min
    T_minmax = [T_min,T_max] 
    return(T_minmax)

[T_min,T_max] = periodo_minmax(x_freq)

muestras_T = tramos_T + 1
# intervalo de t entre [a,b] en pasos dt
a = 0
b = T_max*periodos_graf
dt = T_min/tramos_T # tamaño de paso,tramo, muestreo
ti = np.arange(a,b+dt,dt)

xi_k = [] # muestras de cada señal x
for unasenal in x_senales:
    # x(t) a forma numérica lambda t:
    if unasenal.has(t):
        xk = sym.lambdify(t,unasenal)
    else: # es constante
        xk = lambda t: unasenal + 0*t
    xki = xk(ti) # evalúa en intervalo
    xi_k.append(np.copy(xki))

# graficar
fig_x = plt.figure()
graf_x = fig_x.add_subplot()
for i in range(0,x_conteo,1):
    color_i = _graf_lineacolor_i(i)
    graf_x.plot(ti,xi_k[i], color=color_i,
                label=x_etiqueta[i])
    graf_x.axhline(0,color='black')
    if x_ceros[i]!=sym.nan:
        graf_x.axvline(x_ceros[i],color=color_i,
                       linestyle='dashed')
graf_x.set_title('Señales xk(t)')
graf_x.set_xlabel('t (segundos)')
graf_x.set_ylabel('x_senales')
graf_x.grid()
graf_x.legend()
plt.show()

[ fasores ] [ ejemplo ] [ analítico ] [ algoritmo ] [ gráfica] [ función ]

1.2 Señal – Periodo y Desplazamiento en tiempo con Sympy-Python

[ periodo ] [ desplazamiento ] [ ejercicio ] [ algoritmo ] [ gráfica ] [_telg1034.py _as_dsp ]
..


1. Periodo de x(t)

Referencia: McClellan 2.3.1 p35, 2.3.2 p37

El periodo T0 de una sinusoidal es el tiempo que dura un ciclo de la señal. La frecuencia en Hz determina su periodo, es decir que al aplicar  un desplazamiento de T0, se repite la misma señal.

T_0 = \frac{2\pi}{\omega_0} = \frac{1}{f_0} x(t+T_0) = x(t)

..

2. Desplazamiento en tiempo de x(t)

El efecto de aplicar a t un desplazamiento t0 diferente al periodo T0, presenta un desfase de adelanto o retraso según el signo de t0 como se muestra en la gráfica. Siendo x(t) = A cos(ωt)

x(t+t_0) = A cos(\omega_0(t-t_0)) = A cos(\omega_0 t-\omega_0 t_0) = A cos(\omega_0 t + \varphi) \varphi= -\omega_0 t_0

Desde el punto de vista de la fase de la señal, se puede observar que:

\varphi= -\omega_0 t_0 = -2 pi \Big(\frac{t_0}{T_0} \Big)

El tema también se trata en el curso de Señales y Sistemas, donde se detalla:

Señales con desplazamiento. escalamiento o inversión en tiempo

3. Ceros de x(t)

Para disponer de una marca referencial de la función, se usará el primer cero de x(t), por simplicidad e obtener mediante:

\cos(\omega t + \varphi) = 0

donde se cumple que

\omega t + \varphi = \frac{\pi}{2}

que permtite despejar el tiempo cuando ocurre el cero en un periodo:

\omega t = \frac{\pi}{2} - \varphi t_{cero} = \frac{\frac{\pi}{2} -\varphi}{\omega}=\frac{\frac{\pi}{2} - \varphi}{2\pi f_0}

[ periodo ] [ desplazamiento ] [ ejercicio ] [ algoritmo ] [ gráfica ] [_telg1034.py _as_dsp ]
..


4. Ejercicio

Referencia: McClellan Fig 2.7 p36

Realizar las gráficas para x(t) con varios valores de f0 = [ 200, 100, 0 ] y observar las diferencias entre ellas.

x(t) =5 \cos(2 \pi f_0 t )

Añadir una señal que se desplace en la tercera parte de su periodo (T/3).

Listar los parámetros de cada señal.

[ periodo ] [ desplazamiento ] [ ejercicio ] [ algoritmo ] [ gráfica ] [_telg1034.py _as_dsp ]
..


5. Algoritmo con Sympy

El ejercicio planteado requiere analizar más de una señal, por lo que es conveniente agruparlas en una lista x_senales, para luego usar una funciones como  cos_args_one_term() que obtenga los parámetros de cada señal. La función contiene el algoritmo realizado en Señales – Parámetros de sinusoides con Sympy-Python.

Los parámetros se agrupan en un diccionario, identificando cada uno por el nombre en la entrada. Se incorpora en los parámetros el valor de 'periodo'

Por simplicidad, el análisis del argumento del coseno que es una expresión más simple, se realiza en la función auxiliar _cos_arg_freqfase() con instrucciones para simplificar las expresiones cuando se aplican desplazamientos en el tiempo, según se ha mostrado en la parte teórica.

Con los parámetros de cada señal se agrupan las frecuencias, amplitud y fase, añadiendo el primer cero de x(t) como una referencia en el primer periodo y destacar los desplazamientos de una señal. (ejemplo x3(t))

Resultados del algoritmo

señales revisadas:  4
amplitudes:  [5, 5, 5, 5]
frecuencias: [200, 100, 100, 0]
fases:       [0, 0, -0.333333333333333*pi, 0]
periodo máximo: 0.02
periodo mínimo: 0.005
raíces x_ceros: [0.00125, 0.0025, 0.004166666666667, nan]
>>>

Para realizar las gráficas, se unificará el eje de tiempo usando la frecuencia más lenta o periodo máximo. El procedimiento es semejante al ejercicio anterior,  añadiendo sub-gráficas (subplots) para cada una de las señales.

gráfica de varias señales x(t)

Las funciones presentadas se prueban y se mejoran en cada ejercicio o tema de la unidad, por ejemplo aún no se analizan señales con varios término suma. Las funciones se añaden a un archivo telg1034.py para ser reutilizadas y actualizadas.

Instrucciones en Python

# Figura 2.3.1 p35 parámetros en diccionario de x(t) 
# telg1034 DSP fiec-espol edelros@espol.edu.ec
import sympy as sym
import numpy as np
import matplotlib.pyplot as plt

# variables continuas
t = sym.Symbol('t', real=True,positive=True)
pi = sym.pi
DosPi = sym.UnevaluatedExpr(2*sym.pi)

# INGRESO
x1 = 5*sym.cos(DosPi*(200)*t)
x2 = 5*sym.cos(DosPi*(100)*t)
x3 = x2.subs(t,t-0.005/3) # adelanto en x(t-t0)
x4 = 5*sym.cos(DosPi*(0)*t)

x_senales = [x1,x2,x3,x4]
x_etiqueta = ['x1(t)','x2(t)','x3(t)','x4(t)']

# PROCEDIMIENTO   
def cos_args_one_term(x):
    ''' versión 1. Amplitud, frecuencia en Hz y rad, fase de sin() o cos()
    referenciado a un término coseno.
    Entrega diccionario con 'amplitud', 'freq_Hz','freq_rad','fase','periodo'.
    '''
    parametro={} # parámetros en diccionario
    amplitud = sym.S.One ; T = sym.S.Zero ; fase = sym.S.Zero
    freq_Hz = sym.S.Zero ; freq_rad = sym.S.Zero

    x = sym.sympify(x) # expresión a sympy, por si es constante
    
    SinSumas = True # Analizar un solo cos o sin
    if x.is_Add:
        SinSumas = False

    if SinSumas: # x(t) un solo cos() o sin():
        partes = x.args
        for unaparte in partes:
            cond1 = unaparte.has(sym.cos)
            cond2 = unaparte.has(sym.sin)
            cond3 = unaparte.has(t)
            # revisa argumento_cos(w*t+p)
            if (cond1 or cond2) and cond3:
                argumento_cos = unaparte.args[0]
                [freq_Hz,freq_rad,fase] = _cos_arg_freqfase(argumento_cos)
                if unaparte.has(sym.sin): # convierte a cos(), estandariza
                    fase = fase - sym.pi/2
            else:
                amplitud = amplitud*unaparte
                
        # amplitud, revisión numérica
        if amplitud.is_Number:
            if amplitud<0: # amplitud con signo positivo
                amplitud = np.abs(amplitud)
                fase = fase + sym.pi
            if isinstance(amplitud,sym.Float):
                amplitud = float(amplitud)
            if isinstance(amplitud,sym.Integer):
                amplitud = int(amplitud)
            # mantiene forma racional (a/b) de amplitud
        # periodo T
        if freq_Hz != sym.S.Zero:
            T = 1/freq_Hz
        else: # es una constante
            T = sym.S.NaN
    if not(x.has(t)): # es una constante
        amplitud = x
    parametro['amplitud'] = amplitud
    parametro['freq_Hz']  = freq_Hz
    parametro['freq_rad'] = freq_rad
    parametro['fase'] = fase
    parametro['periodo'] = T
                
    if SinSumas == False:
        msg ='x(t) tiene suma de terminos, usar solo un término de la forma:'
        msg = msg + '\n A*cos(w*t+p) o también A*sin(w*t+p) ... \n'
        print('cos_sin_args:',x,msg)
    return(parametro)

def _cos_arg_freqfase(argumento_cos):
    ''' _función_auxiliar_ que obtiene frecuencia (Hz y rad) y fase de la
    expresion: 2*pi*freq_Hz*t + fase
    entrega: [freq_Hz,freq_rad,fase]
    '''
    freq_Hz = sym.S.Zero ; freq_rad = sym.S.Zero ; fase = sym.S.Zero
    argumento_cos = sym.sympify(argumento_cos) # si es constante
    # simplifica  w*t + fase
    argumento_cos = argumento_cos.doit() # evalua sym.UnevaluatedExpr()
    argumento_cos = sym.expand(argumento_cos) # simplifica polinomio t
    argumento_cos = sym.collect(argumento_cos,t) # agrupa terminos t
    
    if argumento_cos.is_Add: # términos suma: w*t + fase
        term_suma = sym.Add.make_args(argumento_cos)
        for term_k in term_suma:
            if term_k.has(t): # w*t
                freq_rad = sym.simplify(term_k/t)
                freq_Hz = sym.simplify(freq_rad/(2*pi))
            else: # fase suma término sin t
                fase = fase + term_k
    if argumento_cos.is_Mul: # termino único: w*t
        if argumento_cos.has(t):
            freq_rad = sym.simplify(argumento_cos/t)
            freq_Hz = sym.simplify(freq_rad/(2*pi))
        else:
            fase = fase + argumento_cos
    parametros_arg = [freq_Hz,freq_rad,fase]
    return(parametros_arg)

# DESARROLLO -----------------------
x_conteo = len(x_senales)
# lista frecuencias en señales
x_freq = [];  x_amplitud =[]; x_fase = []; x_ceros = []
for unasenal in x_senales:
    parametro = cos_args_one_term(unasenal)
    freq_k = parametro['freq_Hz']
    ampl_k = parametro['amplitud']
    fase_k = parametro['fase']
    x_freq.append(freq_k)
    x_amplitud.append(ampl_k)
    x_fase.append(fase_k)
    # marca un cero cos(wt+p)=0; wt+p=pi/2
    # t_cero = (pi/2-p)/w = (pi/2-p)/(2*pi*f)
    uncero = sym.nan # si es una constante
    if abs(freq_k)>0:
        uncero = (pi/2 - fase_k)/(2*pi*freq_k)
        uncero = np.around(float(uncero),decimals=15)
    x_ceros.append(uncero)

# periodo máximo y mínimo de señales
T_max = 1 ; T_min = 1
freq = np.array(x_freq)
freq_min = float(np.min(np.abs(freq)))
freq_max = float(np.max(np.abs(freq)))
d_freq = np.abs(freq_max - freq_min) # ancho de banda
if freq_max>0:
    T_min = 1/freq_max
if d_freq>0: # min y max no son iguales 
    if freq_min>0:
        T_max = 1/d_freq
    else:  # freq_min==0, selecionar otro T_max
        freq_min = float(np.min(freq[freq>0]))
        T_max = 1/freq_min 
else: # min y max son iguales
    T_max = 1/freq_min

# SALIDA
print('señales revisadas: ', x_conteo)
print('amplitudes: ', x_amplitud)
print('frecuencias:', x_freq)
print('fases:      ', x_fase)
print('periodo máximo:', T_max)
print('periodo mínimo:', T_min)
print('raíces x_ceros:', x_ceros)

[ periodo ] [ desplazamiento ] [ ejercicio ] [ algoritmo ] [ gráfica ] [_telg1034.py _as_dsp ]
..


6. Gráfica con subplots

Para la parte gráfica se amplían las instrucciones anteriores, añadiendo sub-gráficas (subplots) para cada una de las señales.

Referencia: matplotlib.pyplot.subplots

Para esta ocasión se aborda el desplazamiento considerando que la expresión usada es con Sympy. Para la gráfica que requiere la conversión de forma simbólica a numérica con la instrucción lambdify para más de una señal. como se muestra en el ejercicio.

Conversión de forma simbólica a forma numérica Lambda ( curso: métodos numéricos)

Para la gráfica se pueden usar el parámetro del periodo obtenido en el algoritmo. por lo que será necesario solamente añadir la cantidad de periodos a observar de la señal con menor frecuencia. Para mostrar una curva «suave» se indicará la cantidad de tramos por periodo que se usará para las muestras, se sugiere para el coseno al menos 20.

# GRAFICAR x(t) ---------------------------
# INGRESO
tramos_T = 20       # tramos por periodo
periodos_graf = 4   # periodos en una grafica

gráfica de varias señales x(t)Las señales contenidas en la lista x_senales requieren muestras de cada una, que se agrupan en xi_k, reutilizando el proceso usando en la gráfica realizada al inicio de la unidad.

El primer cero de la función se marca con una línea vertical discontínua (dashed) de color rojo en los puntos encontrados en x_ceros. Cuando una señal es una constante, el cero se marca con sym.nan (not a number) que no genera la línea roja.

    if x_ceros[i]!=sym.nan:
        graf_x.axvline(x_ceros[i],color='red',
                       linestyle='dashed')

Instrucciones para la gráfica

# GRAFICAR x(t) ---------------------------
# INGRESO
tramos_T = 20       # tramos por periodo
periodos_graf = 4   # periodos en una grafica

# PROCEDIMIENTO GRAFICA
muestras_T = tramos_T + 1
# intervalo de t entre [a,b] en pasos dt
a = 0
b = T_max*periodos_graf
dt = T_min/tramos_T # tamaño de paso,tramo, muestreo
ti = np.arange(a,b+dt,dt)

xi_k = [] # muestras de cada señal x
for unasenal in x_senales:
    # x(t) a forma numérica lambda t:
    if unasenal.has(t):
        xk = sym.lambdify(t,unasenal)
    else: # es constante
        xk = lambda t: unasenal + 0*t
    xki = xk(ti) # evalúa en intervalo
    xi_k.append(np.copy(xki))

# graficar
fig_x = plt.figure()
for i in range(0,x_conteo,1):
    graf_sub = x_conteo*100+10+i+1
    graf_x = fig_x.add_subplot(graf_sub)
    graf_x.plot(ti,xi_k[i])
    graf_x.set_ylabel(x_etiqueta[i])
    graf_x.grid()
    graf_x.axhline(0,color='black')
    if i == 0:
        graf_x.set_title('Señales xk(t)')
    if x_ceros[i]!=sym.nan:
        graf_x.axvline(x_ceros[i],color='red',
                       linestyle='dashed')
graf_x.set_xlabel('t (segundos)')
# graf_x.legend()
plt.show()

[ periodo ] [ desplazamiento ] [ ejercicio ] [ algoritmo ] [ gráfica ] [_telg1034.py _as_dsp ]
..


7. Funciones de telg1034.py como dsp

Las funciones para uso en otras unidades se agrupan en el archivo telg1034.py que permite simplificar los algoritmos como el del ejercicio anterior. Las funciones se llaman usando el alias dsp

import telg1034 as dsp

La declaración de variables contínuas como símbolo también se simplifican a una línea. Se incorpora la forma equivalentes de funciones Sympy y Numpy al usar lambdify.

# variables continuas
from telg1034 import t,A,w,f,p,pi,DosPi,I,equivalentes

En los siguientes ejercicios se hará uso de telg1034.py para resumir las funciones de cada algoritmo que se incorpore en la unidad.

7.1 Ejemplo de desplazamiento por adelanto o retraso

Referencia: McClellan ejercicio 2.4 p39

Usando x(t), realice la gráfica para x(t-t1) cuando t1 = 0.0075. Luego repita el proceso para t1 = -0.01. Revise la dirección del desplazamiento. Para cada caso calcule la fase de la sinusoide desplazada en el tiempo.

x_1(t) = 20 cos(2\pi(40)t - 0.4\pi) x_2(t) = x_1( t-0.0075) x_3(t) = x_1( t-(-0.01))

Las señales indicadas en el ejercico se ingresan al algoritmo en la lista x_senales. Los desplazamientos se realizan al sustituir t por (t-desplaza).

# INGRESO
x1 = 20*sym.cos(DosPi*40*t-0.4*pi)
x2 = x1.subs(t,t-0.0075)
x3 = x1.subs(t,t-(-0.01))

con los siguientes resultados:

señales revisadas:  3
amplitudes:  [20, 20, 20]
frecuencias: [40, 40, 40]
fases:       [-0.4*pi, -1.0*pi, 0.4*pi]
periodo máximo: 0.05
periodo mínimo: 0.025
raices x_ceros: [0.01125, 0.01875, 0.00125]
>>>

Para la parte gráfica se usa las instrucciones dadas en la sección [ gráfica ]:

señal desplazada en el tiempo

Instrucciones en Python

Instrucciones simplificadas, sin la sección gráfica.

# ejemplo 2.4 p39 x(t) desplazamiento adelanto o retraso
# telg1034 DSP fiec-espol edelros@espol.edu.ec
import sympy as sym
import numpy as np
import matplotlib.pyplot as plt
import telg1034 as dsp

# variables continuas
from telg1034 import t,A,w,f,p,pi,DosPi,I,equivalentes

# INGRESO
x1 = 20*sym.cos(DosPi*40*t-0.4*pi)
x2 = x1.subs(t,t-0.0075)
x3 = x1.subs(t,t-(-0.01))

x_senales = [x1,x2,x3]
x_etiqueta = ['x(t)','x(t-0.0075)','x(t+0.01)']

# PROCEDIMIENTO   
x_conteo = len(x_senales)
# lista frecuencias en señales
x_freq = [];  x_amplitud =[]; x_fase = []; x_ceros = []
for unasenal in x_senales:
    parametro = dsp.cos_args_one_term(unasenal)
    freq_k = parametro['freq_Hz']
    ampl_k = parametro['amplitud']
    fase_k = parametro['fase']
    x_freq.append(freq_k)
    x_amplitud.append(ampl_k)
    x_fase.append(fase_k)
    # marca un cero cos(wt+p)=0; wt+p=pi/2
    # t_cero = (pi/2-p)/w = (pi/2-p)/(2*pi*f)
    uncero = sym.nan # si es una constante
    if abs(freq_k)>0:
        uncero = (pi/2 - fase_k)/(2*pi*freq_k)
        uncero = np.around(float(uncero),decimals=15)
    x_ceros.append(uncero)

# periodo máximo y mínimo de señales
T_max = 1 ; T_min = 1
freq = np.array(x_freq)
freq_min = float(np.min(np.abs(freq)))
freq_max = float(np.max(np.abs(freq)))
d_freq = np.abs(freq_max - freq_min) # ancho de banda
if freq_max>0:
    T_min = 1/freq_max
if d_freq>0: # min y max no son iguales 
    if freq_min>0:
        T_max = 1/d_freq
    else:  # freq_min==0, selecionar otro T_max
        freq_min = float(np.min(freq[freq>0]))
        T_max = 1/freq_min 
else: # min y max son iguales
    T_max = 1/freq_min

# SALIDA
print('señales revisadas: ', x_conteo)
print('amplitudes: ', x_amplitud)
print('frecuencias:', x_freq)
print('fases:      ', x_fase)
print('periodo máximo:', T_max)
print('periodo mínimo:', T_min)
print('raíces x_ceros:', x_ceros)

[ periodo ] [ desplazamiento ] [ ejercicio ] [ algoritmo ] [ gráfica ] [_telg1034.py _as_dsp ]

1.1 Señal – Parámetros de sinusoides con Sympy-Python

Señal [ sinusoidal ] [ ejercicio ] [ analítico ] [ algoritmo ] [ gráfica ]
..


1. Señales Sinusoidales

La forma general de una señal sinusoidal en el dominio del tiempo se obtiene por medio de su argumento de una función coseno de variable independiente t.

x(t) = A\cos(\omega_0 t +\phi) = A\cos(2\pi f_0 t +\phi)

por lo que se deduce que:

\omega_0 = 2\pi f_0

se desatacan tres parámetros independientes: Amplitud (A), fase (φ), frecuencia angular (ω) en radianes o la frecuencia por ciclos (f) en Hz.

Señal [ sinusoidal ] [ ejercicio ] [ analítico ] [ algoritmo ] [ gráfica ]
..


2. Ejercicio

Referencia: McClellan ejemplo 2.1 p35, Downey 1.1 p2.

Para la señal definida por x(t) mostrar los parámetros de Amplitud, Frecuencia y Fase.

x(t) = 20 \cos(2\pi(40)t - 0.4\pi)

señal sinusoidal

Señal [ sinusoidal ] [ ejercicio ] [ analítico ] [ algoritmo ] [ gráfica ]
..


3. Desarrollo analítico

Por observación de x(t) según los conceptos definiciones en la parte teórica:

x(t) = A\cos(\omega_0 t +\phi)

los parámetros buscados son:

Amplitud A = 20
ω0 = 2 π(40) rad/s
f0 = 40 Hz
φ = -0.4π rad.

El tamaño vertical de la señal en el tiempo depende del parámetro Amplitud, su máximo y mínimo +20 a -20.

El tiempo de entre máximos sucesivos es de 0.025 s, que es 1/f0 .

Señal [ sinusoidal ] [ ejercicio ] [ analítico ] [ algoritmo ] [ gráfica ]
..


4. Algoritmo con Sympy

Cuando para x(t) se requiere obtener sus parámetros o características, se usa la forma simbólica con las herramientas de la librería Sympy.

Luego de hacer el llamado a Sympy,  procede a  indicar la variable independiente de tiempo 't' para ser reconocida como un símbolo. También se pueden añadir constantes como π o 2π para facilitar la escritura de la expresión x(t):

# ejemplo 2.1 p35 cos_args grafica senoidal
# telg1034 DSP fiec-espol edelros@espol.edu
import sympy as sym

# variables continuas
t = sym.Symbol('t', real=True,positive=True)
pi = sym.pi
DosPi = sym.UnevaluatedExpr(2*sym.pi)

La expresión x(t) se escribe usando los símbolos definidos en el ejercicio.

# INGRESO
x = 20*sym.cos(DosPi*40*t-0.4*pi)

Como la expresión x(t) es simple y tiene un patrón Amplitud*cos(Freq_rad*t+Fase), se separan sus partes con la instrucción x.args. La separación de argumentos se realiza con expresiones de sumas de términos, factores multiplicadores o argumentos de una función como el cos().

# PROCEDIMIENTO
amplitud = 1; freq_rad = 0; fase = 0 
partes = x.args

El resultado de las partes es una tupla. La Amplitud corresponde a la parte que no tiene sym.cos, el resto debería ser el la parte del coseno. Si la expresión contiene mas de dos partes, se analiza la expresión iterando en cada parte.

parte_cos = sym.S.Zero # busca cos(t)
for unaparte in partes:
    if not(unaparte.has(sym.cos)):
        amplitud = partes[0]
    else:
        parte_cos = unaparte

La expresión interna de cos(Freq_rad*t+Fase) se obtiene también como el  argumento tupla de un solo elemento, obtenida como parte_cos.args[0].

# revisa cos(freq_rad*t+fase)
argumento_cos = parte_cos.args[0]

Los términos de la suma se revisan uno por uno, La frecuencia angular se obtiene como el término que tiene símbolo 't', que al dividir para 't' y se obtiene el valor en en radianes. La revisión contenido de t en un término  se realiza por medio untermino.has(t),  La fase es el término sin el símbolo 't'.

terminos_suma = argumento_cos.args
for untermino in terminos_suma:
    if untermino.has(t):
        freq_rad = sym.simplify(untermino/t)
        freq_Hz = sym.simplify(freq_rad/(2*pi))
    else:
        fase = untermino

La instrucción sym.simplify() permite simplificar los valores como π al calcular freq_Hz

El bloque de resultados muestra en orden los valores obtenidos durante el proceso

# SALIDA
print('x(t): ')
sym.pprint(x)
print('\npartes o argumentos')
print('partes : ', partes)
print('parte_cos: ',parte_cos)
print('argumento_cos: ',argumento_cos)
print('términos_suma: ', terminos_suma)

print('\nparámetros')
print('amplitud: ',amplitud)
print('freq_rad: ',freq_rad)
print('freq_Hz : ',freq_Hz)
print('fase : ',fase)

los resultados se presentan a continuación:

x(t): 
20*cos(40*t*2*pi - 0.4*pi)

partes o argumentos
partes   :  (20, cos(40*t*(2*pi) - 0.4*pi))
parte_cos:  cos(40*t*(2*pi) - 0.4*pi)
argumento_cos:  40*t*(2*pi) - 0.4*pi
términos_suma:  (-0.4*pi, 40*t*(2*pi))

parámetros
amplitud:  20
freq_rad:  80*pi
freq_Hz :  40
fase    :  -0.4*pi
>>> 

Instrucciones en Python

# ejemplo 2.1 p35 parámetros de x(t)
# telg1034 DSP fiec-espol edelros@espol.edu.ec
import sympy as sym

# variables continuas
t = sym.Symbol('t', real=True,positive=True)
pi = sym.pi
DosPi = sym.UnevaluatedExpr(2*sym.pi)

# INGRESO
x = 20*sym.cos(DosPi*40*t-0.4*pi)

# PROCEDIMIENTO
amplitud = 1; freq_rad = 0; fase = 0 
partes = x.args
parte_cos = sym.S.Zero # busca cos(t)
for unaparte in partes:
    if not(unaparte.has(sym.cos)):
        amplitud = partes[0]
    else:
        parte_cos = unaparte

# revisa cos(freq_rad*t+fase)
argumento_cos = parte_cos.args[0]
terminos_suma = argumento_cos.args
for untermino in terminos_suma:
    if untermino.has(t):
        freq_rad = sym.simplify(untermino/t)
        freq_Hz = sym.simplify(freq_rad/(2*pi))
    else:
        fase = untermino

# SALIDA
print('x(t): ')
sym.pprint(x)
print('\npartes o argumentos')
print('partes   : ', partes)
print('parte_cos: ',parte_cos)
print('argumento_cos: ',argumento_cos)
print('términos_suma: ', terminos_suma)

print('\nparámetros')
print('amplitud: ',amplitud)
print('freq_rad: ',freq_rad)
print('freq_Hz : ',freq_Hz)
print('fase    : ',fase)

[ ejercicio ] [ analítico ] [ algoritmo ] [ tarea ]
..


5. Gráfica de la señal sinusoidal

Para realizar la gráfica de la función se debe añadir las librerías para los cálculos (Numpy) y la librería para realizar la gráfica (Matplotlib)

# GRAFICAR x(t) ---------------------------
import numpy as np
import matplotlib.pyplot as plt

Como parámetros de la función periódica se puede indicar el número de periodos a mostrar en la gráfica. Recordando que el periodo T es el inverso de la frecuencia f en Hz.

# INGRESO
muestreo_T = 20+1 # tramos+1, o pasos+1
graf_periodos = 2 # periodos en una grafica

Como la gráfica básica realiza un muestreo de la señal y luego une los puntos con rectas, para observar una curva «suave» se usaría un muestreo de unos 20 tramos o 21 muestras por periodo.

señal sinusoidal

Las muestras de tiempo ti se obtienen con los datos de [a,b] y dt usando la instrucción np.arange().

Para la evaluación de x(t) en el intervalo con el vector ti es preferible usar la forma numérica (lambda) que usa todo el vector ti en una sola instrucción. La forma de evaluación de Sympy es con xt.subs() y sería con valores individuales. Se realiza la conversión de Sympy a Numpy usando la instrucción sym.lambdify(variable, expresión). (lambfify en resumen de sympy)

# x(t) a forma numérica lambda t:
xt = sym.lambdify(t,x)
xi = xt(ti) # evalúa en intervalo

Para revisar o recordar los detalles de la creación de la gráfica con plt.plot(), puede revisar lo indicado en los enlaces siguientes:

Señales Contínuas con Python (curso: Señales y Sistemas)

Gráficas 2D de línea (curso: Fundamentos de Programación)

Observación:  Para la gráfica se requieren valores numéricos reales (float) en los vectores [ti,xi]. Si la expresión de x(t) aún no se definen los valores de los parámetros, la sección de la gráfica mostraría un error. Considerar para el próximo algoritmo.

Instrucciones en Python

# GRAFICAR x(t) ---------------------------
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
muestreo_T = 20+1 # tramos+1, o pasos+1
periodos_graf = 2 # periodos en una grafica

# PROCEDIMIENTO GRAFICA
T = 1/freq_Hz

a = 0    # [a,b] intervalo ti
b = T*periodos_graf
dt = T/(muestreo_T-1)  # tamaño de paso,tramo, muestreo
ti = np.arange(a,b+dt,dt,dtype=float)

# x(t) a forma numérica lambda t:
xt = sym.lambdify(t,x)
xi = xt(ti) # evalúa en intervalo

# graficar
plt.plot(ti,xi,label='x(t)')
plt.plot(ti,xi,'o',label='x(ti)')
plt.axhline(0,color='black')
plt.xlabel('t (segundos)')
plt.ylabel('y')
plt.title('x(t) = ' + str(x))
plt.legend()
plt.grid()
plt.show()

Señal [ sinusoidal ] [ ejercicio ] [ analítico ] [ algoritmo ] [ gráfica ]
..


6. Tarea

Usar el algoritmo con los siguientes ejercicios y observar las diferencias al usar Sympy al determinar los parámetros. Crear nuevos Símbolos de ser necesario

x = 10*sym.cos(2*pi*440*t-0.4*pi)
x = 10*sym.cos(DosPi*440*t-0.4*pi)
x = 20*sym.cos(DosPi*40*t-0.4*pi)
x = 4*sym.cos(40*pi*t+pi/3)
x = A*sym.sin(w*t+p)

Señal [ sinusoidal ] [ ejercicio ] [ analítico ] [ algoritmo ] [ gráfica ]