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

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

def ramlak(N):
    y = np.zeros(N)
    y[0] = 0.25
    fodd = np.arange(1, N, 2)
    fodd = np.minimum(fodd, N - fodd)
    y[1:N:2] = (-1. / (pi * pi)) / (fodd * fodd)
    return y

N = 16
xlim = -6, 6
kernel = ramlak(N)

fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(7, 3.5))

ax0.plot(np.arange(-N//2, N - N//2),
         np.roll(kernel, N//2),
         "o-", label="Ram-Lak")

xa = np.linspace(*xlim, 801)
ya = np.divide(-pi**-2, xa * xa, out=None, where=(xa!=0))
ya[xa==0] = float(nan)
ax0.plot(xa, ya, "--", color="#bbbbbb", zorder=1.9, label=r"$-1/(x \pi)^2$")

annotate_args = {"va":"center", "ha":"left", "xytext":(12, 0), "zorder":1.8,
    "textcoords":"offset points", "fontsize":14,
    "bbox":{"fc":"white", "ec":"none"} }
ax0.annotate(r"$\frac{1}{4}$", (0, kernel[0]), **annotate_args)
ax0.annotate(r"$\frac{-1}{\pi^2}$", (1, kernel[1]), **annotate_args)
annotate_args.update({"va":"top", "ha":"center", "xytext":(0, -8)})
ax0.annotate(r"$\frac{-1}{(3\pi)^2}$", (3, kernel[3]), **annotate_args)
ax0.annotate(r"$\frac{-1}{(5\pi)^2}$", (5, kernel[5]), **annotate_args)


ax0.grid(True)
ax0.set_xlim(*xlim)
ax0.set_ylim(-0.13, 0.28)
ax0.set_xlabel("x")
ax0.legend(loc="upper left", framealpha=1, borderpad=0.6)

kernel_spectrum = np.fft.fft(kernel, n=N).real
xa2 = np.arange(1-N//2, N//2+1)
ax1.plot(xa2, N * np.roll(kernel_spectrum, N//2-1), "o-", label="FT(Ram-Lak)")
ax1.plot(xa2, N * np.roll(np.minimum(np.arange(N), N - np.arange(N)) / N, N//2-1), "--", color="#888888", zorder=1.9, label="|k|")
ax1.grid(True)
ax1.set_xlim(-4, 4)
ax1.set_ylim(-0.03 * N, 0.28 * N)
ax1.set_xlabel("k")
ax1.set_ylabel(r"$a_k\cdot N$")
ax1.legend(loc="upper center", framealpha=1, borderpad=0.6)

plt.tight_layout()
plt.savefig("Ram-Lak kernel + fourier.svg")