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 ±ω₀ con amplitudes de (1/2 A e^±ω₀)

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 ±ω₀/f_s, con todos los alias en
±ω₀/f_s + 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 f_s=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)

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 [-f_s/2, +f_s/2].

En resumen, en sobre-muestreo la frecuencia original f₀ es menor que f_s/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 f_s < 2f₀ 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.

las instrucciones para el bloque de ingreso del algoritmo son las mismas que en el ejercicio anterior, excepto que f_s = 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 [f₀,2f₀] se presenta aliasing denominado «folding» o plegado de señal. Para el ejercicio de prueba y con f_s = 125, se tiene:

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

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(f_s/2π) = f_s/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 f_s=100Hz, el resultado es una constante, pues siempre se obtiene la muestra de la señal con la misma magnitud.

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 f_s 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 ]