In [26]:

    

import numpy as np

def make_gabors(N, width, rng=None, octaves=5, F_N=3, F_max=1.0):
    if rng is None:
        rng = np.random.RandomState()
        
    half_N = N / 2
    if half_N * 2 != N:
        half_N += 1
    
    Fs = []
    for i in range(octaves * F_N):
        Fs.append((2**(-i/float(F_N)))*F_max)
        
    f = rng.choice(Fs, size=half_N)    
        
    theta = rng.randint(0, 20, half_N) * (2 * np.pi / 20)
    psi = rng.uniform(0, 2 * np.pi, half_N)
    phi = 1.5
    gamma = rng.uniform(0.4, 0.8, half_N)
    x_offset = rng.uniform(-0.5, 0.5, half_N)
    y_offset = rng.uniform(-1, 1, half_N)
    
    kappa = sqrt(2*np.log(2))*(2.0**phi+1)/(2.0**phi-1)
    sigma = kappa/(2*np.pi*f)
    f = 2*f
        
    x = np.linspace(-1, 1, width)
    X, Y = np.meshgrid(x, x)
    X.shape = width * width
    Y.shape = width * width
    
    X = X[None,:] - x_offset[:,None]
    Y = Y[None,:] + y_offset[:,None]
    
    cosTheta = np.cos(theta)
    sinTheta = np.sin(theta)
    xTheta = X * cosTheta[:,None]  + Y * sinTheta[:,None]
    yTheta = -X * sinTheta[:,None] + Y * cosTheta[:,None]

    envelope = (np.sqrt(2*np.pi)/sigma[:,None])*np.exp(-(np.pi/(2*sigma[:,None]))*(4*(xTheta**2) + yTheta**2))
    re = np.cos(2*np.pi*f[:,None]*xTheta+psi[:,None])-np.exp(-kappa**2/2.0)
    im = np.sin(2*np.pi*f[:,None]*xTheta+psi[:,None])
    
    gabor_re = envelope * re
    gabor_im = envelope * im
        
    gabor = np.vstack([gabor_re, gabor_im])[:N]    
    gabor = gabor / np.linalg.norm(gabor, axis=1)[:, None]
    return gabor

def plot_gabors(data, columns=None):
    if columns is None:
        columns = int(np.sqrt(len(data)))
    pylab.figure(figsize=(10,10))
    vmin = np.min(data)
    vmax = np.max(data)
    width = int(np.sqrt(data.shape[1]))
    for i, d in enumerate(data):
        w = columns - 1 - (i % columns)
        h = i / columns            
        d.shape = width, width
        pylab.imshow(d, extent=(w+0.025, w+0.975, h+0.025, h+0.975), 
                 interpolation='none', vmin=vmin, vmax=vmax, cmap='gray')
        pylab.xticks([])
        pylab.yticks([])
    pylab.xlim((0, columns))
    pylab.ylim((0, len(data) / columns))


    

In [31]:

    

a = make_gabors(1000, width=50, F_max=50)

plot_gabors(a[:25])
show()
import scipy.sparse.linalg
U, S, V = scipy.sparse.linalg.svds(a, k=300)
plot(S[::-1])
show()


    

In [ ]: