Fun with FFT

This notebook is part of a lecture I gave in my DSP class on Halloween, Oct 31. I'll clean it up later, but need to post it (even in its current minimalist condition) for my students.

%pylab inline
from scipy.fftpack import fft, ifft

Populating the interactive namespace from numpy and matplotlib

x = array([1,-1,1,-1])
fft(x)

array([ 0.+0.j,  0.+0.j,  4.+0.j,  0.-0.j])

fft(x,n=8)

array([ 0.+0.j        ,  1.+0.41421356j,  0.-0.j        ,  1.+2.41421356j,
        4.+0.j        ,  1.-2.41421356j,  0.+0.j        ,  1.-0.41421356j])

xz = r_[x,zeros(4)]
xz

array([ 1., -1.,  1., -1.,  0.,  0.,  0.,  0.])

fft(xz)

array([ 0.+0.j        ,  1.+0.41421356j,  0.-0.j        ,  1.+2.41421356j,
        4.+0.j        ,  1.-2.41421356j,  0.+0.j        ,  1.-0.41421356j])

c_[x,zeros(4)]

array([[ 1.,  0.],
       [-1.,  0.],
       [ 1.,  0.],
       [-1.,  0.]])

y = fft(x)
xhat = ifft(y)

x-xhat

array([ 0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j])

x2 = rand(1024*1024)
y2 = fft(x2)
x2hat = ifft(y2)
r = x2-x2hat
r[:10]

array([  1.11022302e-16 +0.00000000e+00j,
         2.22044605e-16 +5.14996032e-19j,
         2.22044605e-16 +4.07962330e-17j,
         2.22044605e-16 +5.33427469e-17j,
         1.11022302e-16 +4.61328024e-17j,
        -1.94289029e-16 -4.28259858e-17j,
         1.66533454e-16 +2.98801962e-17j,
         2.22044605e-16 -5.40169851e-17j,
         0.00000000e+00 -3.37720897e-17j,   3.33066907e-16 +8.51658483e-18j])

max(abs(x2-x2hat))

1.0431292286976892e-15

x2hat[:10]

array([ 0.41116206 +0.00000000e+00j,  0.33706664 -5.14996032e-19j,
        0.79717000 -4.07962330e-17j,  0.53664741 -5.33427469e-17j,
        0.79604319 -4.61328024e-17j,  0.19359754 +4.28259858e-17j,
        0.11685530 -2.98801962e-17j,  0.50703792 +5.40169851e-17j,
        0.35640277 +3.37720897e-17j,  0.73290768 -8.51658483e-18j])

real(x2hat)[:10]

array([ 0.41116206,  0.33706664,  0.79717   ,  0.53664741,  0.79604319,
        0.19359754,  0.1168553 ,  0.50703792,  0.35640277,  0.73290768])

#x = array([0,0,0,1,0,0,0,0]) + 1j*array([0,0,0,0,0,0,0,0])
x = exp(1j*2.1*2*pi*arange(8)/8)
y = fft(x)
subplot(1,2,1)
stem(real(y))
subplot(1,2,2)
stem(imag(y))
xlim([0,8])

(0, 8)

N=64
omega = 2.3*2*pi/N
n = arange(N)
x = cos(omega*n)
y = fft(x)
stem(abs(y))

<Container object of 3 artists>

Signal Filtering with FFTs

N=256
sigma = 1.0
omega = 2*pi*35.5/N
omega2 = 2*pi*44/N
#x = exp(omega*1j*arange(N))
x = cos(omega*arange(N))+2*cos(omega2*arange(N))
x += sigma*randn(N)
y = fft(x)
plot(abs(y))

[<matplotlib.lines.Line2D at 0x113865d10>]

plot(x)

[<matplotlib.lines.Line2D at 0x1138931d0>]

from scipy.signal import kaiser
window = kaiser(N,beta=14)
plot(abs(fft(x*window)))

[<matplotlib.lines.Line2D at 0x11309d3d0>]

plot(abs(fft(window)))

[<matplotlib.lines.Line2D at 0x113514690>]

plot(window)

[<matplotlib.lines.Line2D at 0x115bf6ad0>]

Dirichlet Kernel, the Discrete Time Sinc Function

The Dirichlet kernel is the discrete time version of the sinc function. For some reason, it gets little of the acclaim of the sinc function. I've even seen textbooks incorrectly use the sinc function when they should be using the Dirichlet function. Even the Wikipedia entry for the sinc function make no mention of the Dirichlet kernel.

The first version uses the piecewise function to take care of the division by 0 problem. However, as written it doesn't work for x being any other integer except 0. I suppose it could be improved, but Dirichlet2 solves the problem a different (and arguably better) way: it allows numpy to compute nan's, then replaced the nan's by the correct limit.

def dirichlet(x,N):
    return piecewise(x, [x==0, x != 0], [N, lambda x: sin(N*pi*x)/sin(pi*x)])

x = linspace(-0.5,0.5,201)
plot(x,abs(dirichlet(x,16)))

[<matplotlib.lines.Line2D at 0x116bdbc10>]

def dirichlet2(x,N):
    seterr(divide='ignore')
    y = sin(N*pi*x)/sin(pi*x)
    y[isnan(y)] = N
    return y

plot(abs(dirichlet2(linspace(-0.5,0.5,201),16)))

[<matplotlib.lines.Line2D at 0x11743b7d0>]