Aliasing Demo¶

In [30]:
import numpy as np
from copy import deepcopy
import matplotlib.pyplot as plt
from numpy.fft import fft, ifft, fft2, ifft2, fftshift, ifftshift
In [31]:
def Gaussian2D(dims, std):
    '''
        Creates an array containing a 2-D Gaussian function, with mean
        at mu, and standard deviation std.

        Inputs:
          dims is a 2-vector containing the number of rows and columns
          std is the standard deviation
    '''
    nr, nc = dims
    rr, cc = np.mgrid[0:nr,0:nc]
    ctr = GetCentre(rr)
    temp1 = rr - ctr[0]
    temp2 = cc - ctr[1]
    temp1 = np.exp(-temp1**2/std)
    temp2 = np.exp(-temp2**2/std)
    blah = temp1*temp2
    return blah/sum(blah.flatten())

def GetCentre(f):
    return np.array(np.floor(np.array(np.shape(f))/2.), dtype=int)

Sample Aliasing¶

Why do indecies above N/2 correspond to lower (negative) frequencies?

In [32]:
N = 128
tt = np.linspace(0, N, 3000)
In [33]:
# Start with a low frequency index
k = 5
F4 = np.zeros(N);
F4[k] = N
f4 = ifft(F4)
plt.figure(figsize=(14,3))
plt.subplot(1,2,1); plt.plot(tt, np.cos(2.*k*np.pi*tt/N), 'r-'); plt.plot(np.real(f4),'.', markersize=8);
plt.subplot(1,2,2); plt.plot(np.real(f4),'.', markersize=8);
In [34]:
k = 7
F4 = np.zeros(N);
F4[k] = N
f4 = ifft(F4)
plt.figure(figsize=(14,3))
plt.subplot(1,2,1); plt.plot(tt, np.cos(2.*k*np.pi*tt/N), 'r-'); plt.plot(np.real(f4),'.', markersize=8);
plt.subplot(1,2,2); plt.plot(np.real(f4),'.', markersize=8);
In [35]:
k = 15
F4 = np.zeros(N);
F4[k] = N
f4 = ifft(F4)
plt.figure(figsize=(14,3))
plt.subplot(1,2,1); plt.plot(tt, np.cos(2.*k*np.pi*tt/N), 'r-'); plt.plot(np.real(f4),'.', markersize=8);
plt.subplot(1,2,2); plt.plot(np.real(f4),'.', markersize=8);
In [36]:
k = 125
F4 = np.zeros(N);
F4[k] = N
f4 = ifft(F4)
plt.figure(figsize=(14,3))
plt.subplot(1,2,1); plt.plot(tt, np.cos(2.*k*np.pi*tt/N), 'r-'); plt.plot(np.real(f4),'.', markersize=8);
plt.subplot(1,2,2); plt.plot(np.real(f4),'.', markersize=8);
In [37]:
k = 3
F4 = np.zeros(N);
F4[k] = N
f4 = ifft(F4)
plt.figure(figsize=(14,3))
plt.subplot(1,2,1); plt.plot(tt, np.cos(2.*k*np.pi*tt/N), 'r-'); plt.plot(np.real(f4),'.', markersize=8);
plt.subplot(1,2,2); plt.plot(np.real(f4),'.', markersize=8);

Aliasing¶

Once your index gets past $\frac{N}{2}$, the Fourier coefficients start to correspond to lower and lower frequencies. The plot below shows what frequency corresponds to each index. This phenomenon of higher frequencies looking like lower frequencies is called "aliasing".

In [18]:
F = np.zeros(8192)
N = len(F)      # Number of samples in time (or space)
Omega = 8192    # Sound is sampled at 'Omega' samples per second (Hz)
L = N / Omega   # The sound clip will be L seconds long

# To get the corresponding freqencies for the Fourier coefficients,
# we can use these lines.
omega_index = np.arange(0, N, dtype=float)
omega = deepcopy(omega_index)
omega[omega>=np.floor(N/2)] -= N

plt.figure(figsize=(10,6));
plt.plot(omega_index, color='lightblue')
plt.axvline(color='lightgray'); plt.axhline(color='lightgray')
plt.plot([0, N-1], [N/2, N/2], '--', color='lightgray')
plt.plot([0, N-1], [-N/2, -N/2], '--', color='lightgray')
plt.plot(omega, 'bo'); plt.xlabel('Index'); plt.ylabel('Frequency (Hz)');
In [19]:
F = np.zeros(8192)
N = len(F)      # Number of samples in time (or space)
Omega = 8192    # Sound is sampled at 'Omega' samples per second (Hz)
L = N / Omega   # The sound clip will be L seconds long

# To get the corresponding freqencies for the Fourier coefficients,
# we can use these lines.
omega_index = np.arange(0, N, dtype=float)
omega = deepcopy(omega_index)
omega[omega>=np.floor(N/2)] -= N

plt.figure(figsize=(10,6));
#plt.plot(fftshift(omega_index), color='lightblue')
plt.axvline(color='lightgray'); plt.axhline(color='lightgray')
plt.plot([np.min(omega), np.max(omega)], [N/2, N/2], '--', color='lightgray')
plt.plot([np.min(omega), np.max(omega)], [-N/2, -N/2], '--', color='lightgray')
plt.plot(omega, omega, 'bo'); plt.xlabel('Shifted Index'); plt.ylabel('Frequency (Hz)');

Let's hear the aliasing¶

In [20]:
from IPython.display import Audio
idx = 440   # frequency idx/L Hz
F = np.zeros(N)
F[idx] = 1.1
f = np.real(ifft(F))  # IDFT
f = f / max(f)  # normalize so we don't blow out our speakers

# Plot it
plt.plot(omega, abs(F) )
plt.title(str(omega[idx])+'Hz Tone')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Modulus')

#PlotFT(f, Omega, 1)
#scipy.io.wavfile.write('tone.wav', Omega, f)
Audio(f, rate=Omega)
Out[20]:
Your browser does not support the audio element.
In [21]:
Omega
Out[21]:
8192
In [22]:
from IPython.display import Audio
idx = 4000   # frequency idx/L Hz
F = np.zeros(N)
F[idx] = 1.1
f = np.real(ifft(F))  # IDFT
f = f / max(f)  # normalize so we don't blow out our speakers

# Plot it
plt.plot(omega, (abs(F)) )
plt.title(str(omega[idx])+'Hz Tone')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Modulus')

#PlotFT(f, Omega, 1)
#scipy.io.wavfile.write('tone.wav', Omega, f)
Audio(f, rate=Omega)
Out[22]:
Your browser does not support the audio element.
In [23]:
from IPython.display import Audio
idx = 6000   # frequency idx/L Hz
F = np.zeros(N)
F[idx] = 1.1
f = np.real(ifft(F))  # IDFT
f = f / max(f)  # normalize so we don't blow out our speakers

# Plot it
plt.plot(omega, abs(F) )
plt.title(str(omega[idx])+'Hz Tone')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Modulus')

#PlotFT(f, Omega, 1)
#scipy.io.wavfile.write('tone.wav', Omega, f)
Audio(f, rate=Omega)
Out[23]:
Your browser does not support the audio element.

Alising in images¶

In [24]:
x, y = np.mgrid[0:256,0:256]
f = np.sin(2.*np.pi*x/256*40)
#f = ndimage.imread('pd.jpg')
plt.figure(figsize=(15,7))
plt.clf()
plt.subplot(1,3,1)
plt.imshow(f, cmap='gray');
c = plt.axis()
plt.title('Original');

g = f[::4,::4]  # every fifth sample
plt.subplot(1,3,2)
plt.imshow(g, cmap='gray'); #plt.axis(c);
plt.title('Subsampled');

plt.subplot(1,3,3);
plt.imshow(g, cmap='gray'); plt.axis(c);
plt.title('Subsampled');
In [26]:
q = plt.imread('images/bricks.jpg')
q = np.array(q[:,:,0], dtype=float)
nr, nc = np.shape(q)
In [27]:
plt.figure(figsize=(15,7));
plt.subplot(1,3,1)
plt.imshow(q, cmap='gray', interpolation='none')
plt.title('Original'); c = plt.axis();

q_small = q[::5,::5]
qsr, qsc = np.shape(q_small)
plt.subplot(1,3,2)
plt.imshow(q_small, cmap='gray', interpolation='nearest')
plt.axis(c);
plt.title('Subsampled, with aliasing');

plt.subplot(1,3,3)
plt.imshow(q_small, cmap='gray', interpolation='nearest')
plt.title('Subsampled, with aliasing');
In [28]:
plt.figure(figsize=(15,7));
plt.subplot(1,3,1)
plt.imshow(q, cmap='gray', interpolation='none')
plt.title('Original');

q_small = q[::5,::5]
qsr, qsc = np.shape(q_small)

g = Gaussian2D(q.shape, 10)   # Create filter
Q = fftshift( fft2(q) )       # DFT of the filter

G = fftshift( fft2(fftshift(g)) )
h = ifft2(ifftshift(Q*G))     # Apply filter in the freq domain

plt.subplot(1,3,2)
plt.imshow(np.real(h[::5,::5]), cmap='gray', interpolation='nearest')
plt.axis(c);
plt.title('Smoothed and subsampled (no aliasing)');

plt.subplot(1,3,3)
plt.imshow(np.real(h[::5,::5]), cmap='gray', interpolation='nearest')
plt.title('Smoothed and subsampled (no aliasing)');
In [ ]: