4.3 FIR – Respuesta de frecuencia y superposición

[ ejercicio bk ] [ algoritmo_bk ] [ ejercicio_xn ] [ algoritmo_xn ]
..

Para sistemas lineales invariantes en el tiempo, LTI, las entradas senoidales complejas tienen como respuesta otra señal senoidal de la misma frecuencia pero con diferente magnitud y fase. TODOS los sistemas LTI tienen la propiedad «entrada-sinusoide» genera una «salida-sinusoide».

Para un filtro FIR LTI, ante una entrada discreta en tiempo tipo señal exponencial compleja, la salida también es una señal discreta en tiempo, tipo exponencial compleja con diferente amplitud, de la misma frecuencia.

|H(e^{j \omega})| = \sum _{k=0}^{M} b_k e^{-j\omega k} = \sum _{k=0}^{M} h[k] e^{-j\omega k} y[n] = |H(e^{j \omega})| x[n] \text{ ; para } x[n] = A e^{j\varphi}e^{-j\omega n}

..


1. Ejercicio bk

Referencia: McClellan Ejemplo 6.1 y 6.2 p217

Considere un sistema LTI con coeficientes de la ecuación de diferencias bk = [1,2,1]. Al convertirlo en la función respuesta de frecuencias del sistema H(ω) se tiene:

H(e^{j \omega}) = 1 + 2 e^{-j\omega} + e^{-j \omega 2}

como bk es simétrico, se obtiene como factor común e(-jω) para simplificar la expresió con un coseno.

= e^{-j \omega} \Big( e^{j \omega} + 2 + e^{-j \omega} \Big) = e^{-j \omega} \Big( 2 +2 \cos(\omega) \Big)

donde la magnitud es:

|H(e^{j \omega})|= \Big( 2 +2 \cos(\omega) \Big)

y la fase:

\angle H(e^{j \omega})= e^{-j \omega}

Al disponer de una señal de entrada con una frecuencia de ω=π/3:

x[n] = 2 e^{j \pi /4} e^{j(\pi /3)n}

se evalúa la magnitud  de |H|:

|H(e^{j \pi /3})|= \Big( 2 +2 \cos(\pi /3) \Big) = 3

y la fase de ∠H:

\angle H(e^{j \pi /3 })= e^{-j \pi/3}

siendo y[n] = (H)(x[n])

y[n] = 3 e^{-j\pi /3}\Big( 2 e^{j \pi /4} e^{j(\pi /3)n} \Big) = (3)(2) e^{j \pi /4-j\pi /3} e^{j\pi /3 n} = 6 e^{-j \pi /12} e^{j\pi /3 n} = 6 e^{j \pi /4} e^{j\pi(n-1) /3}

la salida del sistema es igual a la entrada, multiplicada por 3 y un cambio de fase a -π/3 que es un retraso de una muestra. Recordando que esto aplica cuando bk es simétrico y la entrada es una exponencial compleja.

[ ejercicio bk ] [ algoritmo_bk ] [ ejercicio_xn ] [ algoritmo_xn ]
..


2. Algoritmo bk en Python

Primero se analiza la simetría del vector bk, invirtiendo el vector y restando de su forma original, si todos los resulados son ceros, el vector es simétrico.

Acorde con la simetría de bk se procede a la creación de H y luego yn a partir de xn.

Para xn en el algoritmo, se requiere usar la dsp.euler_args_one_term() para dominio del tiempo contínuo t de la Unidad 1.4 Espectro – Suma de sinusoides y la fórmula de Euler. La función en telg1034.py amplía el uso al dominio de tiempo discreto n al sustituir la variable independiente con las siguientes instrucciones y se obtiene el resultado esperado.

    # revisa si unasenal es discreta con n, Unidad 4 Filtros Digitales
    esdiscreta = False
    varindepend = X.free_symbols
    if (n in varindepend) and not(t in varindepend):
        X = X.subs(n,t)
        esdiscreta = True
    # analiza unasenal con t

El resultado con el algoritmo es:

bk_simetrico: True
H_magn  : exp(I*w) + 2 + exp(-I*w)
H_magcos: 2*cos(w) + 2
H_fase  : -w
  H : (2*cos(w) + 2)*exp(-I*w)
x[n]: 2*exp(I*pi/4)*exp(I*pi*n/3)
wk  : pi/3
y[n]: 6*exp(I*pi/4)*exp(I*pi*(n - 1)/3)
>>>

Instrucciones en Python

# ejemplo 6.1 p217 FIR - Respuesta de frecuencias
# 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
# variables discretas
from telg1034 import n

# INGRESO
bk = [1,2,1]
#bk = [1,-2,4,-2,1]

xn = 2*sym.exp(I*pi/4)*sym.exp(I*(pi/3)*n)

# PROCEDIMIENTO
M = len(bk)-1
# revisa simetria bk
bk = np.array(bk)
H_factor = 0
sumainvertida = np.sum(np.abs(bk-np.flip(bk)))
bk_simetrico = False
if sumainvertida == 0 and M>1:
    bk_simetrico = True
    H_factor = sym.sympify(M)/2

# respuesta de frecuencia
H = sym.S.Zero
for k in range(0,M+1,1):
    H = H + bk[k]*sym.exp(-I*w*(k))

# H cuando bk es simetrico.
H_magn = H ; H_magcos = 0 ; H_fase = 0
if bk_simetrico:
    H_magn = sym.expand(H/sym.exp(-I*w*H_factor))
    H_magcos = H_magn.rewrite(sym.cos)
    H_fase = -w*H_factor
    H = H_magcos*sym.exp(I*H_fase)

# señal de entrada x[n]
xn_arg = dsp.euler_args_one_term(xn)
wk = xn_arg['freq_rad']

# señal de salida y[n]
if bk_simetrico:
    n_veces = H_fase/w
    Ampl = H_magcos.subs(w,wk)*xn_arg['amplitud'] 
    yn = Ampl*sym.exp(I*xn_arg['fase'])*sym.exp(I*wk*(n+n_veces))
else:
    yn = xn*H.subs(w,wk)
    
# SALIDA
print('bk_simetrico:',bk_simetrico)
print('H_magn  :',H_magn)
print('H_magcos:',H_magcos)
print('H_fase  :',H_fase)
print('  H :',H)
print('x[n]:',xn)
print('wk  :',wk)
print('y[n]:',yn)

[ ejercicio bk ] [ algoritmo_bk ] [ ejercicio_xn ] [ algoritmo_xn ]
..


3. Ejercicio x[n] de varias sinusoides

Referencia: McClellan Ejemplo 6.4 p220

Para el filtro FIR con coeficientes b[k] = [1,2,1] , encuentre la señal de salida cuando la señal de entrada tiene varios componentes.

x[n] = 4 + 3 \cos \Big( \frac{\pi}{3} n - \frac{pi}{2} \Big) + 3 \cos \Big(\frac{7\pi}{8}n \Big)

El ejercicio se desarrolla término a término por tener diferentes frecuencias. Se aplica el principio de superposición.

[ ejercicio bk ] [ algoritmo_bk ] [ ejercicio_xn ] [ algoritmo_xn ]
..


4. Algoritmo X[n] en Python

Para el ejercicio, xn tiene una señal senoidal mas general y de varios componentes. Por lo que el procedimiento del algoritmo anterior se amplía para cada uno de los componentes.

La señal de entrada se debe revisar sus componentes sumas, procediendo de forma semejante a lo realizado en la Unidad 1.2 Señal – Periodo y Desplazamiento en tiempo. La función dsp.x_list_term_Add() para separar los términos suma entrega su resultado tanto para la variable t ó n.

Sin embargo la obtención de los parámetros de cada término suma de xn con la función dsp.cos_args_one_term(),  debe ampliar la función de t a n, semejante al ejercicio que precede con las siguientes instrucciones.

    # revisa si unasenal es discreta con n, Unidad 4. Filtros digitales
    esdiscreta = False
    varindepend = x.free_symbols
    if (n in varindepend) and not(t in varindepend):
        x = x.subs(n,t)
        esdiscreta = True
    # analiza unasenal con t

Aplicando H en cada término xn acorde a la simetría de bk, se tiene el siguiente resultado con el algoritmo es:

H_magn  : exp(I*w) + 2 + exp(-I*w)
H_magcos: 2*cos(w) + 2
H_fase  : -w
  H : (2*cos(w) + 2)*exp(-I*w)
x[n]: 3*cos(7*pi*n/8) + 3*cos(pi*n/3 - pi/2) + 4
y[n]: 0.456722804932279*cos(pi*(7*n/8 - 7/8)) + 9*cos(pi*(n - 1)/3 - pi/2) + 16
>>>

Instrucciones en Python

# ejemplo 6.4 p220 FIR - Respuesta de frecuencias a varias senoidales
# 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
# variables discretas
from telg1034 import n

# INGRESO
bk = [1,2,1]
#bk = [1,-2,4,-2,1]

#xn = 3*sym.cos(n*pi/3 - pi/2)
xn = 4 + 3*sym.cos(n*pi/3 - pi/2) + 3*sym.cos(7*pi*n/8)

# PROCEDIMIENTO
M = len(bk)-1
# revisa simetria bk
bk = np.array(bk)
H_factor = 0
sumainvertida = np.sum(np.abs(bk-np.flip(bk)))
bk_simetrico = False
if sumainvertida == 0 and M>1:
    bk_simetrico = True
    H_factor = sym.sympify(M)/2

# respuesta de frecuencia
H = sym.S.Zero
for k in range(0,M+1,1):
    H = H + bk[k]*sym.exp(-I*w*(k))

# H cuando bk es simetrico.
H_magn = H ; H_magcos = 0 ; H_fase = 0
if bk_simetrico:
    H_magn = sym.expand(H/sym.exp(-I*w*H_factor))
    H_magcos = H_magn.rewrite(sym.cos)
    H_fase = -w*H_factor
    H = H_magcos*sym.exp(I*H_fase)

# señal de entrada x[n]
x_senal = dsp.x_list_term_Add(xn) # términos suma
x_conteo = len(x_senal)
yn = sym.S.Zero
for i in range(0,x_conteo,1):
    unasenal = x_senal[i]
    cond1 = unasenal.has(sym.cos) or unasenal.has(sym.sin)
    cond2 = not(unasenal.has(n))
    if cond1 or cond2:
        xn_arg = dsp.cos_args_one_term(unasenal)
        Ak = xn_arg['amplitud']
        wk = xn_arg['freq_rad']
        fk = xn_arg['fase']
        if bk_simetrico:
            H_magnk = dsp._float_is_int(H_magcos.subs(w,wk).evalf())
            H_fasek = H_fase.subs(w,wk)
            if wk!=sym.S.Zero:
                n_veces = H_fasek/wk
                yk = H_magnk*Ak*sym.cos(wk*(n+n_veces) + fk)
            else:
                yk = H_magnk*Ak
        else:
            yk = H.subs(w,wk)*unasenal
    else: # no es senoidal o constante
        yk = H*unasenal
        print(' ____revisar término:',unasenal.subs(t,n))
    # principio de superposición
    yn = yn +yk

# SALIDA
print('H_magn  :',H_magn)
print('H_magcos:',H_magcos)
print('H_fase  :',H_fase)
print('  H :',H)
print('x[n]:',xn)
print('y[n]:',yn)

[ ejercicio bk ] [ algoritmo_bk ] [ ejercicio_xn ] [ algoritmo_xn ]

4.2 FIR – Filtros y Convolución

[ convolución ] [ ejercicio ] [ algoritmo ]
..


Referencia: McClellan 5.4.3 p181

Una expresión general para salidas de Filtros FIR se obtienen de los términos de respuesta al impulso. Si los coeficientes de bk son los mismos valores de respuesta al impulso se puede expresar como una suma de convolución, que es una operación entre dos secuencias.

y[n] = \sum_{k=0}^{M} h[k] x[n-k] y[n] =h[k] * x[n] =\sum_{k=-\infty}^{\infty} h[k] x[n-k]

Lo que expresa que la secuencia de salida y[n] se obtiene al aplicar la convolución de la respuesta al impulso h[n] con la señal de entrada x[n]

[ convolución ] [ ejercicio ] [ algoritmo ]
..


1. Ejercicio

Referencia: McClellan 5-4.3.3 p184

Evaluar la convolución de hn que es una secuencia de 11 puntos, con xn que es una secuencia de 51 puntos. Use el vector hn como respuesta al impulso de un sistema de medias móviles de 11 puntos.

x[n]= \sin(0.7*\pi*n) \text{ ; } 0 \le n \le 50 h[n]= \begin{cases} 1/11 & 0 \le n \le 40 \\0 & n<0\end{cases}

y[n] = x[n]*h[n]

Al aplicar el algoritmo de media móvil se compara el resultado con el algoritmo de convolución.

muestreo: 51  ; tamaño ki: 73
x[n]: (Heaviside(n) - Heaviside(n - 51))*sin(0.07*pi*n)
h[n]: Heaviside(n)/11 - Heaviside(n - 11)/11
>>>

FIR media movil Convolución[ convolución ] [ ejercicio ] [ algoritmo ]

..


2. Algoritmo en Python

Referencia: [1] LTI DT – Convolución de sumas – Python. Blog Señales y Sistemas. [2] https://numpy.org/doc/stable/reference/generated/numpy.convolve.html

# ejemplo 5-4.3.3 p184 FIR - Media Móvil  y convolución
# 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
# variables discretas
from telg1034 import n

# INGRESO
u = sym.Heaviside(n)
# señal como suma de las partes
muestrasN = 50 + 1  # # intervalo [0,50]
gn = (u-u.subs(n,n-muestrasN)) 
xn = sym.sin(0.07*pi*n)*gn

# FIR Media Móvil
ventana = 10 + 1 # intervalo [0,10]
hn = (u-u.subs(n,n-ventana))/ventana 

# PROCEDIMIENTO ----------------
# muestreo x[n]
ki = np.arange(-ventana,muestrasN + ventana,1,dtype=int)
xni = sym.lambdify(n,xn, modules=equivalentes)
xki = xni(ki)

hni = sym.lambdify(n,hn, modules=equivalentes)
hki = hni(ki) 

# FIR Media Móvil
yki = np.convolve(hki,xki)
yki = yki[ventana:len(xki)+ventana] # intervalo ki

# SALIDA
print('muestreo:',muestrasN,' ; tamaño ki:',len(ki))
print('x[n]:',xn)
print('h[n]:',hn)

# GRAFICA FIR
plt.subplot(311) # x[n]
plt.stem(ki,xki,linefmt = 'C0:',
         label='x[n]')
plt.xlabel('n')
plt.ylabel('x[n]')
plt.title('FIR Media Móvil con convolución')
plt.legend()

plt.subplot(312) # h[n]
plt.stem(ki,hki,linefmt = 'C1:',label='h[n]')
plt.xlabel('n')
plt.ylabel('h[n]')
plt.legend()

plt.subplot(313) # y[n]
plt.stem(ki,yki,linefmt = 'C2:',label='y[n]')
plt.xlabel('n')
plt.ylabel('y[n]')
plt.legend()
plt.show()

[ convolución ] [ ejercicio ] [ algoritmo ]

4.1 FIR – Filtro de Medias Móviles

[ media móvil ] | ejercicio: [ rectangular ] [ 2 componentes ] [ algoritmo ]
..


Referencia: McClellan 5 p167

Se introduce el análisis para un sistema discreto en tiempo con respuesta de impulso finita,  FIR (Finite Impulse Response), también conocido como filtros FIR.

Un filtro es un sistema con valores de salida resultantes de la suma ponderada de un número finito de valores de una secuencia de entrada. Se define como una operación también conocida como ecuación de diferencias.

1. El filtro de Medias Móviles

Una señal discreta se modifica por medio de un filtro que se determina como un promedio de un intervalo que se desplaza en el tiempo. Para un sistema tipo causal, se tiene:

y[n] = \sum_{k=n-M}^{n} b_{n-k} x[k] = b_{(M)} x[n-M ] + b_{(M-1)} x[n-M+1] + \text{...} + b_{(0)} x[n]

cuando los coeficientes de bk no son iguales, se dice que el promedio o media  móvil es ponderado.

[ media móvil ] | ejercicio: [ rectangular ] [ 2 componentes ] [ algoritmo ]
..


2. Ejemplo de señal rectangular

Referencia: McClellan ejemplo 5.1 p174

Considere un sistema con media móvil para 3 puntos que se expresa como una ventana que desplaza.

y[n] = \sum_{k=n-2}^{n} \frac{1}{3} x[k]

con entrada rectangular entre [0,10]

x[n]= \begin{cases} 1, & 0 \le n \le 10 \\0 , & en otro caso\end{cases} FIR media móvilRef: Señales Escalón unitario μ(t) e Impulso unitario δ(t). blog de Señales y Sistemas

La señal se define como un escalón que se compone de dos señales escalón:

x[n]= \mu [n] -\mu [n-11]

Para el algoritmo se define a partir de μ[n] que en Sympy es Heaviside. Se indican la cantidad de muestras a observar.

# INGRESO
u = sym.Heaviside(n)
# señal como suma de las partes
muestrasN = 10 + 1  # # intervalo [0,10]
xn = u - u.subs(n,n-muestrasN)

# FIR media móvil
ventana = 2 + 1 # intervalo [0,2]
bk = np.ones(ventana, dtype=float)/ventana
causal = True  # causal o no causal

En la media móvil se requiere M=2, si el sistema es causal y el vector bk que son los valores de ponderación. Se obtienen los siguientes resultados:

muestreo: 11  ; tamaño ki: 17
x[n]: Heaviside(n) - Heaviside(n - 11)
bk: [0.33333333 0.33333333 0.33333333]

[ media móvil ] | ejercicio: [ rectangular ] [ 2 componentes ] [ algoritmo ]
..


3. Ejercicio  con señal de dos componentes

Considere la señal con un componente exponencial y otro componente senoidal que se considera ruido o interferencia que degrada la señal a observar. Se pretende minimizar o remover el efecto del componente senoidal

x[n]= \begin{cases} (1-02)^n + \frac{1}{2}\cos(2 \pi n /8 + \pi/4), & 0 \le n \le 40 \\0 , & n<0\end{cases}

Realice el ejercicio con media móvil que tiene una ventana de 7 puntos o M=6.

Usando el algoritmo se tiene el siguiente resultado:

FIR media movil 02

muestreo: 41  ; tamaño ki: 55
x[n]: (1.02**n + 0.5*cos(pi*n/4 + pi/4))*(Heaviside(n) - Heaviside(n - 40))
bk: [0.14285714 0.14285714 0.14285714 0.14285714 
     0.14285714 0.14285714 0.14285714]

El bloque de ingreso para el algoritmo es:

# INGRESO
u = sym.Heaviside(n)
# señal como suma de las partes
muestrasN = 40 + 1  # # intervalo [0,40]
x0 = u - u.subs(n,n-muestrasN)
x1 = 1.02**n +0.5*sym.cos(2*pi*n/8 + pi/4)
xn = x1*x0

# FIR media móvil
ventana = 6 + 1 # intervalo [0,6]
bk = np.ones(ventana, dtype=float)/ventana
causal=True  # causal o no causal

[ media móvil ] | ejercicio: [ rectangular ] [ 2 componentes ] [ algoritmo ]
..


4. Algoritmo en Python

El modelo de algoritmo para los ejercicios  se presenta para la señal rectangular. Para otros ejercicios se actualiza la señal en el bloque de ingreso.

La variable causal=True establece el signo de la ventana deslizante, si el sistema es causal, el signo es negativo.

# ejemplo 5.1 p174 FIR - Media Móvil general
# 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
# variables discretas
from telg1034 import n

# INGRESO
u = sym.Heaviside(n)
# señal como suma de las partes
muestrasN = 10 + 1  # # intervalo [0,10]
xn = u - u.subs(n,n-muestrasN)

# FIR media móvil
ventana = 2 + 1 # intervalo [0,2]
bk = np.ones(ventana, dtype=float)/ventana
causal = True  # causal o no causal

# PROCEDIMIENTO ----------------
# muestreo x[n]
ki = np.arange(-ventana,muestrasN + ventana,1,dtype=int)
xni = sym.lambdify(n,xn, modules=equivalentes)
xki = xni(ki) 

# FIR Media Móvil
causal_sign = +1 # sistema causal
if not(causal): # sistema no causal
    causal_sign = -1
    
# ajuste de intervalo con ventana
yki = np.zeros(muestrasN + 2*ventana,dtype=float)
for i in range(-ventana,muestrasN + 2*ventana,1):
    unaventana = 0
    for k in range(0,ventana,1):
        if (i-causal_sign*k)>=0 and (i-causal_sign*k)<=(muestrasN+ventana):
            unaventana = unaventana + xki[i-causal_sign*k]*bk[k]
    yki[i] = unaventana

# SALIDA
print('muestreo:',muestrasN,' ; tamaño ki:',len(ki))
print('x[n]:',xn)
print('bk:',bk)

# GRAFICA FIR
plt.subplot(211) # x[n]
plt.stem(ki,xki,linefmt = 'C0:',
         label='x[n]')
plt.xlabel('n')
plt.ylabel('x[n]')
plt.title('FIR Media Móvil')
plt.legend()

plt.subplot(212) # y[n]
plt.stem(ki,yki,linefmt = 'C2:',label='y[n]')
plt.xlabel('n')
plt.ylabel('y[n]')
plt.legend()
plt.show()

[ media móvil ] | ejercicio: [ rectangular ] [ 2 componentes ] [ algoritmo ]

2.3 Reconstrucción, sobre-muestreo, sub-muestreo y Nyquist

[ reconstrucción ] [ sobre-muestreo ] [ sub-muestreo ] [ algoritmo ] [ gráfica ]
..

1. Muestreo y Reconstrucción

La salida de un convertidor contínuo a discreto C-D es una señal discreta en tiempo con un número infinito de alias. Al observar un gráfico de espectro de frecuencias ω es necesario recordar que existen más señales fuera del intervalo mostrado.

Considere una señal x(t) simple, de una sola frecuencia, cuyo espectro contiene dos líneas a ±ω0 con amplitudes de (1/2 A e±ω0)

x(t) = A cos\Big(\omega_0 t + \varphi \Big) x[n] = x(n/f_s) A cos((\omega_0 /f_s)*n + \varphi = \frac{1}{2}A e^{j\varphi}e^{j(\omega_0/f_s)n} + \frac{1}{2}A e^{-j\varphi}e^{j(-\omega_0/f_s)n}

también tiene dos frecuencias ω en ±ω0/fs, con todos los alias en
±ω0/fs + 2πk, siendo k= ±1, ±2, ±3, …

[ reconstrucción ] [ sobre-muestreo ] [ sub-muestreo ] [ algoritmo ] [ gráfica ]
..


2. Sobre-muestreo

Referencia: McClellan ejercicio 4.5 p133

Realice el espectro contínuo de la señal x(t), la señal muestreada a fs=500 y el espectro discreto resultante. Observe que la frecuencia de muestreo es superior a la frecuencia de la señal x(t) para evitar el problema de aliasing. El proceso se conoce como sobre-muestreo. Para el ejercicio se usa una frecuencia 2.5 veces la recomendada por Nyquist.

x(t) = cos\Big(2\pi (100) t + \pi /3 \Big)

espectro n coseno sobremuestreo

Los convertidores D-C trasforman el espectro discreto en tiempo a una salida de espectro contínuo en el tiempo, seleccionando solo un par de las líneas de todas las posibilidades mostradas.

Para ser consistentes con la operación Digital hacia Analógico (D-A) se tomará siempre la frecuencia mas baja de cada grupo de alias o alias principal con |ω|<π.
El resultado siempre se encontrará entre [-fs/2, +fs/2].

En resumen, en sobre-muestreo la frecuencia original f0 es menor que fs/2 permite la reconstrucción mas exacta.

Para el análisis con el algoritmo se obtienen los siguientes parámetros para realizar la gráfica.

señal(t):   cos(100*t*(2*pi) + pi/3)
 espectro:
  freq : [-100.  100.]
  ampl : [1/2 1/2]
  fase : [-pi/3 pi/3]
  etiq : ['$\\frac{1}{2}$$ e^j(- \\frac{\\pi}{3})$'
 '$\\frac{1}{2}$$ e^j(\\frac{\\pi}{3})$']
  freq_max : 100.0
  freq_min : 100.0
  x_expand : cos(2*pi*n/5 + pi/3)

muestreo_tipo:  sobre-muestreo
muestreo: 11
 tamaño ki: 11 ; tamaño ti: 121
fs: 500  ; dt: 0.002  ; fs_veces: 12
x(t): cos(100*t*(2*pi) + pi/3)
x[n]: cos(n*(2*pi)/5 + pi/3)

señal[n]:   cos(n*(2*pi)/5 + pi/3)
 espectro_n:
  freq : [-2.4 -1.6 -0.4  0.4  1.6  2.4]
  ampl : [1/2 1/2 1/2 1/2 1/2 1/2]
  fase : [-pi/3 pi/3 -pi/3 pi/3 -pi/3 pi/3]
  etiq : ['$\\frac{1}{2}$$ e^j(\\frac{\\pi}{3})$'
 '$\\frac{1}{2}$$ e^j(- \\frac{\\pi}{3})$'
 '$\\frac{1}{2}$$ e^j(- \\frac{\\pi}{3})$'
 '$\\frac{1}{2}$$ e^j(\\frac{\\pi}{3})$'
 '$\\frac{1}{2}$$ e^j(- \\frac{\\pi}{3})$'
 '$\\frac{1}{2}$$ e^j(\\frac{\\pi}{3})$']
  freq_factor : pi
  freq_max : 2.4
  freq_min : 0.4
  BW : 2.0
  freq_alias : [ True  True False False  True  True]

señal reconstruida con alias principal:
cos(100*t*(2*pi) + pi/3)

Las instrucciones del bloque de ingreso para el algoritmo son:

# INGRESO
# usar np.pi para evitar frecuencia simplificada por sympy.
x = sym.cos(DosPi*100*t + pi/3)
x_etiqueta = 'x(t)' ; x_etiqueta_n = 'x[n]'

# señal muestreada
fs = 500  # muestreo discreto, freq_sampling
fs_veces = 12 # sobremuestreo para x(t)
muestrasN = 10 + 1  # para x[n]

aliasn = 1 # alias añadidos a x[n]
tolera = 1e-9 # tolerancia a error, números iguales

[ reconstrucción ] [ sobre-muestreo ] [ sub-muestreo ] [ algoritmo ] [ gráfica ]
..


3. Sub-muestreo y aliasing

Cuando fs < 2f0 la señal se encuentra sub-muestreada y se produce el aliasing. Para el ejercicio anterior, si el ejercicio lo realiza con fs = 80Hz, siendo f0=100 Hz, la tasa de muestreo es menor que la frecuencia de Nyquist y se presenta la señal de alias.

espectro n coseno submuestreo

las instrucciones para el bloque de ingreso del algoritmo son las mismas que en el ejercicio anterior, excepto que fs = 80 Hz

# señal muestreada
fs = 80  # muestreo discreto, freq_sampling

el resultado con el algoritmo es una señal reconstruida de tan solo 20Hz

x_0(t) = cos\Big(2\pi (20) t + \pi /3 \Big)
señal(t):   cos(100*t*(2*pi) + pi/3)
 espectro:
  freq : [-100.  100.]
  ampl : [1/2 1/2]
  fase : [-pi/3 pi/3]
  etiq : ['$\\frac{1}{2}$$ e^j(- \\frac{\\pi}{3})$'
 '$\\frac{1}{2}$$ e^j(\\frac{\\pi}{3})$']
  freq_max : 100.0
  freq_min : 100.0
  x_expand : cos(5*pi*n/2 + pi/3)

muestreo_tipo:  sub-muestreo
muestreo: 11
 tamaño ki: 11 ; tamaño ti: 121
fs: 80  ; dt: 0.0125  ; fs_veces: 12
x(t): cos(100*t*(2*pi) + pi/3)
x[n]: cos(5*n*(2*pi)/4 + pi/3)

señal[n]:   cos(5*n*(2*pi)/4 + pi/3)
 espectro_n:
  freq : [-4.5 -2.5 -0.5  0.5  2.5  4.5]
  ampl : [1/2 1/2 1/2 1/2 1/2 1/2]
  fase : [-pi/3 -pi/3 -pi/3 pi/3 pi/3 pi/3]
  etiq : ['$\\frac{1}{2}$$ e^j(\\frac{\\pi}{3})$'
 '$\\frac{1}{2}$$ e^j(- \\frac{\\pi}{3})$'
 '$\\frac{1}{2}$$ e^j(- \\frac{\\pi}{3})$'
 '$\\frac{1}{2}$$ e^j(- \\frac{\\pi}{3})$'
 '$\\frac{1}{2}$$ e^j(\\frac{\\pi}{3})$'
 '$\\frac{1}{2}$$ e^j(\\frac{\\pi}{3})$']
  freq_factor : pi
  freq_max : 4.5
  freq_min : 0.5
  BW : 4.0
  freq_alias : [ True  True False False  True  True]

señal reconstruida con alias principal:
cos(20*t*(2*pi) + pi/3)

[ reconstrucción ] [ sobre-muestreo ] [ sub-muestreo ] [ algoritmo ] [ gráfica ]
..


4. Folding o plegado de señal por sub-muestreo

Para una señal x(t) de 100 Hz se aplica muestreo insuficiente en el intervalo [f0,2f0] se presenta aliasing denominado «folding» o plegado de señal. Para el ejercicio de prueba y con fs = 125, se tiene:

x(t) = cos\Big(2\pi (100) t + \pi /3 \Big)

espectro n coseno folding
Los componentes de frecuencia entre ±π también son ±0.4π, sin embargo  el componente en +0.4π es un alias del componente de la frecuencia negativa -1.6π, que genera el «plegado» o «folding». La reconstrucción de la señal se realiza a f = 0.4(fs/2π) = fs/5 = 25 Hz. También se observa una inversión de fase. El resultado es semejante a una señal muestreada de 25 Hz con inversión de fase.

señal(t):   cos(100*t*(2*pi) + pi/3)
 espectro:
  freq : [-100.  100.]
  ampl : [1/2 1/2]
  fase : [-pi/3 pi/3]
  etiq : ['$\\frac{1}{2}$$ e^j(- \\frac{\\pi}{3})$'
 '$\\frac{1}{2}$$ e^j(\\frac{\\pi}{3})$']
  freq_max : 100.0
  freq_min : 100.0
  x_expand : cos(8*pi*n/5 + pi/3)

muestreo_tipo:  sub-muestreo y folding
muestreo: 11
 tamaño ki: 11 ; tamaño ti: 121
fs: 125  ; dt: 0.008  ; fs_veces: 12
x(t): cos(100*t*(2*pi) + pi/3)
x[n]: cos(4*n*(2*pi)/5 + pi/3)

señal[n]:   cos(4*n*(2*pi)/5 + pi/3)
 espectro_n:
  freq : [-3.6 -1.6 -0.4  0.4  1.6  3.6]
  ampl : [1/2 1/2 1/2 1/2 1/2 1/2]
  fase : [-pi/3 -pi/3 pi/3 -pi/3 pi/3 pi/3]
  etiq : ['$\\frac{1}{2}$$ e^j(\\frac{\\pi}{3})$'
 '$\\frac{1}{2}$$ e^j(- \\frac{\\pi}{3})$'
 '$\\frac{1}{2}$$ e^j(- \\frac{\\pi}{3})$'
 '$\\frac{1}{2}$$ e^j(- \\frac{\\pi}{3})$'
 '$\\frac{1}{2}$$ e^j(\\frac{\\pi}{3})$'
 '$\\frac{1}{2}$$ e^j(\\frac{\\pi}{3})$']
  freq_factor : pi
  freq_max : 3.6
  freq_min : 0.3999999999999999
  BW : 3.2
  freq_alias : [ True  True False False  True  True]

señal reconstruida con alias principal:
cos(25.0*t*(2*pi) - (pi/3))

Otro caso a observar es cuando fs=100Hz, el resultado es una constante, pues siempre se obtiene la muestra de la señal con la misma magnitud.

espectro n coseno folding02

5. Algoritmo en Python

El algoritmo contiene componentes desarrollados en secciones anteriores. Para el espectro de frecuencias en dominio del tiempo se basa en lo descrito en Espectro – Operaciones en dominio de tiempo y frecuencia.

El desarrollo del muestreo usando la señal x(t) y fs toma como referencia lo descrito en Muestreo y aliasing.

La construcción del espectro discreto a partir de x[n] usa lo desarrollado en Espectro x[n] – Señal discreta en tiempo

# ejercicio 4.5 p133 Muestreo y Reconstrucción de señales sinusoidales
# 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
# variables discretas
from telg1034 import n

# INGRESO
# usar np.pi para evitar frecuencia simplificada por sympy.
x = sym.cos(DosPi*100*t + pi/3)
x_etiqueta = 'x(t)' ; x_etiqueta_n = 'x[n]'

# señal muestreada
fs = 500  # muestreo discreto, freq_sampling
fs_veces = 12 # sobremuestreo para x(t)
muestrasN = 10 + 1  # para x[n]

aliasn = 1 # alias añadidos a x[n]
tolera = 1e-9 # tolerancia a error, números iguales

# PROCEDIMIENTO ----------------
# operaciones para espectro x(t)
# Ref: Blog|Espectro-Operaciones en dominio de tiempo y frecuencia
x_term = dsp.x_list_term_Add(x)
Xe_term = dsp.cos_to_euler_one_term(x_term)
x_term_expand = dsp.euler_to_cos_list(Xe_term)
un_espectro = dsp.cos_spectrum_list(x_term)
# espectro de cada señal
freq = un_espectro['freq']
ampl = un_espectro['ampl']
fase = un_espectro['fase']
etiqueta = un_espectro['etiq']

# freq_max y freq_min
freq_max = float(max(freq))
if len(freq[freq>0])>0:
    freq_min = float(min(freq[freq>0]))
else:
    freq_min = 0
un_espectro['freq_max'] = freq_max
un_espectro['freq_min'] = freq_min

# revisa muestreo Nyquist
if freq_max<=(fs/2):
    muestreo_tipo = 'sobre-muestreo'
if freq_max>(fs):
    muestreo_tipo = 'sub-muestreo'
if freq_max>(fs/2) and freq_max<=fs:
    muestreo_tipo = 'sub-muestreo y folding'
    
# Muestreo de señales sinusoidales
# Ref: blog|2.1 Muestreo y aliasing
xn = x.subs(t,n/fs)
# muestreo x[n]
dt = 1/fs 
ki = np.arange(0,muestrasN,1,dtype=int)
# muestreo x(t)
dtc = dt/fs_veces # para x(t)
ti = np.arange(0,ki[muestrasN-1]*dt + dtc, dtc)
tki = ki*dt

xt = sym.lambdify(t,x, modules=equivalentes)
xki = xt(tki) 
xti = xt(ti)

# operaciones para espectro x[n]
# Ref: blog|2.2 Espectro x[n] – Señal discreta en tiempo
x_term = dsp.x_list_term_Add(xn)
Xe_term = dsp.cos_to_euler_one_term(x_term)
x_term_expand = dsp.euler_to_cos_list(Xe_term)
un_espectro_n = dsp.cos_spectrum_list(x_term)
# espectro de cada señal
freq_n = np.array(un_espectro_n['freq'])
ampl_n = np.array(un_espectro_n['ampl'])
fase_n = np.array(un_espectro_n['fase'])
etiq_n = np.array(un_espectro_n['etiq'])
un_espectro['x_expand'] = x_term_expand

# Añadir alias a señal
freqn_conteo = len(freq_n)
freqn_alias = np.zeros(freqn_conteo,dtype=bool)
for i in range(0,freqn_conteo,1):
    freq_i = freq_n[i]
    ampl_i = ampl_n[i]
    fase_i = fase_n[i]
    etiq_i = etiq_n[i]
    for k in range(1,aliasn+1,1):
        # añade freqm[i] + 2*pi*k
        freq_n = np.concatenate([freq_n,[freq_i+2*pi*k]])
        ampl_n = np.concatenate([ampl_n,[ampl_i]])
        fase_n = np.concatenate([fase_n,[fase_i]])
        etiq_n = np.concatenate([etiq_n,[etiq_i]])
        freqn_alias = np.concatenate([freqn_alias,[True]])
        # añade freqm[i] - 2*pi*k
        freq_n = np.concatenate([freq_n,[-(freq_i.evalf()+2*np.pi*k)]])
        ampl_n = np.concatenate([ampl_n,[ampl_i]])
        fase_n = np.concatenate([fase_n,[-fase_i]])
        texto = '$' + sym.latex(ampl_i)+'$'
        if fase_i!=sym.S.Zero:
            texto = texto+f'$ e^j('+sym.latex(fase_i)+')$'
        etiq_n = np.concatenate([etiq_n,[texto]])
        freqn_alias = np.concatenate([freqn_alias,[True]])
# ordena frecuencias
orden = np.argsort(freq_n)
freq_n = freq_n[orden]
ampl_n = ampl_n[orden]
fase_n = fase_n[orden]
etiq_n = etiq_n[orden]
freqn_alias = freqn_alias[orden]

# revisa factor pi en freq_n
unfactor = sym.S.One ; factor_pi = False
freqn_conteo = len(freq_n)
for i in range(0,freqn_conteo,1):
    if freq_n[i].has(pi):
        factor_pi = True
if factor_pi:
    freq_n = np.array(freq_n/pi,dtype=float)
    unfactor = pi
# actualiza espectro de señal con factor pi,alias,orden
un_espectro_n['freq_factor'] = unfactor
un_espectro_n['freq'] = freq_n
un_espectro_n['ampl'] = ampl_n
un_espectro_n['fase'] = fase_n
un_espectro_n['etiq'] = etiq_n

# ancho de banda, freq_max y freq_min
freq_posi = freq_n[freq_n>0]
freq_max_n = float(max(freq_posi))
freq_min_n = 0
if len(freq_posi)>0:
    freq_min_n = float(min(freq_posi))
freq_centro_n = (freq_max_n + freq_min_n)/2
un_espectro_n['freq_max'] = freq_max_n
un_espectro_n['freq_min'] = freq_min_n
un_espectro_n['BW'] = freq_max_n - freq_min_n
    
# aliasing, Revisa espectro de frecuencias  
freq_conteo = len(freq_n)
freq_alias = np.zeros(freq_conteo,dtype=bool)
for i in range (0,freq_conteo,1):
    unafreq = freq_n[i]
    k = 1
    unalias = unafreq + k*2
    while unalias<=freq_max_n:
        for j in range(i+1,freq_conteo,1): 
            if abs(unalias - freq_n[j])<tolera:
                if abs(freq_n[i])>abs(freq_n[j]):
                    freq_alias[i] = True
                else:
                    freq_alias[j] = True
        k = k+1
        unalias = unafreq + k*2
freq_alias_conteo = np.sum(freq_alias)
un_espectro_n['freq_alias'] = freq_alias

# reconstruye señal sin alias
freq_alias0 = np.invert(freq_alias)
freq_x0 = freq_n[freq_alias0]
ampl_x0 = ampl_n[freq_alias0]
fase_x0 = fase_n[freq_alias0]
x0_cuenta = len(freq_x0)
x0t = sym.S.Zero
for i in range(0,x0_cuenta,1):
    if freq_x0[i]>=0:
        una_freq = freq_x0[i]*unfactor*(fs/(2*pi))
        una_freq = dsp._float_is_int(una_freq) # si es entero
        x0t = x0t + 2*ampl_x0[i]*sym.cos(DosPi*una_freq*t+sym.UnevaluatedExpr(fase_x0[i]))        
x0tr = x0t.subs(pi,np.pi)
x0tr = sym.lambdify(t,x0tr,modules=equivalentes)
x0_ti = x0tr(ti)

# SALIDA
print('señal(t):  ',x)
print(' espectro:')
for parametro in un_espectro:
    print(' ',parametro,':',un_espectro[parametro])

print('\nmuestreo_tipo: ',muestreo_tipo)
print('muestreo:',muestrasN)
print(' tamaño ki:',len(ki),'; tamaño ti:',len(ti))
print('fs:',fs,' ; dt:',dt, ' ; fs_veces:',fs_veces)
print('x(t):',x)
print('x[n]:',xn)

print('\nseñal[n]:  ',xn)
print(' espectro_n:')
for parametro in un_espectro_n:
    print(' ',parametro,':',un_espectro_n[parametro])

print('\nseñal reconstruida con alias principal:')
print(x0t)

[ reconstrucción ] [ sobre-muestreo ] [ sub-muestreo ] [ algoritmo ] [ gráfica ]
..


6. Gráfica

Con los componentes descritos en la parte de algoritmos, se ajustan las gráficas de cada sección para cada sub-gráfica en la figura.

# GRAFICAS  ----------------------------------------
# Grafica de espectro de frecuencias
# Ref: blog|2.2 Espectro x[n] – Señal discreta en tiempo
graf_dx = 0.12
fig_espectro = plt.figure()
ampl_max = float(max(ampl))
freq_max = float(max(max(freq),fs))*(1+graf_dx/2)
graf_fasor = fig_espectro.add_subplot(311)
if freq_max!=0:
    graf_fasor.set_xlim([-freq_max*(1+graf_dx),
                         freq_max*(1+graf_dx)])
else:
    graf_fasor.set_xlim([-1,1])
graf_fasor.set_ylim([0,ampl_max*(1+graf_dx*2)])
graf_fasor.grid()
graf_fasor.axhline(0,color='black')
graf_fasor.axvline(0,linestyle='dotted',color='grey')
graf_fasor.stem(freq,ampl) # grafica magnitud
for k in range(0,len(freq),1): # etiquetas de fasor
    graf_fasor.annotate(etiqueta[k],xy=(freq[k],ampl[k]),
        xytext=(0,5),textcoords='offset points',ha='center')
# fs frecuencia de muestreo
graf_fasor.stem(fs,ampl_max/2,linefmt = 'C3:', ) # fs
graf_fasor.annotate('fs',xy=(fs,ampl_max/2),
        xytext=(0,5),textcoords='offset points',ha='center')
graf_fasor.set_ylabel('Magnitud '+x_etiqueta)
graf_fasor.set_xlabel('freq Hz', labelpad=0,)
texto = '$' + sym.latex(x)+'$' +'  ; fs='+str(fs)
texto = texto + '  ; alias='+str(aliasn)
graf_fasor.set_title(texto)

# GRAFICA x(t), x[n], x(t)_alias0
graf_xt = fig_espectro.add_subplot(312)
graf_xt.plot(ti,xti, color='gray',label='x(t)')
graf_xt.plot(ti,x0_ti, color='blue',label='x(t)_alias0',
             linestyle='dotted')
graf_xt.stem(tki,xki,linefmt = 'C0:',label='x[n]')
graf_xt.set_xlabel('t')
graf_xt.set_ylabel('amplitud')
graf_xt.legend()
graf_xt.grid()

# GRAFICAS de espectro de frecuencias discretas---------
# Ref: blog|2.2 Espectro x[n] – Señal discreta en tiempo
freq_n = un_espectro_n['freq']
ampl_n = un_espectro_n['ampl']
etiq_n = un_espectro_n['etiq']
freq_alias = un_espectro_n['freq_alias']
ampl_max = float(max(ampl_n))*(1+graf_dx*2)
freq_max = float(max(un_espectro_n['freq_max'],1))*(1+graf_dx/2)
# grafica
graf_fasorn = fig_espectro.add_subplot(313)
graf_fasorn.set_xlim([-freq_max,freq_max])
graf_fasorn.set_ylim([min([min(ampl_n),0]),ampl_max])
graf_fasorn.grid()
graf_fasorn.axhline(0,color='black')
graf_fasorn.axvline(0,linestyle='dotted',color='grey')
# grafica magnitud
freq_noalias = np.invert(freq_alias) 
graf_fasorn.stem(freq_n[freq_noalias],ampl_n[freq_noalias],
                 linefmt = 'C0--')
graf_fasorn.stem(freq_n[freq_alias],ampl_n[freq_alias],
                 linefmt = 'C1--',label = 'alias') 
if un_espectro_n['freq_factor'] != sym.S.One:
    unfactor = r'$'+sym.latex(un_espectro_n['freq_factor'])+'$'
    freq_etiq = []
    for un_num in freq_n:
        freq_etiq.append(f''+str(round(un_num,4))+unfactor)
    graf_fasorn.set_xticks(ticks=freq_n,labels=freq_etiq)
for k in range(0,len(freq_n),1): # etiquetas de fasor
    graf_fasorn.annotate(etiq_n[k],xy=(freq_n[k],ampl_n[k]),
                 xytext=(0,5),textcoords='offset points',
                 ha='center')
graf_fasorn.set_ylabel('Magnitud '+x_etiqueta_n)
graf_fasorn.set_xlabel('freq rad')
graf_fasorn.legend()

plt.tight_layout()
plt.show()

[ reconstrucción ] [ sobre-muestreo ] [ sub-muestreo ] [ algoritmo ] [ gráfica ]

2.2 Espectro x[n] – Señal discreta en tiempo

[ espectro_n ] [ algoritmo ] [ gráfica ]
..


1. Espectro de frecuencias de x[n]

Referencia: McClellan 4.1.4 p129

De forma semejante al Espectro para Operaciones en dominio de tiempo y frecuencia, se aplica el concepto de frecuencias en radiantes. El espectro permitirá analizar las señales alias: alias principal junto a otros alias generados con ±2π en cada frecuencia. Por ejemplo, para:

x1[n] = cos(0.4π n)

Es señal básica con frecuencia en radianes en 0.4π, sería alias principal, con componentes en ±0.4π. En el dominio discreto ‘n‘, las señales alias se forman con cada componente desplazados en 2πk a la derecha y -2πk a la izquierda. Siendo k un número entero. Para k=1.

Señal / alias principal otros alias, k=1
ω1[n] = 0.4π ω3[n] = 0.4π + 2π = 2.4
ω4[n] = 0.4π – 2π = -1.6
ω2[n] = – 0.4π ω5[n] = -0.4π + 2π = 1.6
ω6[n] = -0.4π – 2π = -2.4

el espectro con señales alias para k=1 es:

espectro n coseno

[ espectro_n ] [ algoritmo ] [ gráfica ]

..


2. Algoritmo en Python

El punto de partida es el algoritmo usado en Espectro para Operaciones en dominio de tiempo y frecuencia, donde las operaciones son semejantes y aplicadas al dominio ‘t‘, que deben cambiarse al domino ‘n‘. Se aplica un cambio de variable y una bandera ‘esdiscreta‘ para ajustar el eje de frecuencias a radianes en el análisis de la función cos_spectrum_list(x_senales) del archivo telg1034.py. con lo que se logra obtener básicamente los resultados del espectro ajustados al nuevo dominio

    #...
    for i in range(0,x_conteo,1):
        unasenal = x_senales[i]
        
        # revisa si unasenal es discreta con n
        esdiscreta = False
        varindepend = unasenal.free_symbols
        if (n in varindepend) and not(t in varindepend):
            unasenal = unasenal.subs(n,t)
            esdiscreta = True
        # analiza unasenal con t
        #...

cuando la señal es discreta, las frecuencias consideran el factor 2π, que se separaba cuando se calculaban en Hz.

            if esdiscreta:
                freq_k = freq_k*2*pi
            x_freq_spectr.append(freq_k)

con lo que el algoritmo se ajusta a la variable y es funcional.

Para la sección gráfica, se puede considerar el factor π como un elemento a simplificar en el vector de frecuencias. Se revisa si existe el factor en el arreglo de frecuencias y se procede a separarlo dentro del algoritmo de espectro para t. Se actualizan los valores en los resultados.

...
# revisa factor pi en freq_n
unfactor = sym.S.One ; factor_pi = False
freqn_conteo = len(freq_n)
for i in range(0,freqn_conteo,1):
    if freq_n[i].has(pi):
        factor_pi = True
if factor_pi:
    freq_n = np.array(freq_n/pi,dtype=float)
    unfactor = pi
# actualiza espectro de señal con factor pi
un_espectro_n['freq_factor'] = unfactor
un_espectro_n['freq'] = freq_n
...

2.1 Revisión de señales «alias» en x_senales en algoritmo para espectro

Para tener un marcador de frecuencias alias en las señales analizadas, se compara cada frecuencia del espectro para diferentes valores de k en la operación ±2π, mientras no sobrepase la frecuencia máxima del intervalo de freq.

...
# aliasing, Revisa espectro de frecuencias 
freq_conteo = len(freq_n)
freq_alias  = np.zeros(freq_conteo,dtype=bool)
for i in range (0,freq_conteo,1):
    unafreq = freq_n[i]
    k = 1
    unalias = unafreq + k*2
    while unalias<=freq_max:
        for j in range(i+1,freq_conteo,1):
            if abs(unalias - freq_n[j])<tolera:
                if abs(freq_n[i])>abs(freq_n[j]):
                    freq_alias[i] = True
                else:
                    freq_alias[j] = True
        k = k+1
        unalias = unafreq + k*2
un_espectro_n['freq_alias'] = freq_alias
...

Se marcan las frecuencias alias en la lista freq_alias para diferenciarlas con otro color en la gráfica de espectro, obteniendo el resultado presentado en el ejercicio.

[ espectro_n ] [ algoritmo ] [ gráfica ]

2.2 Instrucciones en Python

La señal para el espectro se construye con los componentes alias presentados en la parte teórica. Se analiza solo la señal x4[n] para el espectro.

El parámetro tolera en el bloque de ingreso es la tolerancia para redondear los valores de frecuencia en la gráfica cuando tienen varios dígitos decimales.

# ejercicio 4.4 p130 Espectro de señales senoidales discretas
# 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
# variables discretas
from telg1034 import n

# INGRESO
# usar np.pi para evitar frecuencia simplificada por sympy
x1 = sym.cos(0.4*pi*n)
x2 = sym.cos(1.6*np.pi*n)
x3 = sym.cos(2.4*np.pi*n)
xn = x1+x2+x3
x_etiqueta_n = 'x[n]'

fs = 60  # muestreo discreto, freq_sampling
fs_veces = 12 # sobremuestreo para x(t)
muestrasN = 6 + 1  # para x[n]

tolera = 1e-9 # tolerancia a error, números iguales
# PROCEDIMIENTO
# espectro x[n]
# Ref: Blog|Espectro-Operaciones en dominio de tiempo y frecuencia
x_term = dsp.x_list_term_Add(xn)
Xe_term = dsp.cos_to_euler_one_term(x_term)
x_term_expand = dsp.euler_to_cos_list(Xe_term)
un_espectro_n = dsp.cos_spectrum_list(x_term)
un_espectro_n['x_expand'] = x_term_expand
# espectro de cada señal
freq_n = un_espectro_n['freq']
ampl_n = un_espectro_n['ampl']
fase_n = un_espectro_n['fase']
etiq_n = np.array(un_espectro_n['etiq'])

# revisa factor pi en freq_n
unfactor = sym.S.One ; factor_pi = False
freqn_conteo = len(freq_n)
for i in range(0,freqn_conteo,1):
    if freq_n[i].has(pi):
        factor_pi = True
if factor_pi:
    freq_n = np.array(freq_n/pi,dtype=float)
    unfactor = pi
# actualiza espectro de señal con factor pi
un_espectro_n['freq_factor'] = unfactor
un_espectro_n['freq'] = freq_n

# ancho de banda, freq_max y freq_min
freq_posi = freq_n[freq_n>0]
freq_max = float(max(freq_posi))
freq_min = 0
if len(freq_posi)>0:
    freq_min = float(min(freq_posi))
freq_centro = (freq_max+freq_min)/2
un_espectro_n['freq_max'] = freq_max
un_espectro_n['freq_min'] = freq_min
un_espectro_n['BW'] = freq_max-freq_min

# aliasing, Revisa espectro de frecuencias 
freq_conteo = len(freq_n)
freq_alias  = np.zeros(freq_conteo,dtype=bool)
for i in range (0,freq_conteo,1):
    unafreq = freq_n[i]
    k = 1
    unalias = unafreq + k*2
    while unalias<=freq_max:
        for j in range(i+1,freq_conteo,1):
            if abs(unalias - freq_n[j])<tolera:
                if abs(freq_n[i])>abs(freq_n[j]):
                    freq_alias[i] = True
                else:
                    freq_alias[j] = True
        k = k+1
        unalias = unafreq + k*2
un_espectro_n['freq_alias'] = freq_alias

# SALIDA
print('senal[n]:  ',xn)
print(' espectro:')
for parametro in un_espectro_n:
    print(' ',parametro,':',un_espectro_n[parametro])

[ espectro_n ] [ algoritmo ] [ gráfica ]
..


3. Gráfica espectro dominio n

Se añaden la gráfica para el espectro de frecuencias ajustada al dominio n

# GRAFICAS de espectro de frecuencias discretas---------
graf_dx = 0.12
fig_espectro = plt.figure()
freq_alias_conteo = np.sum(freq_alias)
ampl_max = float(max(ampl_n))*(1+graf_dx*2)
freq_max = float(max(un_espectro_n['freq_max'],1))*(1+graf_dx/2)
# grafica
graf_fasorn = fig_espectro.add_subplot(111)
graf_fasorn.set_xlim([-freq_max,freq_max])
graf_fasorn.set_ylim([min([min(ampl_n),0]),ampl_max])
graf_fasorn.grid()
graf_fasorn.axhline(0,color='black')
graf_fasorn.axvline(0,linestyle='dotted',color='grey')
# grafica magnitud
freq_noalias = np.invert(freq_alias) 
graf_fasorn.stem(freq_n[freq_noalias],ampl_n[freq_noalias],
                 linefmt = 'C0--')
if freq_alias_conteo>0: # hay alias
    graf_fasorn.stem(freq_n[freq_alias],ampl_n[freq_alias],
                     linefmt = 'C1--', label='alias') 
if un_espectro_n['freq_factor'] != sym.S.One:
    unfactor = r'$'+sym.latex(un_espectro_n['freq_factor'])+'$'
    freq_etiq = []
    for un_num in freq_n:
        freq_etiq.append(f''+str(round(un_num,4))+unfactor)
    graf_fasorn.set_xticks(ticks=freq_n,labels=freq_etiq)
for k in range(0,len(freq_n),1): # etiquetas de fasor
    graf_fasorn.annotate(etiq_n[k],xy=(freq_n[k],ampl_n[k]),
                         xytext=(0,5),textcoords='offset points',
                         ha='center')
graf_fasorn.set_ylabel(x_etiqueta_n)
graf_fasorn.set_xlabel('freq rad')
graf_fasorn.legend()
texto = '$' + sym.latex(xn)+'$' +'  ; fs='+str(fs)
graf_fasorn.set_title(texto)

plt.show()

[ espectro_n ] [ algoritmo ] [ gráfica ]

2.1 Muestreo y aliasing

[ muestras ] [ alias ] [ ejemplo ] [ algoritmo ] [ ejercicio alias ]
..


1. Muestras de una señal

Referencia: McClellan 4.1 p123

Una señal discreta en tiempo x[n], representa como una secuencia ordenada de números que puede ser almacenada para su uso y procesamiento posterior. Las muestras de una señal contínua o analógica x(t) se encuentran espaciadas en el tiempo separadas por un tiempo Ts.

x[n] = x(nTs) =x(n/fs)   ; -∞ <n < ∞

Cada valor de x[n] se denomina una muestra de la señal contínua x(t). Ts  también se expresa como frecuencia de muestreo fs = 1/Ts  en muestras/s.

1.1 Muestreo de una señal sinusoidal

Para obtener muestras de una señal x(t) = A cos(ωt+φ) se tiene que:

x[n] = x(nTs) = A cos(ωnTs +φ)
= A cos( (ωTs )n+φ)
= A cos(ωrad n+φ)

donde se tiene la frecuencia en radianes normalizada.
La frecuencia normalizada es un parámetro adimensional al quitar la unidad de tiempo de x(t), siendo ‘n‘ el índice de posición en la secuencia de muestras.

1.2. Ejemplo

Referencia: McClellan Figura 4.3 p126

Una señal tipo coseno de 100 Hz se toman muestras con una tasa fs = 500 muestras/s. Realice y observe las gráficas de t y n correspondientes

 x(t) = cos(2π(100)t)

muestreo: 11
 tamaño ki: 11
 tamaño ti: 101
fs: 500  ; dt: 0.002  ; fs_veces: 10
x(t): cos(100*t*(2*pi))
x[n]: cos(n*(2*pi)/5)

De la gráfica se puede observar que sin sobreponer la forma de onda x(t) sobre x[n] no es sencillo discernir la forma exacta de la forma de onda original. Para una mejor observación se sobre-muestrea 10 veces la señal para que sea más semejante a la forma contínua de x(t).

muestreo coseno 100Hz 011.3 Algoritmo

# ejercicio figura 4.1 p126 Muestreo de señales sinusoidales
# 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
# variables discretas
n = sym.Symbol('n', integer=True, nonnegative=True)

# INGRESO
# señal continua
x = sym.cos(100*DosPi*t + 0)

fs = 500  # muestreo discreto
fs_veces = 10 # sobremuestreo para x(t)
muestrasN = 10 + 1  # para x[n]

# PROCEDIMIENTO
xn = x.subs(t,n/fs)
# muestreo x[n]
dt = 1/fs 
ki = np.arange(0,muestrasN,1,dtype=int)
# muestreo x(t)
dtc = dt/fs_veces # para x(t)
ti = np.arange(0,ki[muestrasN-1]*dt + dtc, dtc)
tki = ki*dt

xt = sym.lambdify(t,x, modules=equivalentes)
xki = xt(tki) 
xti = xt(ti)

# SALIDA
print('muestreo:',muestrasN)
print(' tamaño ki:',len(ki))
print(' tamaño ti:',len(ti))
print('fs:',fs,' ; dt:',dt, ' ; fs_veces:',fs_veces)
print('x(t):',x)
print('x[n]:',xn)

# GRAFICA x(t) y x[n]
plt.subplot(211) # x(t) contínua
plt.plot(ti,xti, color='gray',label='x(t)')
plt.stem(tki,xki,linefmt = 'C0:',
         label='x[n]')
plt.xlabel('t')
plt.ylabel('amplitud')
plt.title('x(t) vs x[n]')
plt.legend()

plt.subplot(212) # x[n] discreta
plt.stem(xki,linefmt = 'C0:',label='x[n]')
#plt.xticks(ki)
plt.xlabel('n')
plt.ylabel('amplitud')
plt.legend()
plt.show()

[ muestras ] [ alias ] [ ejemplo ] [ algoritmo ] [ ejercicio alias ]
..


2. Concepto de alias

Una definición simple de «alias» es de dos nombres para dos personas: «Pancho» y «Francisco». Dos señales diferentes en su forma discreta pueden tener la misma secuencia de valores.
..


2.1 Ejemplo

Referencia: McClellan 4.1.2 p127

Considere las siguientes señales y compare sus gráficas.

x1[n] = cos(0.4π n)

x2[n] = cos(2.4π n) = cos(0.4π n + 2π n )  = cos(0.4π n)

dado que 2πn es un número entero de periodos de la señal coseno.

Con Python se realiza las gráficas obteniendo las muestras usando fs para observar los puntos sobre x1(t) de ω=0.4π. Para suavizar la curva x(t) se usa fs_veces como sobremuestreo de la señal y crear la linea azul.
Se compara con x2(t) realizando el mismo proceso, que al unificar las gráficas se observa que los puntos  de muestra son válidos para ambas señales x1(t), x2(t).

muestreo: 13
 tamaño ki: 13
 tamaño ti: 145
fs: 60  ; dt: 0.016666666666666666  ; fs_veces: 12
cos(0.4*pi*n)
cos(7.5398223686155*n)

La gráfica permite mostrar que las muestras son iguales en cada dt=1/fs para ambas señales.

muestreo Alias coseno dos señales
[ muestras ] [ alias ] [ ejemplo ] [ algoritmo ] [ ejercicio alias ]
..


2.2 Instrucciones en Python

 

# ejercicio 4.2 p128 Muestreo de señales sinusoidales
# 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
# variables discretas
from telg1034 import n

# INGRESO
# usar np.pi para evitar frecuencia fundamental automatica de sympy
x1 = sym.cos(0.4*pi*n) 
x2 = sym.cos(2.4*np.pi*n)

x_senales = [x1,x2]
x_etiqueta = ['x1[n]','x2[n]']

fs = 60  # muestreo discreto
fs_veces = 12 # sobremuestreo para x(t)

muestrasN = 12 + 1  # para x[n]

# PROCEDIMIENTO
x_conteo = len(x_senales)
# muestreo x[n]
dt = 1/fs 
ki = np.arange(0,muestrasN,1,dtype=int)
# muestreo x(t)
dtc = dt/fs_veces # para x(t)
ti = np.arange(0,ki[muestrasN-1]*dt+dtc,dtc)
x_muestras = {}
for k in range(0,x_conteo,1):
    unasenal = x_senales[k]
    xtk = unasenal.subs(n,t*fs)
    xt = sym.lambdify(t,xtk, modules=equivalentes)
    xki = xt(ki*dt) 
    xti = xt(ti)
    x_muestras[k] = {'ki':ki,'xki':xki,
                     'ti':ti,'xti':xti,}

# SALIDA
print('muestreo:',muestrasN)
print(' tamaño ki:',len(ki))
print(' tamaño ti:',len(ti))
print('fs:',fs,' ; dt:',dt, ' ; fs_veces:',fs_veces)
for unasenal in x_senales:
    print(unasenal)

# GRAFICA
for i in range(0,len(x_muestras),1):
    color_i = dsp._graf_lineacolor_i(i)
    estilo_i = ':' # estilo de la línea
    if i%2 == 0:  #es impar
        estilo_i = '--'
    # muestreo x[n]
    ki = x_muestras[i]['ki']
    xki = x_muestras[i]['xki']
    puntofmt = color_i + estilo_i
    plt.stem(ki*dt, xki, linefmt = puntofmt,
             label='x'+str(i)+'[n]')
    # muestreo x(t)
    ti = x_muestras[i]['ti']
    xti = x_muestras[i]['xti']
    plt.plot(ti,xti, '-', color = color_i,
             label='x'+str(i)+'(t)')
plt.axhline(0,color='black')
plt.xlabel('t')
plt.ylabel('amplitud')
plt.title('muestreo y alias')
plt.legend()
plt.show()

[ muestras ] [ alias ] [ ejemplo ] [ algoritmo ] [ ejercicio alias ]
..


3. Ejercicios con alias

Referencia: McClellan ejercicio 42 p128

Muestre que x2[n] es un alias de x1[n] y encuentre dos frecuencias que sean alias de 0.4π rad.

x1[n] = 7 cos(0.4π n – 0.2π)

x2[n] = 7 cos(8.4π n – 0.2π) = cos(0.4π n + 4(2π n) – 0.2π )  = cos(0.4π n- 0.2π)

Usando el algoritmo del numeral anterior se realiza la gráfica.

muestreo Alias coseno

Para los alias se usa la frecuencia normalizada en radianes: ωrad , considerando que si se suma varias veces 2π se obtiene un alias.

ωalias= 0.4π + 2π k   ; donde k = 0,1,2,3…

El alias principal se define como la frecuencia única en e intervalo entre
-π <ωalias < π , que para el ejercicio es 0.4π por lo que se propone usar como alias:

x3[n] = 7 cos((0.4π + 2π )n – 0.2π)

x4[n] = 7 cos((0.4π + 2(2π) )n – 0.2π)

muestreoAlias04_cosenoEl bloque de ingreso para el algoritmo del ejercicio es:

# INGRESO
# usar np.pi para evitar frecuencia fundamental automatica
x1 = 7*sym.cos(0.4*np.pi*n - 0.2*pi)
x2 = 7*sym.cos(8.4*np.pi*n - 0.2*pi)
x3 = 7*sym.cos((0.4*np.pi+2*pi*(1))*n - 0.2*pi )
x4 = 7*sym.cos((0.4*np.pi+2*pi*(2))*n - 0.2*pi)

##x_senales = [x1,x2]
##x_etiqueta = ['x1[n]','x2[n]']
x_senales = [x1,x3,x4]
x_etiqueta = ['x1[n]','x3[n]','x4[n]']

fs = 60  # muestreo discreto, freq_sampling
fs_veces = 40 # sobremuestreo para x(t)

muestrasN = 5 + 1  # para x[n]

[ muestras ] [ alias ] [ ejemplo ] [ algoritmo ] [ ejercicio alias ]

1.8 Señales Periódicas – frecuencia fundamental

[ periodo fundamental ] [ frecuencia fundamental ] [ señal a exponente ]
..


1. Periodo Fundamental de una señal periódica

Referencia: McClellan 3.3 p87

Una señal periódica repite  su forma de onda cada intervalo de tamaño T0. El intervalo T0 se conoce como periodo de x(t) y es la repetición más pequeña del intervalo, también llamada periodo fundamental.

x(t+T_0) = x(t)

Una señal periódica se puede formar como una suma de N+1 componentes con frecuencias armónicas que son múltiplos de la frecuencia Fundamental F0.

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

donde fk es la késima frecuencia del componente: fk = k F0

1.2 Ejercicio

Referencia: McClellan Ejemplo 3.8 p87

Determinar la frecuencia fundamental  de la señal x(t), siendo:

x(t) = \cos^2(4\pi t)

1.3 Resultados con el algoritmo

Si se plantea el ejercicio como dos señales, el coseno simple y luego el coseno cuadrado, se observa que la frecuencia fundamental del coseno es 2 Hz.
Sin embargo, la frecuencia del coseno cuadrado es 4Hz.

Algoritmo: Espectro – Operaciones en dominio de tiempo y frecuencia

x_senales: 
senal:   cos(4*pi*t)
  espectro:
   freq : [-2.  2.]
   ampl : [1/2 1/2]
   fase : [0 0]
   etiq : ['$\\frac{1}{2}$' '$\\frac{1}{2}$']
   x_expand : cos(4*pi*t)
   freq_max : 2.0
   freq_min : 2.0
   BW : 0.0
senal:   cos(4*pi*t)**2
  espectro:
   freq : [-4.  0.  4.]
   ampl : [1/4 1/2 1/4]
   fase : [0 0 0]
   etiq : ['$\\frac{1}{4}$' '$\\frac{1}{2}$' '$\\frac{1}{4}$']
   x_expand : cos(8*pi*t)/2 + 1/2
   freq_max : 4.0
   freq_min : 4.0
   BW : 0.0

la gráfica es:

espectroSenales Operación elevado Cuadradobloque de Ingreso:

# INGRESO
x1 = sym.cos(4*pi*t)
x2 = x1**2
# Para espectro de frecuencias y fasores
x_senales = [x1,x2]
x_etiqueta = ['x(t)','(x(t))^2']

[ periodo fundamental ] [ frecuencia fundamental ] [ señal a exponente ]
..


2. Frecuencia fundamental con máximo común divisor

La frecuencia fundamental es el entero F0 más grande, tal que fk = k F0. Se puede obtener con la instrucción de Numpy:

>>> import numpy as np
>>> freq = [12,20,60]
>>> np.gcd.reduce(freq)
4
>>>

de donde se interpreta que la frecuencia fundamental para 12,20 y 60 Hz se pueden crear a partir de 4 Hz.

12 Hz es la 3ra armónica,
20 Hz es la 5ta armónica,
60 Hz es la 15ta armónica.

Recuerde que la operación gcd se realiza con números enteros, en el caso de tener decimales, hay que multiplicar las frecuencias por 10, 100,1000 que conviertan las frecuencias a ser analizadas a números enteros.

[ periodo fundamental ] [ frecuencia fundamental ] [ señal a exponente ]

..


3. Ejercicio, seno al cubo

Referencia: McClellan 3.4.2 p91

Observe el resultado del ejercicio anterior sobre coseno al cuadrado.
Ahora desarrolle el ejercicio para seno al cubo y observe la frecuencia fundamental F0.

x(t) = \sin^3(4\pi t)

Elevar a una potencia una función periódica genera otra función periódica con el mismo o menor periodo.

La expansión en exponenciales usando Euler para el ejercicio muestra los componentes para el espectro de frecuencias:

x(t) = \frac{3 j}{8}(-j \sin(4\pi t) + \cos(4 \pi t)) + \frac{-3j}{8}(j \sin(4\pi t) + cos(4\pi t)) + \frac{-j}{8}(-j \sin(12\pi t) + \cos(12 \pi t)) + \frac{j}{8}(j \sin(12 \pi t) + \cos(12 \pi t))

Las frecuencias resultantes muestran que se tienen la 1ra y 3ra armónicas de la frecuencia fundamental de 2 Hz. Los coeficientes de la expresión en términos k de armónicas son:

ak cuando
0 k = 0
-/+ j 3/8 k = ±1
0 k = ±2
± j 1/8 k = ±3

que se pueden observar en el espectro de frecuencias.

x_senales: 
senal:   sin(4*pi*t)
  espectro:
   freq : [-2.  2.]
   ampl : [1/2 1/2]
   fase : [pi/2 -pi/2]
   etiq : ['$\\frac{1}{2}$$ e^j(\\frac{\\pi}{2})$'
 '$\\frac{1}{2}$$ e^j(- \\frac{\\pi}{2})$']
   x_expand : I*(-I*sin(4*pi*t) + cos(4*pi*t))/2 - I*(I*sin(4*pi*t) + cos(4*pi*t))/2
   freq_max : 2.0
   freq_min : 2.0
   BW : 0.0
senal:   sin(4*pi*t)**3
  espectro:
   freq : [-6. -2.  2.  6.]
   ampl : [1/8 3/8 3/8 1/8]
   fase : [-pi/2 pi/2 -pi/2 pi/2]
   etiq : ['$\\frac{1}{8}$$ e^j(- \\frac{\\pi}{2})$'
 '$\\frac{3}{8}$$ e^j(\\frac{\\pi}{2})$'
 '$\\frac{3}{8}$$ e^j(- \\frac{\\pi}{2})$'
 '$\\frac{1}{8}$$ e^j(\\frac{\\pi}{2})$']
   x_expand : 3*I*(-I*sin(4*pi*t) + cos(4*pi*t))/8 - 3*I*(I*sin(4*pi*t) + cos(4*pi*t))/8 - I*(-I*sin(12*pi*t) + cos(12*pi*t))/8 + I*(I*sin(12*pi*t) + cos(12*pi*t))/8
   freq_max : 6.0
   freq_min : 2.0
   BW : 4.0

se observa que la frecuencia fundamental en ambos casos es 2Hz

gráfica de espectro de frecuencias

espectro Senales Operación elevado Cubo

[ periodo fundamental ] [ frecuencia fundamental ] [ señal a exponente ]

1.7 Espectro – Operaciones en dominio de tiempo y frecuencia

[ Escala o sumar ‘c’ ] [ desplaza t ] [ derivada ] [ desplaza_f ] [ algoritmo ]

Operaciones en dominio de tiempo y frecuencia

Referencia: McClellan 3.3 p81

Uno de los beneficios de uso del espectro de frecuencias es que las operaciones sobre x(t) se muestran como resultados simples en la gráfica de espectro. Estos cambios relacionados en el dominio de t y f se denominan propiedades de la representación espectral. Para los ejemplos se usará el algoritmo de producto de sinusoides.

[ Escala o sumar ‘c’ ] [ desplaza t ] [ derivada ] [ desplaza_f ] [ algoritmo ]
..


1. Escalar o sumar con una constante

Referencia: McClellan 3.3.1 y 3.3.2 p81

Multiplicar una señal por un factor ‘c‘ constante cambia la escala de amplitud en el espectro por el mismo factor ‘c‘, pero deja las frecuencias sin cambios.

c x(t) = c \sum_{k=-M}^{M} a_k e^{j 2 \pi f_k t} = \sum_{k=-M}^{M} (c a_k) e^{j 2 \pi f_k t}

Sumar una constante a una señal x(t)+c, cambia la amplitud compleja de solo del componente constante o frecuencia f=0

x(t) + c = \sum_{f \ne 0} a_k e^{j 2 \pi f_k t} + \Big[a_0 e^{j 2 \pi (0) t} + c \Big]

1.1 Ejemplo

Para la señal x(t) mostrada

x(t) = \frac{3}{2} + 6 \cos ( 6 \pi t - \pi /3) + 4 \cos ( 14 \pi t + \pi/4)

realizar las operaciones

y(t) = 2x(t) +6

La amplitud de las señales en f=±3 y ±7 se duplican en el factor 2,  El valor de la constante a la frecuencia cero, también se duplica por el factor 2 y se añade la constante 6  dando como resultado 2(1.5) +6 = 9

1.2 Resultado con algoritmo

Ingreso

# INGRESO
x1 = 3/2 + 6*sym.cos(6*pi*t - pi/3) + 4*sym.cos(14*pi*t + pi/4)
x2 = 2*x1
x3 = 6
x4 = x2+x3

# Para espectro de frecuencias y fasores
x_senales = [x1,x2,x3,x4]
x_etiqueta = ['x1(t)','(2)x1(t)','x3(t)','(2)x1(t)+x3(t)']

Resultado

senal:   6*sin(6*pi*t + pi/6) + 4*cos(14*pi*t + pi/4) + 1.5
  espectro:
   freq : [-7. -3.  0.  3.  7.]
   ampl : [2.  3.  1.5 3.  2. ]
   fase : [-pi/4 pi/3 0 -pi/3 pi/4]
   etiq : ['$2$$ e^j(- \\frac{\\pi}{4})$' '$3$$ e^j(\\frac{\\pi}{3})$' '$1.5$'
 '$3$$ e^j(- \\frac{\\pi}{3})$' '$2$$ e^j(\\frac{\\pi}{4})$']
   x_expand : 6*sin(6*pi*t + pi/6) + 4*cos(14*pi*t + pi/4) + 1.5
   freq_max : 7.0
   freq_min : 3.0
   BW : 4.0
senal:   12*sin(6*pi*t + pi/6) + 8*cos(14*pi*t + pi/4) + 3.0
  espectro:
   freq : [-7. -3.  0.  3.  7.]
   ampl : [4 6 3 6 4]
   fase : [-pi/4 pi/3 0 -pi/3 pi/4]
   etiq : ['$4$$ e^j(- \\frac{\\pi}{4})$' '$6$$ e^j(\\frac{\\pi}{3})$' '$3$'
 '$6$$ e^j(- \\frac{\\pi}{3})$' '$4$$ e^j(\\frac{\\pi}{4})$']
   x_expand : 12*sin(6*pi*t + pi/6) + 8*cos(14*pi*t + pi/4) + 3
   freq_max : 7.0
   freq_min : 3.0
   BW : 4.0
senal:   6
  espectro:
   freq : [0.]
   ampl : [6]
   fase : [0]
   etiq : ['$6$']
   x_expand : 6
   freq_max : 0.0
   freq_min : 0
   BW : 0.0
senal:   12*sin(6*pi*t + pi/6) + 8*cos(14*pi*t + pi/4) + 9.0
  espectro:
   freq : [-7. -3.  0.  3.  7.]
   ampl : [4 6 9 6 4]
   fase : [-pi/4 pi/3 0 -pi/3 pi/4]
   etiq : ['$4$$ e^j(- \\frac{\\pi}{4})$' '$6$$ e^j(\\frac{\\pi}{3})$' '$9$'
 '$6$$ e^j(- \\frac{\\pi}{3})$' '$4$$ e^j(\\frac{\\pi}{4})$']
   x_expand : 12*sin(6*pi*t + pi/6) + 8*cos(14*pi*t + pi/4) + 9
   freq_max : 7.0
   freq_min : 3.0
   BW : 4.0

gráfica de espectro de frecuencias
espectro Senales Operación escala Suma Constante

[ Escala o sumar ‘c’ ] [ desplaza t ] [ derivada ] [ desplaza_f ] [ algoritmo ]
..


2. Desplazamiento en tiempo o multiplicar por exponencial complejo

Referencia: McClellan 3.3.3 p84

Para desplazamiento en el tiempo, las frecuencias no cambia, solo hay cambio de fase en las amplitudes complejas del espectro.

y(t) = x(t-\tau_d)\leftrightarrow b_k = a_k e^{-j2\pi f_k \tau_d}

dado que

y(t) = x(t-\tau_d) = \sum_k a_k e^{-j2\pi f_k (t-\tau_d)} = \sum_k \Big(a_k e^{-j2\pi f_k \tau_d }\Big) e^{j2\pi f_k t}

2.1 Ejemplo

Considere x(t) con un retraso τd que es 1/4 de su periodo.

x(t) = 6 \cos(250 \pi t) T = 1/125 \tau_d = \frac{T}{4} = \Big( \frac{1}{125} \Big) \frac{1}{4} = \frac{1}{500} = 0.002 e^{-j250 \pi (0.002)} = e^{-j 0.5 \pi} = -j x(t-0.002) = 6 \sin(250 \pi t)

2.2 Resultado con algoritmo

Ingreso

# INGRESO
x0 = 6*sym.cos(250*pi*t)
x1 = x0.subs(t,t-0.002)

# Para espectro de frecuencias y fasores
x_senales = [x0,x1]
x_etiqueta = ['x(t)','x(t-0.002)']

Resultado

x_senales: 
senal:   6*cos(250*pi*t)
  espectro:
   freq : [-125.  125.]
   ampl : [3 3]
   fase : [0 0]
   etiq : ['$3$' '$3$']
   x_expand : 6*cos(250*pi*t)
   freq_max : 125.0
   freq_min : 125.0
   BW : 0.0
senal:   6*cos(pi*(250*t - 0.5))
  espectro:
   freq : [-125.  125.]
   ampl : [3 3]
   fase : [0.5*pi -0.5*pi]
   etiq : ['$3$$ e^j(0.5 \\pi)$' '$3$$ e^j(- 0.5 \\pi)$']
   x_expand : 3*I*cos(250*pi*t) + 3*I*cos(250*pi*t - 1.0*pi) + 6*cos(250*pi*t - 0.5*pi)
   freq_max : 125.0
   freq_min : 125.0
   BW : 0.0

gráfica de espectro de frecuencias

espectro Señales Operación desplazamiento en tiempo

[ Escala o sumar ‘c’ ] [ desplaza t ] [ derivada ] [ desplaza_f ] [ algoritmo ]
..


3. Diferenciación en tiempo

Referencia: McClellan 3.3.4 p84

La derivada de una señal x(t) no realiza cambios en las frecuencias, las amplitudes en el espectro de frecuencia se modifican como

y(t) = \frac{d}{dt}x(t) \leftrightarrow b_k =( j 2 \pi f_k) a_k

dado que

y(t) = \frac{d}{dt} x(t) = \sum_k a_k \frac{d}{dt}e^{ j 2 \pi f_k t} = \sum_k \Big(( j 2 \pi f_k) a_k\Big) e^{j 2 \pi f_k t }

3.1 Ejemplo

Derivar x(t) que tiene sin() y una constante.

x(t) = 7 + 6 \cos(250 \pi t)

La frecuencia f=0 tiene magnitud 7 que es (j 2 π (0)) = 0 que se elimina de la expresión. En la frecuencia de 125 Hz el término se multiplica por (j 2 π (125))

(j 2 \pi (125)) (-3j) = 750 \pi \frac{d}{dt} x(t) = 1500 \pi \cos(250 \pi t)

3.2  Resultados con el algoritmo

Ingreso

# INGRESO
x0 = 7 + 6*sym.sin(250*pi*t)
x1 = sym.diff(x0,t)

# Para espectro de frecuencias y fasores
x_senales = [x0,x1]
x_etiqueta = ['x(t)',"x'(t)"]

Resultado

x_senales: 
senal:   6*cos(250*pi*t)
  espectro:
   freq : [-125.  125.]
   ampl : [3 3]
   fase : [0 0]
   etiq : ['$3$' '$3$']
   x_expand : 6*cos(250*pi*t)
   freq_max : 125.0
   freq_min : 125.0
   BW : 0.0
senal:   6*cos(pi*(250*t - 0.5))
  espectro:
   freq : [-125.  125.]
   ampl : [3 3]
   fase : [0.5*pi -0.5*pi]
   etiq : ['$3$$ e^j(0.5 \\pi)$' '$3$$ e^j(- 0.5 \\pi)$']
   x_expand : 3*I*cos(250*pi*t) + 3*I*cos(250*pi*t - 1.0*pi) + 6*cos(250*pi*t - 0.5*pi)
   freq_max : 125.0
   freq_min : 125.0
   BW : 0.0
>>> 

gráfica de espectro de frecuencias

espectro Señales Operación derivada

[ Escala o sumar ‘c’ ] [ desplaza t ] [ derivada ] [ desplaza_f ] [ algoritmo ]
..


4. Desplazamiento de frecuencia

Referencia: McClellan 3.3.5 p85

Multiplicar una señal x(t) por una sinusoide o una exponencial compleja tiene como resultado que las frecuencias en el espectro se desplazan.

y(t) = x(t) A \cos(2 \pi f_c t + \varphi) = x(t) \frac{1}{2} A e^{j\varphi}e^{j 2 \pi f_c t } + x(t) \frac{1}{2} A e^{-j\varphi}e^{-j 2 \pi f_c t }

que es el resultado observado en el producto de sinusoides y tema tratado para AM

y(t) = A e^{j\varphi}e^{j 2 \pi f_c t } \sum_k a_k e^{j 2 \pi f_k t } y(t) = \sum_k \Big(a_k A e^{j\varphi} \Big) e^{j 2 \pi (f_k+f_c) t }

Que se interpreta en el espectro de frecuencias, cada frecuencia en x(t) se desplaza por fc, suponiendo que fc>0. Todas las amplitudes complejas se multiplican por la amplitud compleja del exponencial complejo.

4.1 Ejemplo

Para una señal x(t) con varios componentes mostrada en la gráfica de espectro de frecuencias, añadir a la gráfica los resultados de las operaciones indicadas:

x_1(t) = x(t) \cos(2 \pi 9 t) x_1(t) = x(t) \sin(2 \pi 9 t)

Siendo

x(t) = 8 + 4 \cos(2 \pi t) + 4\cos(2\pi(2)t)+ 12\cos( 2 \pi (3)t + \pi/2)

4.2  Resultados con el algoritmo

Ingreso

# INGRESO
x0 = 8 + 4*sym.cos(DosPi*t) + 4*sym.cos(DosPi*2*t)+ 12*sym.cos(DosPi*3*t + pi/2)
x1 = x0*sym.cos(DosPi*9*t)
x2 = x0*sym.sin(DosPi*9*t)

# Para espectro de frecuencias y fasores
x_senales = [x0,x1,x2]
x_etiqueta = ['x(t)','x(t)*cos(2*pi*9*t)','x(t)*sin(2*pi*9*t)']

Resultado

x_senales: 
senal:   8 - 12*sin(3*t*(2*pi)) + 4*cos(t*(2*pi)) + 4*cos(2*t*(2*pi))
  espectro:
   freq : [-3. -2. -1.  0.  1.  2.  3.]
   ampl : [6 2 2 8 2 2 6]
   fase : [-pi/2 0 0 0 0 0 pi/2]
   etiq : ['$6$$ e^j(- \\frac{\\pi}{2})$' '$2$' '$2$' '$8$' '$2$' '$2$'
 '$6$$ e^j(\\frac{\\pi}{2})$']
   x_expand : -6*I*(-I*sin(6*pi*t) + cos(6*pi*t)) + 6*I*(I*sin(6*pi*t) + cos(6*pi*t)) + 4*cos(2*pi*t) + 4*cos(4*pi*t) + 8
   freq_max : 3.0
   freq_min : 1.0
   BW : 2.0
senal:   (8 - 12*sin(3*t*(2*pi)) + 4*cos(t*(2*pi)) + 4*cos(2*t*(2*pi)))*cos(9*t*(2*pi))
  espectro:
   freq : [-12. -11. -10.  -9.  -8.  -7.  -6.   6.   7.   8.   9.  10.  11.  12.]
   ampl : [3 1 1 4 1 1 3 3 1 1 4 1 1 3]
   fase : [-pi/2 0 0 0 0 0 pi/2 -pi/2 0 0 0 0 0 pi/2]
   etiq : ['$3$$ e^j(- \\frac{\\pi}{2})$' '$1$' '$1$' '$4$' '$1$' '$1$'
 '$3$$ e^j(\\frac{\\pi}{2})$' '$3$$ e^j(- \\frac{\\pi}{2})$' '$1$' '$1$'
 '$4$' '$1$' '$1$' '$3$$ e^j(\\frac{\\pi}{2})$']
   x_expand : 3*I*(-I*sin(12*pi*t) + cos(12*pi*t)) - 3*I*(I*sin(12*pi*t) + cos(12*pi*t)) - 3*I*(-I*sin(24*pi*t) + cos(24*pi*t)) + 3*I*(I*sin(24*pi*t) + cos(24*pi*t)) + 8*cos(18*pi*t)
   freq_max : 12.0
   freq_min : 6.0
   BW : 6.0
senal:   (8 - 12*sin(3*t*(2*pi)) + 4*cos(t*(2*pi)) + 4*cos(2*t*(2*pi)))*sin(9*t*(2*pi))
  espectro:
   freq : [-12. -11. -10.  -9.  -8.  -7.  -6.   6.   7.   8.   9.  10.  11.  12.]
   ampl : [3 1 1 4 1 1 3 3 1 1 4 1 1 3]
   fase : [0 pi/2 pi/2 pi/2 pi/2 pi/2 pi pi -pi/2 -pi/2 -pi/2 -pi/2 -pi/2 0]
   etiq : ['$3$' '$1$$ e^j(\\frac{\\pi}{2})$' '$1$$ e^j(\\frac{\\pi}{2})$'
 '$4$$ e^j(\\frac{\\pi}{2})$' '$1$$ e^j(\\frac{\\pi}{2})$'
 '$1$$ e^j(\\frac{\\pi}{2})$' '$3$$ e^j(\\pi)$' '$3$$ e^j(\\pi)$'
 '$1$$ e^j(- \\frac{\\pi}{2})$' '$1$$ e^j(- \\frac{\\pi}{2})$'
 '$4$$ e^j(- \\frac{\\pi}{2})$' '$1$$ e^j(- \\frac{\\pi}{2})$'
 '$1$$ e^j(- \\frac{\\pi}{2})$' '$3$']
   x_expand : I*(-I*sin(14*pi*t) + cos(14*pi*t)) - I*(I*sin(14*pi*t) + cos(14*pi*t)) + I*(-I*sin(16*pi*t) + cos(16*pi*t)) - I*(I*sin(16*pi*t) + cos(16*pi*t)) + 4*I*(-I*sin(18*pi*t) + cos(18*pi*t)) - 4*I*(I*sin(18*pi*t) + cos(18*pi*t)) + I*(-I*sin(20*pi*t) + cos(20*pi*t)) - I*(I*sin(20*pi*t) + cos(20*pi*t)) + I*(-I*sin(22*pi*t) + cos(22*pi*t)) - I*(I*sin(22*pi*t) + cos(22*pi*t)) - 6*cos(12*pi*t) + 6*cos(24*pi*t)
   freq_max : 12.0
   freq_min : 6.0
   BW : 6.0

gráfica de espectro de frecuencias

espectro Señales Operación desplaza freq
[ Escala o sumar ‘c’ ] [ desplaza t ] [ derivada ] [ desplaza_f ] [ algoritmo ]
..


5. Algoritmo en Python

La instrucciones en Python para cada ejercicio a observar en el dominio de la frecuencia requiere actualiar la sección de INGRESO que fué proporcinada en cada sección anterior. El algoritmo se presenta como ejemplo la escala o suma de una constante.

# ejemplo 3-3.1 p82 escala o suma con una constante
# 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 = 3/2 + 6*sym.cos(6*pi*t - pi/3) + 4*sym.cos(14*pi*t + pi/4)
x2 = 2*x1
x3 = 6
x4 = x2+x3

# Para espectro de frecuencias y fasores
x_senales = [x1,x2,x3,x4]
x_etiqueta = ['x1(t)','(2)x1(t)','x3(t)','(2)x1(t)+x3(t)']

# PROCEDIMIENTO
x_conteo = len(x_senales)
x_espectro = [] ; freq_max_graf = 1; freq_min_graf = 0
for unasenal in x_senales:
    # operaciones para espectro
    x_term = dsp.x_list_term_Add(unasenal)
    Xe_term = dsp.cos_to_euler_one_term(x_term)
    x_term_expand = dsp.euler_to_cos_list(Xe_term)
    Xe_term_spectr = dsp.cos_spectrum_list(x_term)
    
    # espectro de cada señal
    un_espectro = {}
    freq = Xe_term_spectr['freq']
    ampl = Xe_term_spectr['ampl']
    fase = Xe_term_spectr['fase']
    un_espectro['freq'] = freq
    un_espectro['ampl'] = ampl
    un_espectro['fase'] = fase
    un_espectro['etiq'] = Xe_term_spectr['etiq']
    un_espectro['x_expand'] = x_term_expand
    
    # ancho de banda
    freq_max = float(max(freq))
    if len(freq[freq>0])>0:
        freq_min = float(min(freq[freq>0]))
    else:
        freq_min = 0
    un_espectro['freq_max'] = freq_max
    un_espectro['freq_min'] = freq_min
    freq_centro = (freq_max+freq_min)/2
    un_espectro['BW'] = freq_max-freq_min

    x_espectro.append(un_espectro)
    # freq_max para grafica, eje unificado x_senales
    if freq_max>freq_max_graf:
        freq_max_graf = freq_max

# SALIDA
print('x_senales: ')
for i in range(x_conteo):
    print('senal:  ',x_senales[i])
    print('  espectro:')
    unespectro = x_espectro[i]
    for unparam in unespectro:
        print('  ',unparam,':',unespectro[unparam])

# GRAFICAS de espectro de frecuencias ---------
graf_dx = 0.14
fig_espectro = plt.figure()
for i in range(0,x_conteo,1):
    unespectro = x_espectro[i]
    freq = unespectro['freq']
    ampl = unespectro['ampl']
    etiqueta = unespectro['etiq']
    ampl_max = float(max(ampl))
    # grafica
    graf_sub = x_conteo*100+10+i+1
    graf_fasor = fig_espectro.add_subplot(graf_sub)
    if freq_max_graf!=0:
        graf_fasor.set_xlim([-freq_max_graf*(1+graf_dx),
                             freq_max_graf*(1+graf_dx)])
    else:
        graf_fasor.set_xlim([-1,1])
    graf_fasor.set_ylim([0,ampl_max*(1+graf_dx*2)])
    graf_fasor.grid()
    graf_fasor.axhline(0,color='black')
    graf_fasor.axvline(0,linestyle='dotted',color='grey')
    graf_fasor.stem(freq,ampl) # grafica magnitud
    for k in range(0,len(freq),1): # etiquetas de fasor
        plt.annotate(etiqueta[k],xy=(freq[k],ampl[k]),
                     xytext=(0,5),textcoords='offset points',
                     ha='center')
    graf_fasor.set_ylabel(x_etiqueta[i])
    if i == 0:
        graf_fasor.set_title('Espectro: x_senales')
    if i ==(x_conteo-1):
        graf_fasor.set_xlabel('freq Hz')

plt.show()

[ Escala o sumar ‘c’ ] [ desplaza t ] [ derivada ] [ desplaza_f ] [ algoritmo ]