#! /usr/bin/env python3
# -*- coding:utf8 -*-

import matplotlib.pyplot as plt
import numpy as np
from math import *

plt.rcParams['font.sans-serif'] = 'DejaVu Sans'
np.random.seed(1)

"""
Note: Blue noise is actually not so well depictable. As more higher frequency
componentes are added, the amplitude goes to infinity and will be dominated
by high-frequency noise. Thus, the image changes a lot with the cutoff sampling
density. We choose roughtly half a linewidth for sampling for a decent
appearance.
"""

nsamples = 501
t0 = 1.
t = np.linspace(0, t0, nsamples)

dt = t[1:] - t[:-1]
white_noise = np.random.normal(0, 1, nsamples)
fourier_amplitudes = np.fft.rfft(white_noise)
frequencies = np.fft.rfftfreq(nsamples, d=t[1] - t[0])

fourier_amplitudes *= np.sqrt(frequencies) # sqrt(f) amplitude spectrum
X = np.fft.irfft(fourier_amplitudes, n=nsamples, norm='ortho')
X /= np.std(X)

fig = plt.figure(figsize=(520 / 90.0, 340 / 90.0), dpi=72)
plt.plot(t, X)
plt.grid(True)
plt.xlim(t[0], t[-1])
plt.ylim(-3.1, 3.1)
plt.xlabel('t')
plt.ylabel('X')
plt.tight_layout()
plt.savefig('Blue-noise-trace.svg')