Autor: Edison Del Rosario

Pares DFT

Referencia: McClellan Tabla 8.1 p327-328

x[n]	X[k]
δ[n]	1
δ[n-nd]	$e^{-j(2\pi k/N)n_d}$
$r_L[n] = \mu [n] - \mu [n-L]$	$\frac{\sin\Big(\frac{1}{2}L(2\pi k/N)\Big)}{\sin\Big(\frac{1}{2}(2\pi k/N)\Big)}e^{-j(2\pi k/N)(L-1)/2}$
	$D_L(2\pi k/N) = \frac{\sin\Big(\frac{1}{2}L(2\pi k/N)\Big)}{\sin\Big(\frac{1}{2}(2\pi k/N)\Big)}$
$r_L[n] e^{j(2\pi k_0/N)n}$	$D_L(2 \pi (k-k0)/N)e^{-j(2\pi (k-k_0)/N)(L-1)/2}$

DFT Propiedades

 

Propiedad	dominio tiempo x[n]	dominio frecuencia X[k]
Periódica	x[n] = x[n+N]	X[k] = X[k+N]
Linealidad	ax1[n] +bx2[n]	aX1[k] +bX2[k]
Simetría Conjugada	x[n] Real	X[N-k] = x*[k]
Conjugación	x*[n]	X*[N-K]
Reversible en tiempo	x[ ((N-n))N ]	X[N-k]
Retraso	x[ ((n-nd))N ]	$e^{ -j (2\pi k/N)n_d} X[k]$
Desplazamiento en frecuencia	$x[n] e^{ j (2\pi k_0/N)n}$	X[k-k0]
Modulación	$x[n] \cos\Big((2π k_0/N)n\Big)$	$\frac{1}{2}X[k-k_0] + \frac{1}{2}X[k+k_0]$
Convolución	$\sum_{m=0}^{N-1} h[m]x[((n-m))_N ]$	H[k]X[k]
Teorema de Parseval	$\sum_{n=0}^{N-1} |x[n]|^2 = \frac{1}{N}\sum_{k=0}^{N-1} |X[k]|^2$

IDFT: [ periodicidad ] [ convolución ]
..

---

1. Periódicidad

Referencia: McClellan 8.3.1 p319

La sumatoria de IDFT también tiene propiedad e periodicidad. Considera evaluar una función en n+N, con n en el intervalo 0≤n≤N-1:

$x[n+N] = \frac{1}{N}\sum_{k=0}^{N-1} X[k] e^{j(2\pi /N)k(n+N)}$ $= \frac{1}{N}\sum_{k=0}^{N-1} X[k] e^{j(2\pi /N)k(n)} \cancel{e^{j(2\pi /N)k(N)} }$

un factor exponencial se convierte en 1, quedando

$= \frac{1}{N}\sum_{k=0}^{N-1} X[k] e^{j(2\pi /N)k(n)} = x[n]$

1. 1 Ejemplo: IDFT de impulso desplazado

Para una señal de L=N=10 puntos x[n]=δ[n-14] tendrá una DFT:

$X[k] = e^{-j(0.2\pi (14)k}$

El resultado de la IDFT tendrá un valor diferente de cero en n=4.

con el algoritmo para no.fft.ifft() se pueden obtener los valores para

N: 10
X[k]:
[ 1.        +0.00000000e+00j -0.80901699-5.87785252e-01j
  0.30901699+9.51056516e-01j  0.30901699-9.51056516e-01j
 -0.80901699+5.87785252e-01j  1.        +1.71450552e-15j
 -0.80901699-5.87785252e-01j  0.30901699+9.51056516e-01j
  0.30901699-9.51056516e-01j -0.80901699+5.87785252e-01j]
Re(x[n]):
[-4.88498131e-16 -7.98670570e-16 -9.76996262e-16 -9.51130884e-16
  1.00000000e+00  9.49240686e-16  9.76996262e-16  8.05446317e-16
  5.32907052e-16 -4.88554849e-18]
Im(x[n]):
[-7.16727868e-16 -3.76270034e-16  2.39540920e-17  4.01492493e-16
  2.23107982e-15  3.94763191e-16  4.54309871e-17 -3.94454259e-16
 -7.26484270e-16 -8.82784150e-16]

considerando como bloque de ingreso:

# INGRESO
X_k = lambda k: np.exp(-1j*0.2*np.pi*14*k)
N  = 10 # Tramos en intervalo
k0 = 0 # inicio de gráfica

IDFT: [ periodicidad ] [ convolución ]
..

---

2.  Propiedad de la Convolución

La convolución en el domino del tiempo se puede desarrollar como una multiplicación en el dominio de la frecuencia.

$y[n] = h[n] * x[n] \leftrightarrow Y(e^{j \omega} )= H(e^{j \omega}) X(e^{j \omega})$

Siendo la DFT una versión discreta  de tamaño finito, de la DTFT contínua:

$y[n] = \sum_{m=0}^{M-1} h[m] x[n-m]$

donde se supone que h[n]= 0 fuera del intervalo n= 0, 1, … , M-1,
x[n]=0 fuera del intervalo n= 0, 1, … , L-1
h[n] es una secuencia d M puntos y x[n] es una secuencia de L puntos.

$Y[k] = H[k] X[k]$

obteniendo y[n]

$y[n] = \frac{1}{N}\sum_{k=0}^{N-1} Y[k] e^{j(2 \pi /N)kn}$

2. 1 Ejemplo: Convolución de pulsos

Suponga una señal h[n] compuesta de 16 puntos de los que 6 se encuentran en valor 1 y el resto en cero. La señal x[n] es una señal de pulsos de tamaño 16 puntos, 10 en estado uno, el resto en cero. Considerando que Y[k] = H[k] X[k], use np.fft.fft y np.fft.ifft para mostrar la convolución resultante en el dominio del tiempo.

k0x: 0 ; i: 8
freq_Hz: 0.0 ; Xk[k0x]: (10+0j)
[ i, freq[ki],|X[k]| ,X[k]]
[0, 0.0, 10.0, (10+0j)]
[1, 0.2, 4.735650251450756, (-0.9238795325112867-4.6446560597607585j)]
[2, 0.4, 1.8477590650225735, (1.7071067811865475-0.7071067811865476j)]
[3, 0.6000000000000001, 0.688811980233627, (-0.38268343236508984-0.5727262301542022j)]

k0h: 0 ; i: 8
freq_Hz: 0.0 ; Hk[k0h]: (6+0j)
[ i, freq[ki],|H[k]| ,H[k]]
[0, 0.0, 6.0, (6+0j)]
[1, 0.2, 4.735650251450756, (2.630986313697834-3.937549278574211j)]
[2, 0.4, 1.8477590650225735, (-0.7071067811865476-1.7071067811865475j)]
[3, 0.6000000000000001, 0.688811980233627, (0.6755766511785423+0.1343805510323453j)]

Algoritmo en Python

# ejemplo 8.3.3 p324 DFT convolucion
# telg1034 DSP fiec-espol edelros@espol.edu.ec
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
L  = 16   # Tramos en un periodo

# señal es rectangular periódico rect[n]
f0 = 0.2  # frecuencia fundamental de señal
ancho = (10/16)*(1/f0) # ancho del impulso rectangular
xt = lambda t: np.heaviside(t,1)-np.heaviside(t-ancho,1)

# señal H es rectangular periódico rect[n]
f0 = 0.2  # frecuencia fundamental de señal
anchoh = (6/16)*(1/f0) # ancho del impulso rectangular
ht = lambda t: np.heaviside(t,1)-np.heaviside(t-anchoh,1)

# para DFT
N = 1*L # N par mayor que L, tramos DFT

# PROCEDIMIENTO
# muestreo de señal x[n]
t0 = 0 ; tn = 1/f0 # intervalo [t0,tn]
dt = (tn-t0)/L
ti = np.arange(t0,tn,dt)

xn = xt(ti)
hn = ht(ti)

# verificar que L<=N
if L<=N:
    # DFT transformada rapida de Fourier
    Xk = np.fft.fft(xn,n=N)
    Xkmagn = np.abs(Xk) # para gráfica
    Xkfase = np.angle(Xk)
    k0x = np.argmax(Xk) # Pico |X[k]|
    
    # DFT transformada rapida de Fourier
    Hk = np.fft.fft(hn,n=N)
    Hkmagn = np.abs(Hk) # para gráfica
    Hkfase = np.angle(Hk)
    k0h = np.argmax(Hk) # Pico |X[k]|

    freq_Hz = np.fft.fftfreq(N,dt) # frecuencias en Hz
    mitad = int(N/2)
    

    # Convolucion en dominio de frecuencia
    Yk = Hk*Xk
    yn = np.fft.ifft(Yk,n=N)
    yn_re = np.real(yn)
    yn_im = np.imag(yn)

# SALIDA
if N>=L:
    print('k0x:',k0x,'; i:',mitad+k0x)
    print('freq_Hz:',freq_Hz[k0x],
          '; Xk[k0x]:',Xk[k0x])
    print('[ i, freq[ki],|X[k]| ,X[k]]')
    kmin = k0x-4 ; kmax = k0x+4
    if kmin<0:
        kmin = 0
    if kmax>N:
        kmax = N
    for i in range(kmin,kmax,1):
        print([i,freq_Hz[i],Xkmagn[i],Xk[i]])
    print()
    print('k0h:',k0h,'; i:',mitad+k0h)
    print('freq_Hz:',freq_Hz[k0h],
          '; Hk[k0h]:',Hk[k0h])
    print('[ i, freq[ki],|H[k]| ,H[k]]')
    kmin = k0h-4 ; kmax = k0h+4
    if kmin<0:
        kmin = 0
    if kmax>N:
        kmax = N
    for i in range(kmin,kmax,1):
        print([i,freq_Hz[i],Hkmagn[i],Hk[i]])
        
else:
    print('Revisar: L debe ser menor o igual a N')
    print('L:',L,'; N:',N)

# GRAFICA
if N>=L:
    LN_texto= ', L='+str(L)+', N='+str(N)
    #zp_texto = 'zero padding:'+str(N-L)
    plt.subplot(311)
    plt.stem(np.real(xn),label='Re(x[n])',
             markerfmt ='C2.',linefmt='C2-')
    plt.stem(L,0,markerfmt ='C3.')
    plt.xlabel('ti (s)')
    
    plt.title('Convolución entre x[n],H[n]'+LN_texto)
    plt.legend()

    plt.subplot(312)
    plt.stem(np.real(hn),label='Re(h[n])',
             markerfmt ='C1.',linefmt='C1-')
    plt.stem(L,0,markerfmt ='C1.')
    plt.xlabel('ti (s)')
    plt.legend()

    plt.subplot(313)
    plt.stem(np.real(yn),label='Re(y[n])',
             markerfmt ='C0.',linefmt='C0-')
    plt.stem(L,0,markerfmt ='C0.')
    plt.xlabel('ti (s)')
    plt.legend()
    plt.tight_layout()
    plt.show()

 

DFT: [ periódicas ] [ freq negativas ] [ simetría conjugada]  [ zero_padding ] [ ejemplo_cos ]
..

---

1. Periódicidad

Referencia: McClellan 8.2.1 p311

La DFT de una secuencia finita de N-puntos es una versión discreta de DTFT:

$X[k] = \sum_{n=0}^{N-1} x[n] e^{-j(2\pi k/N)}$ $= X\Big( e^{j\omega} \Big) \Big|_{\omega=(2\pi k/N)}$ $= X\Big( e^{j (2\pi k/N)}\Big)$

k= 0, 1, …, N-1

Dado que DTFT X(e^jω) es siempre periódica con periodo 2π, la DFT X[k] debe ser por lo tanto también periódica. Esto implica que el índice k se mantiene siempre entre 0 y N-1, sin afectar ser evaluado para k≥N o k<0.

Considere que ω + 2π = 2πk/N + 2π(N/N) = 2π(k+N)/N = ω_k+N
que implica que los coeficientes X[k] tienen periodo igual a N. Por lo que no se considera necesario más cálculos X[k] fuera del intervalo 0≤k<N.

X[k] = X[k+N]

DFT: [ periódicas ] [ freq negativas ] [ simetría conjugada]  [ zero_padding ] [ ejemplo_cos ]
..

---

2. Frecuencias Negativas y DFT

Referencia: McClellan 8.2.2 p311, 8.2.3.1 p313

Todos los resultados de DFT y IDFT tienen índices no negativos por ser convenientes para los cálculos y expresiones matemáticas. Para la gráfica de espectro de frecuencias se espera utilizar frecuencias positivas y negativas, con simetría en las conjugadas cuando la señal es tipo real.

Cuando la DFT de tamaño N, siendo N par, la transformada X[k] tiene un valor a k=N/2. Cuando N es impar, lo indicado no aplica. El índice k=N/Se  corresponde a una frecuencia normalizada ω=2π(N/2)/N)=π. Sin embargo los valores en el espectro son de periodos 2π, siendo ω=-π un alias de ω=π.

Si se realiza la gráfica de forma centrada, la frecuencia esta entre [-π, π]

En Python el desplazamiento de los índices se puede realizar usando np.fft.fftshift() que ubica el centro en la muestra N/2, siempre que N sea par.

DFT: [ periódicas ] [ freq negativas ] [ simetría conjugada]  [ zero_padding ] [ ejemplo_cos ]
..

---

3. Simetría conjugada de DFT

Para una señal x[n] de tipo real, existe la simetría conjugada en la DTFT, los coeficientes deben satisfacer la propiedad:

X[-1]  = X^*[1]
X[-2]  = X^*[2]
…
X[-k]  = X^*[k]

La DFT de una señal real satisface:

X[N-k] = X^*[k]
k = 0, 1, …, N-1

DFT: [ periódicas ] [ freq negativas ] [ simetría conjugada]  [ zero_padding ] [ ejemplo_cos ]
..

---

4. Relleno de ceros – Zero Padding

Referencia: McClellan 8.2.4 p314

Para suavizar una gráfica de una DTFT como respuesta de frecuencia, se necesita muestras de frecuencia con tamaños de paso pequeños. De forma implícita, la DFT asume que la longitud o tamaño de la transformada es igual al tamaño  o longitud de la señal L, por lo que la DFT de L puntos calcula las muestras de DTFT para (ω=2π/L)k, siendo k=0,1,…,L-1. Las L muestras de frecuencia son con espacios iguales en el intervalo 0≤ω≤2π, con resultado:

$H[k] = H \Big( e^{j(2\pi /L)k} \Big) = H \Big( e^{j\omega}\Big) \Big|_{\omega=(2\pi/L)/k}$

k = 0, 1, …, L-1

Para disponer de  muestras de frecuencia a espacios menores que (2π/N) donde N>L, una forma sencilla es rellenar con ceros o zero-padding. Con lo que la señal se alarga antes de realizar el cálculo para la DFT de N puntos.

$h_{zp}[n] = \begin{cases} h[n] & n=\text {0, 1, ..., L-1} \\0 & n= L, L+1, ... , N-1 \end{cases}$

4.1 Ejemplo: Respuesta de frecuencias con DFT con zero-padding

Referencia: McClellan 8.4 p315

Suponer que requiere una gráfica de respuesta de frecuencia de un filtro FIR de 8 puntos. Para evaluar la respuesta de frecuencia con mejor resolución mediante tamaños de paso mas pequeños entre frecuencias, se usará relleno con ceros o zero-padding, usando un número par en N, por ejemplo N=80.

En el caso que N=40, el resultado tiene menor resolución considerando que el intervalo ω se mantiene entre [-π, π]

4.2 Algoritmo en Python

El filtro FIR se puede realizar usando un rectángulo con funciones escalón o Heaviside(n), luego evaluando para las muestras ti.

Para obtener la DFT se usa np.fft.fft(xn, n=N), si n es menor que el tamaño de la señal de entrada, la entrada se trunca. Si  es mas grande, se usa zero-padding para completar las muestras. Por lo que la instrucción de Numpy contempla el uso de zero-padding, lo que simplifica el uso del algoritmo presentado de forma didáctica paso a paso.

El centrado del espectro de frecuencias se puede realizar usando un desplazamiento de índices ki_shift.

El resultado del algoritmo, con X[k] centrada, se presenta como:

k0: 0 ; i: 40
freq_Hz: 0.0 ; Xk[k0]: (8+0j)
[ i, freq[ki],|X[k]| ,X[k]]
[0, 0.0, 8.0, (8+0j)]
[1, 0.0125, 7.871076020105053, (7.575558332074842-2.1365284158195137j)]
[2, 0.025, 7.491613901992368, (6.387650908672175-3.914357511197061j)]
[3, 0.037500000000000006, 6.883060301425795, (4.672226464156024-5.054386113140136j)]
...

Instrucciones en Python

# ejemplo 8.4 p315 DFT zero-padding
# telg1034 DSP fiec-espol edelros@espol.edu.ec
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
L = 16   # Tramos en un periodo

# señal es rectangular periódico rect[n]
f0 = 1/16  # frecuencia fundamental de señal
ancho = 0.5*(1/f0) # ancho del impulso rectangular
xt = lambda t: np.heaviside(t,1)-np.heaviside(t-ancho,1)

# para DFT
N = 80 #2*L # N par mayor que L, tramos DFT

# PROCEDIMIENTO
# muestreo de señal x[n]
t0 = 0 ; tn = 1*(1/f0) # intervalo [t0,tn]
dt = (tn-t0)/L
ti = np.arange(t0,tn,dt)
xn = xt(ti)

# verificar que L<=N
if L<=N:
    # DFT transformada rapida de Fourier
    Xk = np.fft.fft(xn,n=N)
    Xkmagn = np.abs(Xk) # para gráfica
    Xkfase = np.angle(Xk)
    freq_Hz = np.fft.fftfreq(N,dt) # frecuencias en Hz

    mitad = int(N/2)
    k0 = np.argmax(Xk) # Pico |X[k]|

    # ki indices de orden
    ki = np.arange(0,N,1)
    ki_shift = np.copy(ki)
    ki_shift[mitad:] = ki[mitad:]-N

# SALIDA
if N>=L:
    print('k0:',k0,'; i:',mitad+k0)
    print('freq_Hz:',freq_Hz[k0],
          '; Xk[k0]:',Xk[k0])
    print('[ i, freq[ki],|X[k]| ,X[k]]')
    kmin = k0-4 ; kmax = k0+4
    if kmin<0:
        kmin = 0
    if kmax>N:
        kmax = N
    for i in range(kmin,kmax,1):
        print([i,freq_Hz[i],Xkmagn[i],Xk[i]])
else:
    print('Revisar: L debe ser menor o igual a N')
    print('L:',L,'; N:',N)

# GRAFICA
if N>=L:
    LN_texto= ', L='+str(L)+', N='+str(N)
    zp_texto = 'zero padding:'+str(N-L)
    plt.subplot(311)
    plt.stem(ki[0:len(xn)],np.real(xn),label='x[n]',
             markerfmt ='C2.',linefmt='C2-')
    if N>L:
        kj = np.arange(L,N,dt)
        zp = np.zeros(N-L)
        plt.stem(kj,zp,label = 'zero padding',
                 markerfmt ='C3.')
    plt.xlabel('ki')
    
    plt.title('x[n] vs |x[k]|'+LN_texto+', '+zp_texto)
    plt.legend()

    plt.subplot(312)
    plt.stem(ki[0:mitad], Xkmagn[0:mitad],
             label='X[k], k_=[0:N/2]',
             markerfmt ='C0.', linefmt='C0-')
    plt.stem(ki[mitad:], Xkmagn[mitad:],
             label='X[k], k_=[N/2:]',
             markerfmt ='C1.',linefmt='C1-')
    plt.stem(N,0,markerfmt ='C3.')
    plt.xlabel('ki_')
    plt.legend()

    plt.subplot(313)
    plt.stem(ki_shift[0:mitad], Xkmagn[0:mitad],
             label='H[k] centrada]',
             markerfmt ='C0.',linefmt='C0-')
    plt.stem(ki_shift[mitad:], Xkmagn[mitad:],
             markerfmt ='C1.',linefmt='C1-')
    plt.stem(mitad,0,markerfmt ='C3.')
    plt.xlabel('ki')
    plt.legend()
    plt.legend()
    plt.tight_layout()
    plt.show()

DFT: [ periódicas ] [ freq negativas ] [ simetría conjugada]  [ zero_padding ] [ ejemplo_cos ]
..

---

5. DFT de una señal Real tipo Coseno

Referencia: McClellan 8.2.5 p316

Considere la señal de entrada x[n] coseno, de frecuencia 0.2 de forma discreta para L=50 tramos en tiempo en un intervalo de dos periodos. Para la DFT considere el caso que N=L, no se requiere relleno de ceros.

$x[n] = 5 cos(2\pi(0.1)n]$

n = 0, 1, …, L-1

Una señal senoidal se puede escribir como la suma de exponenciales complejos y usar la propiedad de linealidad para encontrar el resultado.
El resultado de la DFT muestra la característica de simetría donde ocurre un máximo (pico), en este caso en k=k₀=1 y k=N-k₀=50-1=49. Las correspondientes frecuencias en Hertz son -0.1 y 0.1. La altura de los picos son aproximadamente AL/2 = 5(50)/2 = 125.

Resultados numéricos con el algoritmo:

k0: 1 ; i: 26
freq_Hz: 0.1 ; Xk[k0]: (125+1.509903313490213e-14j)
[ i, freq[ki],|X[k]| ,X[k]]
[0, 0.0, 9.492406860545088e-15, (9.492406860545088e-15+0j)]
[1, 0.1, 125.0, (125+1.509903313490213e-14j)]
[2, 0.2, 4.842919552269065e-15, (-4.774353709997803e-15-8.120446056592065e-16j)]
[3, 0.30000000000000004, 3.4364294309926203e-15, (2.3207649852950894e-15-2.5343829855056534e-15j)]
[4, 0.4, 2.032850942959038e-15, (8.108973628378316e-16+1.8641159897474457e-15j)]

Para el caso que L<N, el valor de N se aumenta a 8*L, N=400, el resultado del pico de frecuencia se mantiene respecto al caso anterior.

los resulttados numéricos con el algoritmo son:

k0: 8 ; i: 208
freq_Hz: 0.1 ; Xk[k0]: (125.00000000000001+1.554312234475219e-14j)
[ i, freq[ki],|X[k]| ,X[k]]
[4, 0.05, 53.33891994969286, (4.999999999999986+53.10405240092081j)]
[5, 0.0625, 75.60689606012141, (33.14721982316793+67.95339984018891j)]
[6, 0.07500000000000001, 96.58278363275512, (70.7485680957908+65.74856809579083j)]
[7, 0.08750000000000001, 113.74232040562568, (105.80156476944595+41.75337523593449j)]
[8, 0.1, 125.00000000000001, (125.00000000000001+1.554312234475219e-14j)]
[9, 0.1125, 128.94638956201757, (119.85004661115146-47.572446945521044j)]
[10, 0.125, 125.00127087574828, (90.85388423989122-85.85388423989122j)]
[11, 0.1375, 113.46542771722382, (47.653106752530654-102.97370879930693j)]

Instrucciones en Python

# ejemplo 8.7 p318 DFT zero-padding cos()
# telg1034 DSP fiec-espol edelros@espol.edu.ec
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
L = 50   # Tramos en un periodo

# señal es rectangular periódico rect[n]
f0 = 0.1  # frecuencia fundamental de señal
xt = lambda t: 5*np.cos(2*np.pi*f0*t)

# para DFT
N = 8*L # N par mayor que L, tramos DFT

# PROCEDIMIENTO
# muestreo de señal x[n]
t0 = 0 ; tn = 1*(1/f0) # intervalo [t0,tn]
dt = (tn-t0)/L
ti = np.arange(t0,tn,dt)
xn = xt(ti)

# verificar que L<=N
if L<=N:
    # DFT transformada rapida de Fourier
    Xk = np.fft.fft(xn,n=N)
    Xkmagn = np.abs(Xk) # para gráfica
    Xkfase = np.angle(Xk)
    freq_Hz = np.fft.fftfreq(N,dt) # frecuencias en Hz

    mitad = int(N/2)
    k0 = np.argmax(Xk) # Pico |X[k]|

# SALIDA
if N>=L:
    print('k0:',k0,'; i:',mitad+k0)
    print('freq_Hz:',freq_Hz[k0],
          '; Xk[k0]:',Xk[k0])
    print('[ i, freq[ki],|X[k]| ,X[k]]')
    kmin = k0-4 ; kmax = k0+4
    if kmin<0:
        kmin = 0
    if kmax>N:
        kmax = N
    for i in range(kmin,kmax,1):
        print([i,freq_Hz[i],Xkmagn[i],Xk[i]])
else:
    print('Revisar: L debe ser menor o igual a N')
    print('L:',L,'; N:',N)

# GRAFICA
if N>=L:
    LN_texto= ', L='+str(L)+', N='+str(N)
    zp_texto = 'zero padding:'+str(N-L)
    plt.subplot(211)
    plt.stem(ti,np.real(xn),label='Re(x[n])',
             markerfmt ='C2.',linefmt='C2-')
    plt.stem(ti,np.imag(xn),label='Im(x[n])',
             markerfmt ='C4.',linefmt='C4:')
    if N>L:
        tj = np.arange(L*dt,N*dt,dt)
        zp = np.zeros(N-L)
        plt.stem(tj,zp,label = 'zero padding',
                 markerfmt ='C3.')
    plt.stem(L*dt,0,markerfmt ='C4.')
    plt.xlabel('ti (s)')
    
    plt.title('x[n] vs |x[k]|'+LN_texto+', '+zp_texto)
    plt.legend()

    plt.subplot(212)
    plt.stem(freq_Hz[0:mitad], Xkmagn[0:mitad],
             label='|X[k]| centrada]',
             markerfmt ='C0.',linefmt='C0-')
    plt.stem(freq_Hz[mitad:], Xkmagn[mitad:],
             markerfmt ='C1.',linefmt='C1-')
    plt.xlabel('freq (Hz)')
    plt.legend()
    plt.legend()
    plt.tight_layout()
    plt.show()

DFT: [ periódicas ] [ freq negativas ] [ simetría conjugada]  [ zero_padding ] [ ejemplo_cos ]

DFT: [ impulso unitario ] [ pulso rectangular ] [ exponencial complejo ] [ algoritmo ]

..

---

La DFT es una versión discreta del espectro de frecuencias de  tipo X(e^jω ) que es contínua y conocida como DTFT. Se puede construir pares de DFT haciendo la sustitución ω=(2π/N)/k en un par DTFT

1. DFT de un impulso unitario desplazado

Referencia: McClellan 8.1.2.1 p306. Transformada de Fourier – Señales aperiódicas continuas

Tomando la DFT de x[n] =δ[n], la sumatoria de DFT se simplifica a un término

$X[k] = \sum_{n=0}^{N-1} \delta [n] e^{-j(2\pi k/N)n}$ $= \delta [0] e^{-j(2\pi k/N)(0)} + \sum_{n=1}^{N-1} \cancel{\delta [n]} e^{-j(2\pi k/N)n}$ $= \delta [0] (1) + 0 = 1$

Usando la respuesta anterior, para x[n] = δ[n-n_d]m se puede escribir

$X[k] = \delta [n_d] e^{-j(2\pi k/N)n_d} = e^{-j\omega n_d}$ $\delta [n-n_d] \leftrightarrow e^{-j\omega n_d}$

Se crea a partir de una función rectangular con ancho igual al tamaño de paso.

# INGRESO
L  = 50   # Tramos en un periodo

# señal es impulso unitario periódico d[n], rectangular ancho "pequeño"
f0 = 0.2  # frecuencia fundamental de señal
ancho = 1/f0/L # ancho del impulso rectangular
n0 = 0.5 # desplazamiento en tiempo
xt = lambda n: np.heaviside(n-n0,1)-np.heaviside(n-n0-ancho,1)

# para DFT
N = 1*L # N par mayor que L, tramos DFT, FFT

DFT: [ impulso unitario ] [ pulso rectangular ] [ exponencial complejo ] [ algoritmo ]
..

---

2. DFT de un pulso rectangular

Referencia: Transformada de Fourier – Señales aperiódicas continuas

Para un pulso rectangular en el domino de tiempo discreto, expresado como:

$r [n] = \begin {cases} 1 & 0\leq n \leq L-1 \\ 0 & L \leq n \leq N-1 \end {cases}$

La DFT tiene un término de la forma

$R_L[k] = R_L \Big( e^{j\omega} \Big)$ $= \frac{\sin \Big( \frac{1}{2} L(2\pi k/N)\Big)}{\sin \Big( \frac{1}{2} (2\pi k/N)\Big)} e^{-j(2 \pi k/N)(L-1)/2}$

Se crea a partir de una función rectangular con ancho igual a la mitad de un periodo.

# INGRESO
L  = 50   # Tramos en un periodo

# señal es rectangular periódico rect[n]
f0 = 0.2  # frecuencia fundamental de señal
ancho = 0.5*(1/f0) # ancho del impulso rectangular
xt = lambda t: np.heaviside(t,1)-np.heaviside(t-ancho,1)

# para DFT
N = 1*L # N par mayor que L, tramos DFT

DFT: [ impulso unitario ] [ pulso rectangular ] [ exponencial complejo ] [ algoritmo ]
..

---

3. DFT de un exponencial complejo

Una señal exponencial compleja de tamaño finito se expresa como

$x [n] = r_L[n] e^{j \omega_0 n}= \begin {cases} e^{j\omega_0 n} & 0\leq n \leq L-1 \\ 0 & L \leq n \leq N-1 \end {cases}$

La DFT tiene un término de la forma

$X[k] = R_L \Big( e^{j\omega-\omega_0} \Big)$ $= \frac{\sin \Big( \frac{1}{2} L(2\pi k/N - \omega_0)\Big)}{\sin \Big( \frac{1}{2} (2\pi k/N- \omega_0)\Big)} e^{-j(2 \pi k/N- \omega_0)(L-1)/2}$

expresión que también se representa en términos de Dirichlet

$X[k] = D_L(2\pi k/N - \omega_0) e^{-j(2 \pi k/N- \omega_0)(L-1)/2}$

Dado que el exponencial solo contribuye a la fase, el resultado muestra que la magnitud depende solo del termino D_L

La DFT es muy simple cuando la frecuencia de la señal exponencial es un múltiplo entero de 2π/N y la cantidad de muestras DFT son iguales a la cantidad de muestras de la señal N=L.

$X[k] = N \delta [k-k_0]$

Para N=2L, la gráfica adquiere mayor resolución en el dominio de la frecuencia, donde se requiere completar los puntos L para que sean iguales a N y realizar el cálculo de fft.

Se crea a partir de una función exponencial multiplicada por una rectangular con ancho igual a un periodo. Para la gráfica se considera solo la parte real de la función.

# INGRESO
L  = 20   # Tramos en un periodo

# señal es exponencial complejo periódico
f0 = 0.2  # frecuencia fundamental de señal
ancho = 1*(1/f0) # ancho del impulso rectangular
rect = lambda t: np.heaviside(t,1)-np.heaviside(t-ancho,1)
xt = lambda t: rect(t)*np.exp(1j*2*np.pi*f0*t)

# para DFT
N = 2*L # N par mayor que L, tramos DFT

DFT: [ impulso unitario ] [ pulso rectangular ] [ exponencial complejo ] [ algoritmo ]
..

---

4. Algoritmo en Python

Se adjunta el algoritmo para los ejercicios, el bloque de ingreso presentado corresponde al ejercicio de exponencial complejo. Para los otros ejercicios actualizar el bloque de INGRESO

# ejemplo 8.1.2.2 p307 DFT x[n] zero-padding
# telg1034 DSP fiec-espol edelros@espol.edu.ec
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
L  = 20   # Tramos en un periodo

# señal es exponencial complejo periódico
f0 = 0.2  # frecuencia fundamental de señal
ancho = 1*(1/f0) # ancho del impulso rectangular
rect = lambda t: np.heaviside(t,1)-np.heaviside(t-ancho,1)
xt = lambda t: rect(t)*np.exp(1j*2*np.pi*f0*t)

# para DFT
N = 2*L # N par mayor que L, tramos DFT

# PROCEDIMIENTO
# muestreo de señal x[n]
t0 = 0 ; tn = 1/f0 # intervalo [t0,tn]
dt = (tn-t0)/L
ti = np.arange(t0,tn,dt)
xn = xt(ti)

# verificar que L<=N
if L<=N:
    # DFT transformada rapida de Fourier
    Xk = np.fft.fft(xn,n=N)
    Xkmagn = np.abs(Xk) # para gráfica
    Xkfase = np.angle(Xk)
    freq_Hz = np.fft.fftfreq(N,dt) # frecuencias en Hz

    mitad = int(N/2)
    k0 = np.argmax(Xk) # Pico |X[k]|

# SALIDA
if N>=L:
    print('k0:',k0,'; i:',mitad+k0)
    print('freq_Hz:',freq_Hz[k0],
          '; Xk[k0]:',Xk[k0])
    print('[ i, freq[ki],|X[k]| ,X[k]]')
    kmin = k0-4 ; kmax = k0+4
    if kmin<0:
        kmin = 0
    if kmax>N:
        kmax = N
    for i in range(kmin,kmax,1):
        print([i,freq_Hz[i],Xkmagn[i],Xk[i]])
else:
    print('Revisar: L debe ser menor o igual a N')
    print('L:',L,'; N:',N)

# GRAFICA
if N>=L:
    LN_texto= ', L='+str(L)+', N='+str(N)
    zp_texto = 'zero padding:'+str(N-L)
    plt.subplot(211)
    plt.stem(ti,np.real(xn),label='Re(x[n])',
             markerfmt ='C2.',linefmt='C2-')
    plt.stem(ti,np.imag(xn),label='Im(x[n])',
             markerfmt ='C4.',linefmt='C4:')
    plt.stem(L*dt,0,markerfmt ='C3.')
    plt.xlabel('ti (s)')
    
    plt.title('x[n] vs |x[k]|'+LN_texto+', '+zp_texto)
    plt.legend()

    plt.subplot(212)
    plt.stem(freq_Hz[0:mitad], Xkmagn[0:mitad],
             label='|X[k]| centrada]',
             markerfmt ='C0.',linefmt='C0-')
    plt.stem(freq_Hz[mitad:], Xkmagn[mitad:],
             markerfmt ='C1.',linefmt='C1-')
    plt.xlabel('freq (Hz)')
    plt.legend()
    plt.tight_layout()
    plt.show()

DFT: [ impulso unitario ] [ pulso rectangular ] [ exponencial complejo ] [ algoritmo ]

 

[ DFT ] [ Algoritmo DFT ] [ IDFT ] [ Algoritmo IDFT ]
..

---

1. La Transformada de Fourier Discreta – DFT

Referencia: McClellan 8.1 p302

La ecuación para la transformada de Fourier discreta o DFT, discreta en tiempo y discreta en frecuencia

$X[k] = \sum_{n=0}^{N-1} x[n] e^{-j(2\pi/N) k n}$

k = 0,1, … N-1

que toma N muestras en el dominio del tiempo y la transforma en N valores de X[k] en el dominio de la frecuencia. Los valores de X[k] son de tipo complejo, mientras que los valores de x[n] son reales la mayoría de las veces.

1.1. DFT – ejemplo para  x[n] de pocas muestras

Referencia: McClellan ejemplo 8.1 p303

Para calcular la DFT de la secuencia x[n] = [1,1,0,0] se aplica la expresión para k=0,1,2,3,  siendo N=4, todos los exponentes son múltiplos de π/2.

$X[0] = x[0] e^{-j(2\pi/4) (0) (0)} + x[1] e^{-j(2\pi/4) (0) (1)}$ $+ x[2] e^{-j(2\pi/4) (0) (2)} + x[3] e^{-j(2\pi/4) (0) (3)}$ $= 1 e^{-j0} + 1 e^{-j0} + 0 e^{-j0} + 0 e^{-j0} = 2$ $X[1] = x[0] e^{-j(2\pi/4) (1) (0)} + x[1] e^{-j(2\pi/4) (1) (1)}$ $+ x[2] e^{-j(2\pi/4) (1) (2)} + x[3] e^{-j(2\pi/4) (1) (3)}$ $= 1 e^{-j0} + 1 e^{-j\pi/2} + 0e^{-j\pi} + 0 e^{-j(3\pi/2) }$ $= 1-j$ $X[2] = x[0] e^{-j(2\pi/4) (2) (0)} + x[1] e^{-j(2\pi/4) (2) (1)}$ $+ x[2] e^{-j(2\pi/4) (2) (2)} + x[3] e^{-j(2\pi/4) (2) (3)}$ $= 1 e^{-j0} + 1 e^{-j\pi} + 0e^{-j2\pi} + 0 e^{-j(3\pi) }$ $= 1+(-1) = 0$ $X[3] = x[0] e^{-j(2\pi/4) (3) (0)} + x[1] e^{-j(2\pi/4) (3) (1)}$ $+ x[2] e^{-j(2\pi/4) (3) (2)} + x[3] e^{-j(2\pi/4) (3) (3)}$ $= 1 e^{-j0} + 1 e^{-j3\pi/2} + 0e^{-j3\pi} + 0 e^{-j(9\pi/2) }$ $= 1+j$

Coeficientes resultantes en el dominio de la frecuencia:

X[k] = [2, 1-j, 0,1+j]

Se presenta el resultado en forma polar:

$X[k] = [2,\sqrt{2} e^{-j\pi/4}, 0 , \sqrt{2} e^{-j\pi/4}]$

[ DFT ] [ Algoritmo DFT ] [ IDFT ] [ Algoritmo IDFT ]

..

---

1.2 Algoritmo DFT con Python

Para el algoritmo, el bloque de ingreso de datos es el vector con las muestras numéricas:

# INGRESO
xn = [1,1,0,0]

En Python se obtiene los coeficientes de la DFT con la instrucción np.fft.fft(xn). El resultado para el ejercicio anterior se muestra como:

X[k]:
[2.+0.j 1.-1.j 0.+0.j 1.+1.j]
Magnitud X[k]:
[2.         1.41421356 0.         1.41421356]
Fase X[k]:
[ 0.         -0.78539816  0.          0.78539816]
>>>

detalle de las instrucciones:

# ejemplo 8.1 p302 DFT de x[n]
# telg1034 DSP fiec-espol edelros@espol.edu.ec
import numpy as np

# INGRESO
xn = [1,1,0,0]

# PROCEDIMIENTO
X_k = np.fft.fft(xn)
magnk = np.abs(X_k)
fasek = np.angle(X_k)

# SALIDA
print('X[k]:')
print(X_k)
print('Magnitud X[k]:')
print(magnk)
print('Fase X[k]:')
print(fasek)

[ DFT ] [ Algoritmo DFT ] [ IDFT ] [ Algoritmo IDFT ]
..

---

2. Transformada Inversa de Fourier Discreta – IDFT

La transformada inversa de Fourier discreta convierte la secuencia X[k] en el dominio de la frecuencia a x[n] en el dominio del tiempo

$x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{j(2\pi/N) k n}$

n = 0, 1, … N-1

2.1 IDFT – Ejemplo para  X[k] de pocas muestras

Referencia: McClellan ejemplo 8.2 p305

Para calcular la IDFT de la secuencia X[k] obtenida en el ejercicio anterior, se aplica la fórmula con n=0,1,2,3,  siendo N=4, todos los exponentes son múltiplos de π/2.

$X[k] = [2,\sqrt{2} e^{-j\pi/4}, 0 , \sqrt{2} e^{-j\pi/4}]$

siendo:

$x[0] = \frac{1}{4} \Big( X[0] e^{j(2\pi/4) (0) (0)} + X[1] e^{j(2\pi/4) (1) (0)}$ $+ X[2] e^{j(2\pi/4) (2) (0)} + X[3] e^{j(2\pi/4) (3) (0)} \Big)$ $= \frac{1}{4} \Big( 2 e^{j0} + \sqrt{2} e^{-j\pi/4} e^{j0}$ $+ 0 e^{j0} + \sqrt{2} e^{-j\pi/4} e^{j0} \Big)$ $= \frac{1}{4} \Big( 2 +(1-j)+0+(1+j) \Big) = 1$ $x[1] = \frac{1}{4} \Big( X[0] e^{j(2\pi/4) (0) (1)} + X[1] e^{j(2\pi/4) (1) (1)}$ $+ X[2] e^{j(2\pi/4) (2) (1)} + X[3] e^{j(2\pi/4) (3) (1)} \Big)$ $= \frac{1}{4} \Big( 2 e^{j0} + \sqrt{2} e^{-j\pi/4} e^{j\pi/2}$ $+ 0 e^{j\pi/2} + \sqrt{2} e^{-j\pi} e^{j3\pi/2} \Big)$ $= \frac{1}{4} \Big( 2 +(1+j)+0+(1-j) \Big) = 1$ $x[2] = \frac{1}{4} \Big( X[0] e^{j(2\pi/4) (0) (2)} + X[1] e^{j(2\pi/4) (1) (2)}$ $+ X[2] e^{j(2\pi/4) (2) (2)} + X[3] e^{j(2\pi/4) (3) (2)} \Big)$ $= \frac{1}{4} \Big( 2 e^{j0} + \sqrt{2} e^{-j\pi/4} e^{j\pi}$ $+ 0 e^{j2\pi} + \sqrt{2} e^{-j\pi} e^{j3\pi} \Big)$ $= \frac{1}{4} \Big( 2 +(-1+j)+0+(-1-j) \Big) = 0$ $x[3] = \frac{1}{4} \Big( X[0] e^{j(2\pi/4) (0) (3)} + X[1] e^{j(2\pi/4) (1) (3)}$ $+ X[2] e^{j(2\pi/4) (2) (3)} + X[3] e^{j(2\pi/4) (3) (3)} \Big)$ $= \frac{1}{4} \Big( 2 e^{j0} + \sqrt{2} e^{-j\pi/4} e^{j3\pi/2}$ $+ 0 e^{j3\pi} + \sqrt{2} e^{-j\pi} e^{j9\pi/2} \Big)$ $= \frac{1}{4} \Big( 2 +(-1-j)+0+(-1+j) \Big) = 0$

Coeficientes resultantes en el dominio del tiempo:

x[n] = [1, 1, 0,0]

Con lo que se verifica que es la inversa del ejercicio anterior.

[ DFT ] [ Algoritmo DFT ] [ IDFT ] [ Algoritmo IDFT ]

..

---

2.2 Algoritmo IDFT con Python

Para el ejercicio, la secuencia de entrada es X[k] = [2, 1 + j, 0, 1 – j]. Usando la librería Numpy, la parte compleja de un número se expresa como 1j  el vector de ingreso se expresa como:

# INGRESO
Xk = [2, 1 - 1j, 0, 1 + 1j]

el resultado del algoritmo obtenido es:

x[n]:
[1. 1. 0. 0.]

Las instrucciones se desarrollan de forma semejante a las usadas para DFT. La inversa de la transformada, se obtiene con:

# ejemplo 8.2 p306 IDFT de X[k]
# telg1034 DSP fiec-espol edelros@espol.edu.ec
import numpy as np

# INGRESO
Xk = [2, 1 - 1j, 0, 1 + 1j]

# PROCEDIMIENTO
xn = np.fft.ifft(Xk)

# simplifica xn si es solo real
xn_re = np.real(xn)
xn_im = np.imag(xn)

if np.sum(xn_im)==0:
    xn= xn_re

# SALIDA
print('x[n]:')
print(xn)

[ DFT ] [ Algoritmo DFT ] [ IDFT ] [ Algoritmo IDFT ]

DFT/IDFT – Transformada de Fourier Discreta con Python

DFT – Pares DFT ↔ DTFT con Python

DFT – Propiedades – FFT con Python

DFT – Tablas de Transformadas de Fourier Discretas