DFT

DFT using complex sinusoids



In [1]:

    
import numpy as np
import matplotlib.pyplot as plt

N-point DFT of an input $x[n]$ is given by:

$X[k] = \sum\limits_{0}^{N-1} x(n) \ \mathit{e}^{-j2 \pi n k/N} , \ \ k = 0, \cdots, N-1$

where:

n = discrete time instances
k = discrete frequency index
$\omega_k = 2 \pi k / N $

Also, $X[k] = \sum\limits_{0}^{N-1} x(n) \ \mathit{e}^{-j2 \pi n k/N} = \langle x, s^{*}_{k} \rangle$

i.e : DFT is the projection of I/P signal onto a finite set of conjugate complex sinusoids

where:

$s^{*}_{k} = \mathit{e}^{-j2 \pi n k/N} = \cos(2 \pi n k/N) - j \sin(2 \pi n k/N)$, is the complex exponential.
$s^{*}_{k}$ forms the $k^{th}$ complex basis function, and is of the same size as $x(n)$
A non-zero $X[k]$ implies presence of $k^{th}$-frequency in $x(n)$

NOTE:

$N$ = (DFT-size) is the period of the DFT signal
$k$ specifies the number of periods

DFT of real-valued sinusoids

$ X[k] = N \frac{A_{0}}{2}, \ \ \text{for} \ k = -k_{0}, k_{0}; \text{zero elsewhere}$



In [42]:

    
N = 512
n = np.arange(-N/2, N/2)
Fs = 16000
f = 100
A = 0.75

# I/P signal
x1 = A*np.exp(1j*2*np.pi*f*n/Fs)
x2 = A*np.sin(2*np.pi*f*n/Fs)

plt.subplot(221)
plt.plot(n,np.real(x))
plt.grid('on')
plt.title("I/P signal 1 (complex)")

plt.subplot(222)
plt.plot(n,np.real(x))
plt.grid('on')
plt.title("I/P signal 2 (real)")

# DFT
def dft(x):
    X = []
    for k in np.arange(N):
        s = np.exp(-1j*2*np.pi*k*n/N)
        X = np.append(X, np.sum(x*s))
    return X

# get DFT of x
X1 = dft(x1)
X2 = dft(x2)

plt.subplot(223)
plt.plot(n,np.abs(X1)/N)
plt.grid('on')
plt.title("DFT1")

plt.subplot(224)
plt.plot(n,np.abs(X2)/N)
plt.grid('on')
plt.title("DFT2")

plt.tight_layout()
plt.show()

Inverse DFT

$ x(n) = \langle X, s_{k} \rangle = \frac{1}{N} \sum \limits_{k = 0}^{N-1} X[k] s_{k}[n] \ \ \text{for} \ k = 0, \cdots, N-1$

IDFT is dot product of spectral components and complex sinusoids
Comparing DFT and IDFT, IDFT does not include conjugates of sinunoids, and has normalizing factor

Sinusoids (time-domain)

Real sinusoid:

$ x[n] = A \sin (\omega n T_{s} + \phi) = A \sin(2 \pi f n T_{s} + \phi)$



In [27]:

    
A = 0.8
f = 1000
Fs = 44100
phi = np.pi/2
n = np.arange(-0.002,0.002,1.0/Fs)
print n[:10], n.size

# signal
x = A*np.cos(2*np.pi*f*n + phi)

plt.plot(n,x)
plt.xlabel('time')
plt.ylabel('amplitude')
plt.grid('on')
plt.show()









    



[-0.002      -0.00197732 -0.00195465 -0.00193197 -0.0019093  -0.00188662
 -0.00186395 -0.00184127 -0.00181859 -0.00179592] 177

Complex sinusoid

$ x[n] = A e^{(j \omega n T_{s} + \phi)} = A e^{(j 2 \pi f n T_{s} + \phi)} $

where:

$x[n]$ = array of complex values of the sinusoid
$n$ = discrete time index
$A$ = amplitude
$T_{s}$ = sampling interval
$\phi$ = initial phase
$f$ = frequency of the sinusoid

Complex sinusoid basis:

$ s_{k}[n] = e^{2 \pi k n/N} = \cos(2 \pi k n/N) + j \sin (2 \pi k n/N)$



In [24]:

    
N = 512
A = 0.8
k = 3
n = np.arange(-N/2, N/2)

# signal
x = A*np.exp(1j* 2*np.pi*k*n/N)

plt.subplot(121)
plt.plot(n,np.real(x))
plt.xlabel('time')
plt.ylabel('amplitude')
plt.title('Real part')
plt.grid('on')

plt.subplot(122)
plt.plot(n,np.imag(x))
plt.xlabel('time')
plt.ylabel('amplitude')
plt.title('Imaginary part')
plt.grid('on')

plt.tight_layout()
plt.axis([-N/2, N/2, -1, 1])
plt.show()



In [10]:

    
N = 8
n = np.arange(-N/2, N/2)
Fs = 16000
f = 0.5
A = 1

# I/P signal
x1 = A*np.sin(2*np.pi*f*n/Fs)

plt.subplot(211)
plt.plot(n,np.abs(x1))
plt.grid('on')
plt.title("I/P signal 1 (complex)")

# DFT
def dft(x):
    X = []
    for k in np.arange(N):
        s = np.exp(-1j*2*np.pi*k*n/N)
        X = np.append(X, np.sum(x*s))
    return X

# get DFT of x
X1 = dft(x1)

plt.subplot(212)
plt.plot(n,np.abs(X1)/N)
plt.grid('on')
plt.title("DFT1")

plt.tight_layout()
plt.show()

IDFT



In [18]:

    
k = np.arange(0,4)
X = np.array([4,0,0,0])
N = np.size(X)

def idft(X):
    x = []
    for n in np.arange(np.size(X)):
        s = np.exp(1j*2*np.pi*k*n/N)
        x = np.append(x, np.sum(X*s))
    return x/N

print idft(X)









    



[ 1.+0.j  1.+0.j  1.+0.j  1.+0.j]