Fourier Compression Demo¶

In [1]:
# Imports
import numpy as np
from numpy.fft import fft, ifft, fft2, ifft2, fftshift, ifftshift
import matplotlib.pyplot as plt
from scipy import ndimage, misc
from scipy.signal import gaussian
gray = plt.cm.gray

Compress a 1-D Signal¶

In [4]:
# Read a signal in
img = plt.imread('images/pd.jpg')
f = np.array(img[128,:,1], dtype=float)
N = len(f)
t = range(N)
In [5]:
plt.plot(f)
plt.title('Original Signal');
In [7]:
# Compute the DFT
F = fft(f)
In [8]:
# Let's look at the first few Fourier coefficients.
F[3]
Out[8]:
(1773.0144397941126-145.8992067869071j)
In [9]:
F[-3]
Out[9]:
(1773.0144397941128+145.89920678690703j)
In [10]:
# Input was real-valued, so our Fourier coefficients are conjugate-symmetric.
idx = 5
print(F[idx])
print(F[256-idx])

(-368.9844443497931-351.85254750210095j)
(-368.9844443497932+351.8525475021009j)
In [11]:
plt.stem(abs(F))
plt.title('DFT of Signal')
plt.xlabel('Freq. Index')
plt.ylabel('Modulus');
In [12]:
ctr = int(np.floor(N/2.))
omega = range(-128,128)
plt.stem(omega, fftshift(abs(F)))
plt.title('DFT of Signal')
plt.xlabel('Freq. Index')
plt.ylabel('Modulus');

Filter out many of the coefficients¶

In [13]:
T = 10
G = F.copy()
G[abs(ifftshift(omega))>T] = 0.
In [14]:
plt.stem(omega, fftshift(abs(G)))
plt.title('DFT of Signal')
plt.xlabel('Freq. Index')
plt.ylabel('Modulus');
In [15]:
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¶

In [16]:
T = 40
G = F.copy()
G[abs(ifftshift(omega))>T] = 0.
In [17]:
plt.stem(omega, fftshift(abs(G)))
plt.title('DFT of Signal')
plt.xlabel('Freq. Index')
plt.ylabel('Modulus');
In [37]:
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¶

In [18]:
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);
In [19]:
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¶

In [20]:
r = 40.  # 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');

Remaining coefs: 7.66754150390625%
In [21]:
# What do only the high frequencies of an image look like?
r = 40.  # Keep coefficients of frequency more than this
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');

Remaining coefs: 92.35076904296875%
In [ ]: