import std.stdio;

extern(C) __gshared string[] rt_options =
    [ "gcopt=minPoolSize:2K maxPoolSize:1M incPoolSize:4K" ];

auto fastAbs(float x) @system
{
	uint p = *(cast(uint*) &x);
	p &=0x7fffffff;
	return *(cast(float*) &p);
}

auto renormalize(Range)(Range r)
{
	import std.algorithm;

	auto maximalAbsoluteValue = r
								.map!(a => fastAbs(a))
								.fold!max;
	return r
			.map!(a => a / maximalAbsoluteValue);
}

void main()
{
	import std.algorithm;
	import std.complex;
	import std.math;
	import std.numeric;
	import std.random;
	import std.range;

	import ppmformats;

	auto rng = Random(unpredictableSeed);

	auto signalWithNoise(float x)
	{
		return (3.3 * sin(2*x)) + uniform(-0.3, 0.3, rng);
	}
	
	alias squareWave = delegate(float x) {
	  return
          ((4 * 10) / (PI)) * (sin(x) + (sin(3 * x) / (3 * x)) + (sin(5 * x) /(5 * x)) + (sin(7 * x) / (7 * x)) + (sin(9 * x) / (9 * x)) + (sin(11 * x) / (11 * x))) + (sin(13 * x) / (13 * x)) + (sin(15 * x) / (15 * x)) + (sin(17 * x) / (17 * x));
	};

	auto signal = iota(0, 2_048)
						.map!squareWave;

	auto normalizedFFT(Range)(Range r)
	{
		return 
				r
				  .fft
				  .map!(a => 20.0f * log10(abs(a ^^ 2)))
				  .renormalize;
	}

	auto pd = signal
					.slide(1_024)
					.map!(a => normalizedFFT(a));

	auto img = new P6Image(1024, 1024);

	foreach (i, e; pd.enumerate)
	{
		foreach (j, s; e.enumerate)
		{
			img[i, j] = new RGBColor(
				cast(int) (0.8 * 255.0 * abs(s)),
				0,
				255 - cast(int) (255 * abs(s))
			);
		}
	}
	img.save(`fft.ppm`);		
}