Skip to content

Instantly share code, notes, and snippets.

Embed
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)
@dmishin

This comment has been minimized.

Copy link
Owner Author

commented Oct 23, 2014

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.