Fourier Compression Demo

# Imports
import numpy as np
from numpy.fft import fft, ifft, fft2, ifft2, fftshift, ifftshift
import matplotlib.pyplot as plt
%matplotlib inline
from scipy import ndimage, misc
from scipy.signal import gaussian
gray = plt.cm.gray

Compress a 1-D Signal

# Read a signal in
img = plt.imread('pd.jpg')
f = np.array(img[128,:,1], dtype=float)
N = len(f)
t = range(N)

plt.plot(f)
plt.title('Original Signal');

F = fft(f)

plt.plot(abs(F))
plt.title('DFT of Signal')
plt.xlabel('Freq. Index')
plt.ylabel('Modulus');

ctr = int(np.floor(N/2.))
omega = range(-128,128)
plt.plot(omega, fftshift(abs(F)))
plt.title('DFT of Signal')
plt.xlabel('Freq. Index')
plt.ylabel('Modulus');

Filter out many of the coefficients

T = 10
G = F.copy()
G[abs(ifftshift(omega))>T] = 0.

plt.plot(omega, fftshift(abs(G)))
plt.title('DFT of Signal')
plt.xlabel('Freq. Index')
plt.ylabel('Modulus');

g = ifft(G)
plt.plot(t, f, 'b')
plt.plot(t, np.real(g), 'r')
plt.title('Reconstruction of compressed signal')
plt.xlabel('Spatial Index');

Filter out fewer of the coefficients

T = 40
G = F.copy()
G[abs(ifftshift(omega))>T] = 0.

plt.plot(omega, fftshift(abs(G)))
plt.title('DFT of Signal')
plt.xlabel('Freq. Index')
plt.ylabel('Modulus');

g = ifft(G)
plt.plot(t, f, 'b')
plt.plot(t, np.real(g), 'r')
plt.title('Reconstruction of compressed signal')
plt.xlabel('Spatial Index');

Compress a 2-D Image

f = np.array(img[:,:,0])
f = f + np.random.normal(scale=5.0, size=np.shape(f))
plt.figure(figsize=[7,7])
plt.imshow(f,cmap=gray);

F = fftshift(fft2(f))
plt.figure(figsize=(15,7))
plt.subplot(1,2,1); plt.imshow(abs(F), cmap='gray')
plt.title('Modulus');
plt.subplot(1,2,2); plt.imshow(np.log(abs(F)+1), cmap='gray')
plt.title('log Modulus');

Only keep a circle of Fourier coefficients

r = 10.  # Radius to KEEP
rr, cc = np.mgrid[-128:128,-128:128]
d = np.sqrt(rr**2 + cc**2)
G = F.copy()
G[d>r] = 0.
print('Remaining coefs: '+str((256.**2-sum(G.flatten()==0))/256**2*100)+'%')
g = np.fft.ifft2( np.fft.ifftshift(G) )
plt.figure(figsize=[15,7])
plt.subplot(1,2,1); plt.imshow(np.log(abs(G)+1), cmap='gray')
plt.subplot(1,2,2); plt.imshow(np.real(g), cmap='gray');

r = 80.  # Radius to KEEP
G = F.copy()
G[d>r] = 0.
print('Remaining coefs: '+str((256.**2-sum(G.flatten()==0))/256**2*100)+'%')
g = np.fft.ifft2( np.fft.ifftshift(G) )
plt.figure(figsize=[15,7])
plt.subplot(1,2,1); plt.imshow(np.log(abs(G)+1), cmap='gray')
plt.subplot(1,2,2); plt.imshow(np.real(g), cmap='gray');