Generate Fourier transform of the coprime integers map

coprime_fft.py

#Requires: numpy, PIL
from fractions import gcd
from numpy import array, meshgrid, frompyfunc, clip, histogram
from numpy.fft import fft2
import numpy as np
from PIL import Image

def value_diapason(x, percent=0.95, nbins=100):
    """Use histogram to determine area, covering 95% of values"""
    counts, bins = histogram(x.ravel(),nbins)
    total = sum(counts)
    accum = 0
    low = bins[-1]
    high = bins[0]
    for i, cnt in sorted( enumerate(counts), 
                          key = (lambda i_c: i_c[1]),
                          reverse=True ):
        accum += cnt
        low = min(low, bins[i])
        high = max(high, bins[i+1])
        if accum > percent * total:
            break
    return low, high

############## Parameters #################
N = 2048
gamma = 1.9
output = r"coprime-fft.png"

############## Calculation ###############
t = array(range(N), dtype=int)
x,y = meshgrid(t,t)
g = frompyfunc(gcd, 2, 1)(x,y) == 1
lfg = np.log(np.abs(fft2(g))) - np.log(N**2)
low, high  = value_diapason( lfg )
img = Image.fromarray( (clip((lfg-low)/(high-low),0,1)**gamma*255).astype(np.uint8), "P")
img.save(output)

See the explanation here: http://dmishin.blogspot.ru/2014/10/fourier-transform-of-coprime-numbers-map.html