Простая реализация преобразования Фурье

Многие программисты используют в своих приложения алгоритмы цифровой обработки сигналов, наиболее распространенным из которых служит быстрое преобразование Фурье (Fast Fourier Transform, FFT).

Несмотря на то, что в стандартной библиотеке D есть шаблон fft (в модуле std.numeric), тем не менее иногда возникает потребность в самописном модуле, который проще контролировать и расширять под свои нужды.
Чтобы удовлетворить эту потребность, а также чтобы показать классическую реализацию FFT, я создал небольшой класс SpecialFFT.

Данный класс выглядит так:

module sfft;

import std.complex;
import std.math;

// псевдоним для комплексных чисел двойнй точности
alias ComplexDouble = Complex!double;

class SpecialFFT
{
	public static
	{
		// поворачивающий множитель
		ComplexDouble twiddleFactor(int k, int N)
		{
			if (k % N == 0)
			{
				return complex(1.0);
			}
			else
			{
				double argument = - 2.0 * PI * (k / N);
				return complex(
					cos(argument),
					sin(argument)
					);
			}
		}

		// быстрое преобразование Фурье (для комплексных чисел)
		ComplexDouble[] fft(ComplexDouble[] x)
		{
			auto N = cast(int) x.length;
			ComplexDouble[] fftResult;
			
			// определяем преобразование Фурье для массива из двух данных
			if (N == 2)
			{
				fftResult = new ComplexDouble[2];
				fftResult[0] = (x[0] + x[1]);
				fftResult[1] = (x[0] - x[1]);
			}
			else
			{
				// собираем преобразование рекурсивно для массивов длины большей чем 2
				ComplexDouble[] 
							evenPart = new ComplexDouble[N / 2], 
							oddPart =  new ComplexDouble[N / 2];

				for (int i = 0; i < N / 2; i++)
				{
					evenPart[i] = x[2 * i];
					oddPart[i] = x[2 * i + 1];
				}

				ComplexDouble[] 
							fftEven = fft(evenPart),
							fftOdd = fft(oddPart);

				fftResult = new ComplexDouble[N];

				for (int i = 0; i < N / 2; i++)
				{
					fftResult[i] = fftEven[i] + twiddleFactor(i, N) * fftOdd[i];
					fftResult[i + N / 2] = fftEven[i] - twiddleFactor(i, N) * fftOdd[i];
				}
			}
			return fftResult;
		}

		// центрирование значений (спектральная составляющая при нулевой частоте будет в центре массива)
		ComplexDouble[] center(ComplexDouble[] x)
		{
			int N = cast(int) x.length;
			ComplexDouble[] centerResult = new ComplexDouble[N];

			for (int i = 0; i < N / 2; i++)
			{
				centerResult[i] = x[(N / 2) + i];
				centerResult[(N / 2) + i] = x[i];
			}

			return centerResult;
		}

		// оконный вариант преобразования Фурье
		ComplexDouble[] takeWindow(ComplexDouble[] x, int start, int size)
		{
			ComplexDouble[] window = new ComplexDouble[size];

			for (int i = 0; i < size; i++)
			{
				window[i] = ComplexDouble(x[i + start].re, x[i + start].im);
			}

			return window;
		}
	}
}

В данный класс входит само преобразование на комплексных числах, а также ряд необходимых функций, таких как: взятие произвольного окна от преобразования и группировка значений полученного спектра вокруг нулевой частоты (центрирование спектра).