Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
Generate Fourier transform of the coprime integers map
#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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment