Numerical Explorations

#pragma once
#include <complex>
#include <vector>

class Transform
{
public:
	static void VandermondeDFT(
		int n,
		std::vector<std::complex<double>>& a,
		std::vector<std::complex<double>>& y);
	static void InverseVandermondeDFT(
		int n,
		std::vector<std::complex<double>>& a,
		std::vector<std::complex<double>>& y);
	static std::vector<std::complex<double>> DFT(
		std::vector<double>& x, std::vector<double>& f);
	static std::vector<double> InverseDFT(
		std::vector<double>& f,
		std::vector<std::complex<double>>& X);
	/*
	 * Reference: "Elementary Numerical Analysis:
	 * An Algorithmic Approach Third Edition" (c)
	 * 1980 by S. D. Conte and Carl de Boor
	 * Section 6.5 pages 268 - 277 and Section 6.6
	 * pages 277 - 283
	 * Input to FFT
	 * Z1, Z2 complex n-vectors
	 * n the length of the vectors
	 * inzee
	 *	= 1 transform in Z1
	 *   = 2 transform in Z2
	 * Constructs the discrete Fourier transform in the Cooley-
	 * Tukey way, but with a twist.
	 */
	static void FFT(
		std::vector<std::complex<double>>& Z1,
		int& after, int& now, int& before, int& inzee,
		std::vector<std::complex<double>>& Z2);
	/*
	 * This computes an in - place complex - to - complex FFT
	 * x and y are the real and imaginary arrays of 2^m points.
	 * dir =  1 gives forward transform
	 * dir = -1 gives reverse transform
	 * see http://astronomy.swin.edu.au/~pbourke/analysis/dft/
	 * Website no longer exists
	 */
	static void FFT(short dir, int m,
		std::vector<double>& x, std::vector<double>& y);
	/*
	 * Reference: "Introduction to Algorithms" by
	 * Thomas H. Cormen, Charles E. Leiserson, and
	 * Ronald L. Rivest, pages 794 - 795
	 */
	static void IterativeFFT(
		std::vector<std::complex<double>>& a,
		std::vector<std::complex<double>>& A);
	/*
	 * Reference: "Introduction to Algorithms" by
	 * Thomas H. Cormen, Charles E. Leiserson, and
	 * Ronald L. Rivest, page 788
	 */
	static std::vector<std::complex<double>> RecursiveFFT(
		std::vector<std::complex<double>>& a);
};

#include "Transform.h"

void Transform::VandermondeDFT(
	int n,
	std::vector<std::complex<double>>& a,
	std::vector<std::complex<double>>& y)
{
	double pi = 4.0 * atan(1.0);
	std::complex<double> z(0.0, 2.0 * pi / n);
	std::complex<double> omegaN = exp(z);
	std::vector<std::vector<std::complex<double>>> V(n);

	for (int k = 0; k < n; k++)
	{
		V[k].resize(n);

		for (int j = 0; j < n; j++)
		{
			V[k][j] = std::pow(omegaN, k * j);
		}
	}

	for (int k = 0; k < n; k++)
	{
		std::complex<double> sum = 0.0;

		for (int j = 0; j < n; j++)
		{
			sum += V[k][j] * a[j];
		}

		y[k] = sum;
	}
}

void Transform::InverseVandermondeDFT(
	int n,
	std::vector<std::complex<double>>& a,
	std::vector<std::complex<double>>& y)
{
	double pi = 4.0 * atan(1.0);
	std::complex<double> nc = { static_cast<double>(n), 0.0 };
	std::complex<double> z(0.0, 2.0 * pi / n);
	std::complex<double> omegaN = exp(z);
	std::vector<std::vector<std::complex<double>>> invV(n);

	for (int k = 0; k < n; k++)
	{
		invV[k].resize(n);

		for (int j = 0; j < n; j++)
		{
			invV[k][j] = std::pow(omegaN, -k * j);
		}
	}

	for (int k = 0; k < n; k++)
	{
		std::complex<double> sum = 0.0;

		for (int j = 0; j < n; j++)
		{
			sum += invV[k][j] * y[j];
		}

		a[k] = sum / nc;
	}
}

std::vector<std::complex<double>> Transform::DFT(
	std::vector<double>& x, std::vector<double>& f)
{
	int length = static_cast<int>(x.size());
	double pi = 4.0 * atan(1.0);
	double pi2oN = 2.0 * pi / length;
	int k, n;
	std::vector<double> X(length);
	std::vector<double> Y(length);
	std::vector<std::complex<double>> Z(length);

	f.resize(length);

	for (k = 0; k < length; k++)
	{
		X[k] = Y[k] = 0;

		for (n = 0; n < length; n++)
		{
			X[k] += x[n] * cos(pi2oN * k * n);
			Y[k] -= x[n] * sin(pi2oN * k * n);
		}

		f[k] = pi2oN * k;
		X[k] /= length;
		Y[k] /= length;
		Z[k] = { X[k], Y[k] };
	}

	return Z;
}

std::vector<double> Transform::InverseDFT(
	std::vector<double>& f,
	std::vector<std::complex<double>>& X)
{
	double imag = 0.0;
	int length = static_cast<int>(X.size());
	std::vector<double> x(length);

	for (int n = 0; n < length; n++)
	{
		imag = x[n] = 0.0;

		for (int k = 0; k < length; k++)
		{
			x[n] += X[k]._Val[0] * cos(f[k] * n)
				- X[k]._Val[1] * sin(f[k] * n);
			imag += X[k]._Val[0] * sin(f[k] * n)
				+ X[k]._Val[1] * cos(f[k] * n);
		}
	}

	return x;
}

static void FFTStep(
	std::vector<std::complex<double>>& Zinp,
	int after, int now, int before,
	std::vector<std::complex<double>>& Zout)
{
	double angle = 0.0, ratio = 0.0;
	double twoPi = 2.0 * 4.0 * atan(1.0);
	int ia = 0, ib = 0, inp = 0, j = 0;
	std::complex<double> arg = 1.0, omega = 0, value = 0;
	angle = twoPi / ((now + 1) * (after + 1));
	omega = std::complex<double>(cos(angle), -sin(angle));
	int address = 1;

	for (int i = 1; i <= now; i++)
	{
		for (int j = 1; j <= after; j++)
		{
			for (int k = 1; k <= before; k++)
			{
				address = i * j * k;

				if (address < Zout.size())
					Zout[address] = { 0.0, 0.0 };
			}
		}
	}

	address = 1;

	for (int j = 1; j <= now; j++)
	{
		for (ia = 1; ia <= after; ia++)
		{
			for (ib = 1; ib <= before; ib++)
			{
				int address = j * ia * ib;

				if (address < Zinp.size())
					value = Zinp[address];

				for (inp = now - 1; inp >= 1; inp--)
				{
					address = ia * ib * inp;

					if (address < Zinp.size())
						value = value * arg + Zinp[address];
				}

				address = ia * j * ib;

				if (address < Zout.size())
					Zout[address] = value;
			}

			arg *= omega;
		}
	}
}

void Transform::FFT(
	std::vector<std::complex<double>>& Z1,
	int& after, int& now, int& before, int& inzee,
	std::vector<std::complex<double>>& Z2)
{
	std::vector<int> prime =
	{ 0, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43,
		47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97 };
	int next = 1, nextmx = 25;

	after = 1;
	before = (int)Z1.size();
	now = 1;

Label10:

	if (before / prime[next] * prime[next] < before)
	{
		next++;

		if (next <= nextmx)
			goto Label10;

		else
		{
			now = before;
			before = 1;
		}
	}

	else
	{
		now = prime[next];
		before /= prime[next];
	}

	if (inzee == 1)
		FFTStep(Z1, after, now, before, Z2);
	else
		FFTStep(Z2, after, now, before, Z1);

	inzee = 3 - inzee;

	if (before == 1)
		return;

	after *= now;
	goto Label10;
}

void Transform::FFT(short dir, int m,
	std::vector<double>& x, std::vector<double>& y)
{
	int n, i, i1, j, k, i2, l, l1, l2;
	double c1, c2, tx, ty, t1, t2, u1, u2, z;

	// Calculate the number of points

	n = 1;

	for (i = 0; i < m; i++)
		n *= 2;

	// Do the bit reversal

	i2 = n >> 1;
	j = 0;

	for (i = 0; i < n - 1; i++)
	{
		if (i < j)
		{
			tx = x[i];
			ty = y[i];
			x[i] = x[j];
			y[i] = y[j];
			x[j] = tx;
			y[j] = ty;
		}

		k = i2;

		while (k <= j)
		{
			j -= k;
			k >>= 1;
		}

		j += k;
	}

	// Compute the FFT

	c1 = -1.0;
	c2 = 0.0;
	l2 = 1;

	for (l = 0; l < m; l++)
	{
		l1 = l2;
		l2 <<= 1;
		u1 = 1.0;
		u2 = 0.0;

		for (j = 0; j < l1; j++)
		{
			for (i = j; i < n; i += l2)
			{
				i1 = i + l1;
				t1 = u1 * x[i1] - u2 * y[i1];
				t2 = u1 * y[i1] + u2 * x[i1];
				x[i1] = x[i] - t1;
				y[i1] = y[i] - t2;
				x[i] += t1;
				y[i] += t2;
			}

			z = u1 * c1 - u2 * c2;
			u2 = u1 * c2 + u2 * c1;
			u1 = z;
		}

		c2 = sqrt((1.0 - c1) / 2.0);

		if (dir == 1)
			c2 = -c2;

		c1 = sqrt((1.0 + c1) / 2.0);
	}

	// Scaling for forward transform

	if (dir == 1)
	{
		for (i = 0; i < n; i++)
		{
			x[i] /= n;
			y[i] /= n;
		}
	}
}

static void FFTBase(
	std::vector<std::complex<double>> a,
	std::vector<std::complex<double>> A)
{
	double pi = 4.0 * atan(1.0);
	int n = static_cast<int>(a.size());

	for (int s = 1; s <= log2(n); s++)
	{
		int m = static_cast<int>(pow(2, s));
		std::complex<double> z(0.0, 2.0 * pi / m);
		std::complex<double> omegaM = exp(z);

		for (int k = 0; k <= n - 1; k += m)
		{
			std::complex<double>  omega = { 1.0, 0.0 };

			for (int j = 0; j <= m / 2 - 1; j++)
			{
				std::complex<double> t = omega * A[k + j + m / 2];
				std::complex<double> u = A[k + j];
				std::complex<double> jc = { static_cast<double>(j), 0.0 };
				A[k + j] = u + jc;
				A[k + j + m / 2] = u - t;
				omega *= omegaM;
			}
		}
	}
}

static int Reverse(int k)
{
	int digits[32] = { 0 }, i = 0;

	while (k > 0)
	{
		int digit = k & 1;
		k >>= 1;
		digits[i++] = digit;
	}

	int result = digits[0];

	for (int j = 1; j < i; j++)
		result = result * 2 + digits[j];

	return result;
}

static void BitReverseCopy(
	std::vector<std::complex<double>>& a,
	std::vector<std::complex<double>>& A)
{
	int n = static_cast<int>(a.size());

	for (int k = 0; k <= n - 1; k++)
		A[Reverse(k)] = a[k];
}

void Transform::IterativeFFT(
	std::vector<std::complex<double>>& a,
	std::vector<std::complex<double>>& A)
{
	BitReverseCopy(a, A);

	double pi = 4.0 * atan(1.0);
	int n = static_cast<int>(a.size());

	for (int s = 1; s <= static_cast<int>(log2(n)); s++)
	{
		int m = static_cast<int>(pow(2.0, s));
		std::complex<double> z(0.0, 2.0 * pi / m);
		std::complex<double> omegaM = exp(z);
		std::complex<double> omega = { 1.0, 0.0 };

		for (int j = 0; j <= m / 2 - 1; j++)
		{
			for (int k = j; k <= n - 1; k += m)
			{
				std::complex<double> t = omega * A[k + m / 2];
				std::complex<double> u = A[k];
				A[k] = u + t;
				A[k + m / 2] = u - t;
				omega *= omegaM;
			}
		}
	}
}

std::vector<std::complex<double>> Transform::RecursiveFFT(
	std::vector<std::complex<double>>& a)
{
	int n = static_cast<int>(a.size());

	if (n == 1)
		return a;

	std::vector<std::complex<double>> a0;
	std::vector<std::complex<double>> a1;
	std::vector<std::complex<double>> y0;
	std::vector<std::complex<double>> y1;

	for (int i = 0; i <= n - 2; i++)
		a0.push_back(a[i]);

	for (int i = 1; i <= n - 1; i++)
		a1.push_back(a[i]);

	y0 = RecursiveFFT(a0);
	y1 = RecursiveFFT(a1);

	double pi = 4.0 * atan(1.0);
	std::complex<double> z(0.0, 2.0 * pi / n);
	std::complex<double> omegaN = exp(z);
	std::complex<double> omega(1.0, 0.0);
	std::vector<std::complex<double>> y(n, 0.0);

	for (int k = 0; k <= n / 2 - 1; k++)
	{
		y[k] = y0[k] + omega * y1[k];
		y[(long long)k + n / 2] = y0[k] - omega * y1[k];
		omega *= omegaN;
	}

	return y;
}

// FilteringNoisySignal.cpp : Defines the entry point for the application.
// Copyright (c) Monday, June 15, 2026 by James Pate Williams, Jr.
// Reference: "Numerical Computation 2: Methods, Software and Analysis"
// (c) 1997 by Christoph W. Ueberhuber pages 52-53.

#include "framework.h"
#include "Resource.h"
#include "FilteringNoisySignal.h"
#include "Transform.h"
#include <float.h>
#include <cmath>
#include <vector>

#define MAX_LOADSTRING 100

typedef struct tagPoint2d
{
    double t, f;
} Point2d, * PPoint2d;

// Global Variables:
HINSTANCE hInst;                                // current instance
WCHAR szTitle[MAX_LOADSTRING];                  // The title bar text
WCHAR szWindowClass[MAX_LOADSTRING];            // the main window class name
char thresholdText[128];                        // threshold buffer
char noisePCText[128];                          // noise % buffer
double threshold;                               // noise threshold
double noisePercent;                            // noise parameter
int Npts = 1024;                                // number of data points
std::vector<double> originalSignal;             // original signal
std::vector<double> perturbedSignalReal;        // perturbed signal real part
std::vector<double> perturbedSignalImag;        // perturbed signal imag part
std::vector<double> recoveredSignalReal;        // recovered signal after filtering
std::vector<double> recoveredSignalImag;        // recovered signal after filtering
std::vector<Point2d> points;                    // plotting data

// Forward declarations of functions included in this code module:
ATOM                MyRegisterClass(HINSTANCE hInstance);
BOOL                InitInstance(HINSTANCE, int);
LRESULT CALLBACK    WndProc(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK    About(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK    InputDialog(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK    DrawGraphDialog(HWND, UINT, WPARAM, LPARAM);

int APIENTRY wWinMain(_In_ HINSTANCE hInstance,
                     _In_opt_ HINSTANCE hPrevInstance,
                     _In_ LPWSTR    lpCmdLine,
                     _In_ int       nCmdShow)
{
    UNREFERENCED_PARAMETER(hPrevInstance);
    UNREFERENCED_PARAMETER(lpCmdLine);

    // TODO: Place code here.

    // Initialize global strings
    LoadStringW(hInstance, IDS_APP_TITLE, szTitle, MAX_LOADSTRING);
    LoadStringW(hInstance, IDC_FILTERINGNOISYSIGNAL, szWindowClass, MAX_LOADSTRING);
    MyRegisterClass(hInstance);

    // Perform application initialization:
    if (!InitInstance (hInstance, nCmdShow))
    {
        return FALSE;
    }

    HACCEL hAccelTable = LoadAccelerators(hInstance, MAKEINTRESOURCE(IDC_FILTERINGNOISYSIGNAL));

    MSG msg;

    // Main message loop:
    while (GetMessage(&msg, nullptr, 0, 0))
    {
        if (!TranslateAccelerator(msg.hwnd, hAccelTable, &msg))
        {
            TranslateMessage(&msg);
            DispatchMessage(&msg);
        }
    }

    return (int) msg.wParam;
}

//
//  FUNCTION: MyRegisterClass()
//
//  PURPOSE: Registers the window class.
//
ATOM MyRegisterClass(HINSTANCE hInstance)
{
    WNDCLASSEXW wcex = { 0 };

    wcex.cbSize = sizeof(WNDCLASSEX);

    wcex.style          = CS_HREDRAW | CS_VREDRAW;
    wcex.lpfnWndProc    = WndProc;
    wcex.cbClsExtra     = 0;
    wcex.cbWndExtra     = 0;
    wcex.hInstance      = hInstance;
    wcex.hIcon          = LoadIcon(hInstance, MAKEINTRESOURCE(IDI_FILTERINGNOISYSIGNAL));
    wcex.hCursor        = LoadCursor(nullptr, IDC_ARROW);
    wcex.hbrBackground  = (HBRUSH)(COLOR_WINDOW+1);
    wcex.lpszMenuName   = MAKEINTRESOURCEW(IDC_FILTERINGNOISYSIGNAL);
    wcex.lpszClassName  = szWindowClass;
    wcex.hIconSm        = LoadIcon(wcex.hInstance, MAKEINTRESOURCE(IDI_SMALL));

    return RegisterClassExW(&wcex);
}

//
//   FUNCTION: InitInstance(HINSTANCE, int)
//
//   PURPOSE: Saves instance handle and creates main window
//
//   COMMENTS:
//
//        In this function, we save the instance handle in a global variable and
//        create and display the main program window.
//
BOOL InitInstance(HINSTANCE hInstance, int nCmdShow)
{
   hInst = hInstance; // Store instance handle in our global variable

   HWND hWnd = CreateWindowW(szWindowClass, szTitle, WS_OVERLAPPEDWINDOW,
      CW_USEDEFAULT, 0, CW_USEDEFAULT, 0, nullptr, nullptr, hInstance, nullptr);

   if (!hWnd)
   {
      return FALSE;
   }

   ShowWindow(hWnd, nCmdShow);
   UpdateWindow(hWnd);

   return TRUE;
}

static double f(double t)
{
    double pi = 4.0 * atan(1.0);
    double arg = 2.0 * pi * t;
    return 2.0 * sin(arg / 500.0) + cos(arg / 200.0) -
        0.5 * sin(arg / 50.0);
}

//
//  FUNCTION: WndProc(HWND, UINT, WPARAM, LPARAM)
//
//  PURPOSE: Processes messages for the main window.
//
//  WM_COMMAND  - process the application menu
//  WM_PAINT    - Paint the main window
//  WM_DESTROY  - post a quit message and return
//
//
LRESULT CALLBACK WndProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam)
{
    switch (message)
    {
    case WM_COMMAND:
        {
            int wmId = LOWORD(wParam);
            // Parse the menu selections:
            switch (wmId)
            {
            case IDM_ABOUT:
                DialogBox(hInst, MAKEINTRESOURCE(IDD_ABOUTBOX), hWnd, About);
                break;
            case IDM_INPUT:
                DialogBox(hInst, MAKEINTRESOURCE(IDD_INPUT_DIALOG), hWnd, InputDialog);
                break;
            case IDM_EXIT:
                DestroyWindow(hWnd);
                break;
            default:
                return DefWindowProc(hWnd, message, wParam, lParam);
            }
        }
        break;
    case WM_PAINT:
        {
            PAINTSTRUCT ps;
            HDC hdc = BeginPaint(hWnd, &ps);
            // TODO: Add any drawing code that uses hdc here...
            EndPaint(hWnd, &ps);
        }
        break;
    case WM_DESTROY:
        PostQuitMessage(0);
        break;
    default:
        return DefWindowProc(hWnd, message, wParam, lParam);
    }
    return 0;
}

// Message handler for about box.
INT_PTR CALLBACK About(HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
    UNREFERENCED_PARAMETER(lParam);
    switch (message)
    {
    case WM_INITDIALOG:
        return (INT_PTR)TRUE;

    case WM_COMMAND:
        if (LOWORD(wParam) == IDOK || LOWORD(wParam) == IDCANCEL)
        {
            EndDialog(hDlg, LOWORD(wParam));
            return (INT_PTR)TRUE;
        }
        break;
    }
    return (INT_PTR)FALSE;
}

static void AddWhiteNoise(
    int n, unsigned int seed)
{
    double signalMin = DBL_MAX;
    double signalMax = DBL_MIN;
    std::vector<double> noise(n, 0.0);
    
    perturbedSignalReal.resize(n, 0.0);
    srand(seed);
    
    for (int i = 0; i < n; i++)
    {
        if (originalSignal[i] < signalMin)
            signalMin = originalSignal[i];
        if (originalSignal[i] > signalMax)
            signalMax = originalSignal[i];
    }

    for (int i = 0; i < n; i++)
        noise[i] = (signalMax - signalMin) * rand() / 
            RAND_MAX + signalMin;

    for (int i = 0; i < n; i++)
    {
        double newSignal = originalSignal[i] +
            noisePercent * noise[i];

        if (newSignal < signalMin)
            newSignal = signalMin;
        if (newSignal > signalMax)
            newSignal = signalMax;

        perturbedSignalReal[i] = newSignal;
    }
}

static void Filter(
    double threshold, int n, unsigned int seed, HWND hDlg)
{
    double pi = 4.0 * atan(1.0);
    double step = 2.0 * pi * 200.0 / n, t = 0;
    points.clear();

    for (int i = 0; i < n; i++)
    {
        double ft = f(t);
        Point2d pt = { t, ft };
        points.push_back(pt);
        t += step;
    }

    DialogBox(hInst, MAKEINTRESOURCE(IDD_DRAW_GRAPH_DIALOG), hDlg,
        DrawGraphDialog);

    originalSignal.resize(n, 0.0);

    for (int i = 0; i < n; i++)
        originalSignal[i] = points[i].f;

    AddWhiteNoise(n, seed);
    points.clear();
    t = 0;

    for (int i = 0; i < n; i++)
    {
        double ft = perturbedSignalReal[i];
        Point2d pt = { t, ft };
        points.push_back(pt);
        t += step;
    }

    DialogBox(hInst, MAKEINTRESOURCE(IDD_DRAW_GRAPH_DIALOG), hDlg,
        DrawGraphDialog);

    perturbedSignalImag.resize(n, 0.0);
    recoveredSignalReal.resize(n, 0.0);
    recoveredSignalImag.resize(n, 0.0);

    Transform::FFT(+1, static_cast<int>(log2(n)),
        perturbedSignalReal, perturbedSignalImag);

    for (int i = 0; i < n; i++)
    {
        double magnitude2 = 
            perturbedSignalReal[i] *
            perturbedSignalReal[i] +
            perturbedSignalImag[i] *
            perturbedSignalImag[i];

        if (magnitude2 > threshold)
        {
            recoveredSignalReal[i] = perturbedSignalReal[i];
            recoveredSignalImag[i] = perturbedSignalImag[i];
        }
    }

    Transform::FFT(+1, static_cast<int>(log2(n)),
        recoveredSignalReal, recoveredSignalImag);
    
    points.clear();
    t = 6 * step;

    for (int i = 6; i < n - 6; i++)
    {
        double ft = recoveredSignalReal[i];
        Point2d pt = { t, ft };
        points.push_back(pt);
        t += step;
    }

    DialogBox(hInst, MAKEINTRESOURCE(IDD_DRAW_GRAPH_DIALOG), hDlg,
        DrawGraphDialog);
}

INT_PTR CALLBACK InputDialog(HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
    UNREFERENCED_PARAMETER(lParam);
    switch (message)
    {
    case WM_INITDIALOG:
        SetDlgItemText(hDlg, IDC_EDIT_THRESHOLD, L"1.0e-2");
        SetDlgItemText(hDlg, IDC_EDIT_NOISE, L"0.25");
        SetDlgItemText(hDlg, IDC_EDIT_N, L"1024");
        return (INT_PTR)TRUE;

    case WM_COMMAND:
        if (LOWORD(wParam) == IDC_BUTTON_DRAW)
        {
            GetDlgItemTextA(hDlg, IDC_EDIT_THRESHOLD, thresholdText, 128);
            threshold = atof(thresholdText);
            GetDlgItemTextA(hDlg, IDC_EDIT_NOISE, noisePCText, 128);
            noisePercent = atof(noisePCText);
            Npts = GetDlgItemInt(hDlg, IDC_EDIT_N, FALSE, FALSE);

            if (LOWORD(wParam) == IDC_BUTTON_DRAW)
            {
                Filter(threshold, Npts, 1, hDlg);
                return (INT_PTR)TRUE;
            }
        }
        if (LOWORD(wParam) == IDCANCEL)
        {
            EndDialog(hDlg, LOWORD(wParam));
            return (INT_PTR)TRUE;
        }
        break;
    }
    return (INT_PTR)FALSE;
}

static void FindMinMax(
    double& tMin, double& tMax,
    double& fMin, double& fMax)
{
    // uses global 2D double points structure

    tMin = fMin = DBL_MAX;
    tMax = fMax = DBL_MIN;

    for (size_t i = 0; i < points.size(); i++)
    {
        Point2d pt = points[i];
        double t = pt.t;
        double f = pt.f;

        if (t < tMin)
            tMin = t;
        if (t > tMax)
            tMax = t;
        if (f < fMin)
            fMin = f;
        if (f > fMax)
            fMax = f;
    }
}

static void DrawFormattedText(HDC hdc, char text[], RECT rect)
{
    // Draw the text with formatting options
    DrawTextA(hdc, text, -1, &rect, DT_SINGLELINE | DT_NOCLIP);
}

INT_PTR CALLBACK DrawGraphDialog(
    HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
    UNREFERENCED_PARAMETER(lParam);
    WCHAR line[256] = { };
    switch (message)
    {
    case WM_INITDIALOG:
        SetWindowText(hDlg, L"Graph Dialog");
        return (INT_PTR)TRUE;
    case WM_COMMAND:
        if (LOWORD(wParam) == IDOK || LOWORD(wParam) == IDCANCEL)
        {
            EndDialog(hDlg, LOWORD(wParam));
            return (INT_PTR)TRUE;
        }
        break;
    case WM_PAINT:
        double h = 0;
        double tMax = 0, tMin = 0, fMax = 0, fMin = 0;
        FindMinMax(tMin, tMax, fMin, fMax);
        float tSpan = (float)(tMax - tMin);
        float fSpan = (float)(fMax - fMin);
        RECT rect = { };
        GetClientRect(hDlg, &rect);
        float width = (float)(rect.right - rect.left + 1);
        float height = (float)(rect.bottom - rect.top - 32 + 1);
        float st0 = 2.0f * width / 16.0f;
        float st1 = 14.0f * width / 16.0f;
        float sf0 = 2.0f * height / 16.0f;
        float sf1 = 14.0f * height / 16.0f;
        float deltaT = tSpan / 8.0f;
        float deltaF = fSpan / 8.0f;
        float tSlope = (st1 - st0) / tSpan;
        float tInter = (float)(st0 - tSlope * tMin);
        float fSlope = (sf0 - sf1) / fSpan;
        float fInter = (float)(sf0 - fSlope * fMax);
        float pt = 0, pf = 0, st = 0, sf = 0;
        PAINTSTRUCT ps;
        POINT wPt = { };
        HDC hdc = BeginPaint(hDlg, &ps);
        int i = 0;
        float t = (float)tMin;
        float f = (float)fMax;
        pt = t;
        pf = f;
        st = tSlope * pt + tInter;
        sf = fSlope * pf + fInter;
        MoveToEx(hdc, (int)st, (int)sf0, &wPt);
        char buffer[128] = { };

        while (i <= 8)
        {
            st = tSlope * t + tInter;
            wPt.x = wPt.y = 0;
            MoveToEx(hdc, (int)st, (int)sf0, &wPt);
            LineTo(hdc, (int)st, (int)sf1);
            sprintf_s(buffer, "%5.4f", t);
            SIZE size = { };
            GetTextExtentPoint32A(
                hdc,
                buffer,
                (int)strlen(buffer),
                &size);
            RECT textRect = { };
            textRect.left = (long)(st - size.cx / 2.0f);
            textRect.right = (long)(st + size.cx / 2.0f);
            textRect.top = (long)sf1;
            textRect.bottom = (long)(sf1 + size.cy / 2.0f);
            DrawFormattedText(hdc, buffer, textRect);
            t += deltaT;
            i++;
        }

        i = 0;
        f = (float)fMin;

        while (i <= 8)
        {
            sf = fSlope * f + fInter;
            wPt.x = wPt.y = 0;
            MoveToEx(hdc, (int)st0, (int)sf, &wPt);
            LineTo(hdc, (int)st, (int)sf);

            if (i != 0)
            {
                sprintf_s(buffer, "%+5.3lf", f);
                SIZE size = { };
                GetTextExtentPoint32A(
                    hdc,
                    buffer,
                    (int)strlen(buffer),
                    &size);
                RECT textRect = { };
                textRect.left = (long)(st0 - size.cx - size.cx / 5.0f);
                textRect.right = (long)(st0 - size.cx / 2.0f);
                textRect.top = (long)(sf - size.cy / 2.0f);
                textRect.bottom = (long)(sf + size.cy / 2.0f);
                DrawFormattedText(hdc, buffer, textRect);
            }

            f += deltaF;
            i++;
        }

        HGDIOBJ bPenNew = NULL;
        HGDIOBJ hPenOld = NULL;

        bPenNew = CreatePen(PS_SOLID, 2, RGB(0, 0, 255));
        hPenOld = SelectObject(hdc, bPenNew);

        pt = (float)points[0].t;
        pf = (float)points[0].f;
        st = tSlope * pt + tInter;
        sf = fSlope * pf + fInter;
        wPt.x = wPt.y = 0;
        MoveToEx(hdc, (int)st, (int)sf, &wPt);

        for (size_t j = 1; j < points.size(); j++)
        {
            pt = (float)points[j].t;
            pf = (float)points[j].f;
            st = tSlope * pt + tInter;
            sf = fSlope * pf + fInter;
            LineTo(hdc, (int)st, (int)sf);
        }

        SelectObject(hdc, hPenOld);
        DeleteObject(bPenNew);

        return (INT_PTR)FALSE;
    }

    return (INT_PTR)FALSE;
}

//Microsoft Visual C++ generated resource script.
//
#include "resource.h"

#define APSTUDIO_READONLY_SYMBOLS
/////////////////////////////////////////////////////////////////////////////
//
// Generated from the TEXTINCLUDE resource.
//
#ifndef APSTUDIO_INVOKED
#include "targetver.h"
#endif
#define APSTUDIO_HIDDEN_SYMBOLS
#include "windows.h"
#undef APSTUDIO_HIDDEN_SYMBOLS
/////////////////////////////////////////////////////////////////////////////
#undef APSTUDIO_READONLY_SYMBOLS

#if !defined(AFX_RESOURCE_DLL) || defined(AFX_TARG_ENU)
LANGUAGE 9, 1

/////////////////////////////////////////////////////////////////////////////
//
// Icon
//

// Icon with lowest ID value placed first to ensure application icon
// remains consistent on all systems.

IDI_FILTERINGNOISYSIGNAL       ICON         "FilteringNoisySignal.ico"
IDI_SMALL               ICON         "small.ico"

/////////////////////////////////////////////////////////////////////////////
//
// Menu
//

IDC_FILTERINGNOISYSIGNAL MENU
BEGIN
    POPUP "&Start"
    BEGIN
        MENUITEM "&Input",               IDM_INPUT
        MENUITEM SEPARATOR
        MENUITEM "E&xit",                IDM_EXIT
    END
    POPUP "&Help"
    BEGIN
        MENUITEM "&About ...",           IDM_ABOUT
    END
END

/////////////////////////////////////////////////////////////////////////////
//
// Accelerator
//

IDC_FILTERINGNOISYSIGNAL ACCELERATORS
BEGIN
    "?",            IDM_ABOUT,              ASCII,  ALT
    "/",            IDM_ABOUT,              ASCII,  ALT
END


/////////////////////////////////////////////////////////////////////////////
//
// Dialog
//

IDD_ABOUTBOX DIALOGEX 0, 0, 170, 62
STYLE DS_SETFONT | DS_MODALFRAME | DS_FIXEDSYS | WS_POPUP | WS_CAPTION | WS_SYSMENU
CAPTION "About FilteringNoisySignal"
FONT 8, "MS Shell Dlg"
BEGIN
    ICON            IDI_FILTERINGNOISYSIGNAL,IDC_STATIC,14,14,21,20
    LTEXT           "FilteringNoisySignal, Version 1.0",IDC_STATIC,42,14,114,8,SS_NOPREFIX
    LTEXT           "Copyright (c) 2026",IDC_STATIC,42,26,114,8
    DEFPUSHBUTTON   "OK",IDOK,113,41,50,14,WS_GROUP
END

IDD_INPUT_DIALOG DIALOGEX 0, 0, 310, 120
STYLE DS_SETFONT | WS_POPUP | WS_VISIBLE | WS_CAPTION | WS_SYSMENU
CAPTION "Input Dialog"
FONT 10, "Courier New", 700
BEGIN
LTEXT       "Threshold", IDC_STATIC, 10, 0, 40, 12
EDITTEXT    IDC_EDIT_THRESHOLD, 55, 0, 40, 14, ES_AUTOHSCROLL
LTEXT       "Noise %", IDC_STATIC, 10, 15, 40, 12
EDITTEXT    IDC_EDIT_NOISE, 55, 15, 40, 14, ES_AUTOHSCROLL
LTEXT       "N", IDC_STATIC, 10, 30, 40, 12
EDITTEXT    IDC_EDIT_N, 55, 30, 40, 14, ES_AUTOHSCROLL
PUSHBUTTON  "Draw", IDC_BUTTON_DRAW, 10, 75, 50, 16
PUSHBUTTON  "Cancel", IDCANCEL, 180, 75, 50, 16
END

IDD_DRAW_GRAPH_DIALOG DIALOGEX 0, 0, 410, 310
STYLE DS_SETFONT | WS_POPUP | WS_VISIBLE | WS_CAPTION | WS_SYSMENU
CAPTION "Draw Graph Dialog"
FONT 10, "Courier New", 700
BEGIN
END

/////////////////////////////////////////////////////////////////////////////
//
// DESIGNINFO
//

#ifdef APSTUDIO_INVOKED
GUIDELINES DESIGNINFO
BEGIN
    IDD_ABOUTBOX, DIALOG
    BEGIN
        LEFTMARGIN, 7
        RIGHTMARGIN, 163
        TOPMARGIN, 7
        BOTTOMMARGIN, 55
    END
END
#endif    // APSTUDIO_INVOKED

#ifdef APSTUDIO_INVOKED
/////////////////////////////////////////////////////////////////////////////
//
// TEXTINCLUDE
//
1 TEXTINCLUDE
BEGIN
    "resource.h\0"
END

2 TEXTINCLUDE
BEGIN
    "#ifndef APSTUDIO_INVOKED\r\n"
    "#include ""targetver.h""\r\n"
    "#endif\r\n"
    "#define APSTUDIO_HIDDEN_SYMBOLS\r\n"
    "#include ""windows.h""\r\n"
    "#undef APSTUDIO_HIDDEN_SYMBOLS\r\n"
    "\0"
END

3 TEXTINCLUDE
BEGIN
    "\r\n"
    "\0"
END

#endif    // APSTUDIO_INVOKED

/////////////////////////////////////////////////////////////////////////////
//
// String Table
//

STRINGTABLE
BEGIN
   IDC_FILTERINGNOISYSIGNAL   "FILTERINGNOISYSIGNAL"
   IDS_APP_TITLE       "FilteringNoisySignal"
END

#endif
/////////////////////////////////////////////////////////////////////////////



#ifndef APSTUDIO_INVOKED
/////////////////////////////////////////////////////////////////////////////
//
// Generated from the TEXTINCLUDE resource.
//

/////////////////////////////////////////////////////////////////////////////
#endif    // not APSTUDIO_INVOKED

//{{NO_DEPENDENCIES}}
// Microsoft Visual C++ generated include file.
// Used by FilteringNoisySignal.rc

#define IDS_APP_TITLE			103

#define IDR_MAINFRAME			128
#define IDD_FILTERINGNOISYSIGNAL_DIALOG	102
#define IDD_ABOUTBOX			103
#define IDM_ABOUT				104
#define IDM_EXIT				105
#define IDI_FILTERINGNOISYSIGNAL			107
#define IDI_SMALL				108
#define IDC_FILTERINGNOISYSIGNAL			109
#define IDC_MYICON				2
#ifndef IDC_STATIC
#define IDC_STATIC				-1
#endif
// Next default values for new objects
//
#ifdef APSTUDIO_INVOKED
#ifndef APSTUDIO_READONLY_SYMBOLS

#define _APS_NO_MFC					130
#define _APS_NEXT_RESOURCE_VALUE	129
#define _APS_NEXT_COMMAND_VALUE		32771
#define _APS_NEXT_CONTROL_VALUE		1000
#define _APS_NEXT_SYMED_VALUE		110
#endif
#endif

#define IDD_INPUT_DIALOG			1000
#define IDC_EDIT_THRESHOLD			1010
#define IDC_EDIT_NOISE				1020
#define IDC_EDIT_N					1030
#define IDC_BUTTON_DRAW				1040
#define IDD_DRAW_GRAPH_DIALOG		2000
#define IDM_INPUT					32771

#pragma once
#include <complex>
#include <vector>

class Transform
{
public:
	static void VandermondeDFT(
		int n,
		std::vector<std::complex<double>>& a,
		std::vector<std::complex<double>>& y);
	static void InverseVandermondeDFT(
		int n,
		std::vector<std::complex<double>>& a,
		std::vector<std::complex<double>>& y);
	static std::vector<std::complex<double>> DFT(
		std::vector<double>& x, std::vector<double>& f);
	static std::vector<double> InverseDFT(
		std::vector<double>& f,
		std::vector<std::complex<double>>& X);
	/*
	 * Reference: "Elementary Numerical Analysis:
	 * An Algorithmic Approach Third Edition" (c)
	 * 1980 by S. D. Conte and Carl de Boor
	 * Section 6.5 pages 268 - 277 and Section 6.6
	 * pages 277 - 283
	 * Input to FFT
	 * Z1, Z2 complex n-vectors
	 * n the length of the vectors
	 * inzee
	 *	= 1 transform in Z1
	 *   = 2 transform in Z2
	 * Constructs the discrete Fourier transform in the Cooley-
	 * Tukey way, but with a twist.
	 */
	static void FFT(
		std::vector<std::complex<double>>& Z1,
		int& after, int& now, int& before, int& inzee,
		std::vector<std::complex<double>>& Z2);
	/*
	 * This computes an in - place complex - to - complex FFT
	 * x and y are the real and imaginary arrays of 2^m points.
	 * dir =  1 gives forward transform
	 * dir = -1 gives reverse transform 
	 * see http://astronomy.swin.edu.au/~pbourke/analysis/dft/
	 */
	static void FFT(short dir, int m,
		std::vector<double>& x, std::vector<double>& y);
	/*
	 * Reference: "Introduction to Algorithms" by
	 * Thomas H. Cormen, Charles E. Leiserson, and
	 * Ronald L. Rivest, pages 794 - 795
	 */
	static void IterativeFFT(
		std::vector<std::complex<double>>& a,
		std::vector<std::complex<double>>& A);
	/*
	 * Reference: "Introduction to Algorithms" by
	 * Thomas H. Cormen, Charles E. Leiserson, and
	 * Ronald L. Rivest, page 788
	 */
	static std::vector<std::complex<double>> RecursiveFFT(
		std::vector<std::complex<double>>& a);
};

#include "Transform.h"

void Transform::VandermondeDFT(
	int n,
	std::vector<std::complex<double>>& a,
	std::vector<std::complex<double>>& y)
{
	double pi = 4.0 * atan(1.0);
	std::complex<double> z(0.0, 2.0 * pi / n);
	std::complex<double> omegaN = exp(z);
	std::vector<std::vector<std::complex<double>>> V(n);

	for (int k = 0; k < n; k++)
	{
		V[k].resize(n);

		for (int j = 0; j < n; j++)
		{
			V[k][j] = std::pow(omegaN, k * j);
		}
	}

	for (int k = 0; k < n; k++)
	{
		std::complex<double> sum = 0.0;

		for (int j = 0; j < n; j++)
		{
			sum += V[k][j] * a[j];
		}

		y[k] = sum;
	}
}

void Transform::InverseVandermondeDFT(
	int n,
	std::vector<std::complex<double>>& a,
	std::vector<std::complex<double>>& y)
{
	double pi = 4.0 * atan(1.0);
	std::complex<double> nc = { static_cast<double>(n), 0.0 };
	std::complex<double> z(0.0, 2.0 * pi / n);
	std::complex<double> omegaN = exp(z);
	std::vector<std::vector<std::complex<double>>> invV(n);

	for (int k = 0; k < n; k++)
	{
		invV[k].resize(n);

		for (int j = 0; j < n; j++)
		{
			invV[k][j] = std::pow(omegaN, -k * j);
		}
	}

	for (int k = 0; k < n; k++)
	{
		std::complex<double> sum = 0.0;

		for (int j = 0; j < n; j++)
		{
			sum += invV[k][j] * y[j];
		}

		a[k] = sum / nc;
	}
}

std::vector<std::complex<double>> Transform::DFT(
	std::vector<double>& x, std::vector<double>& f)
{
	int length = static_cast<int>(x.size());
	double pi = 4.0 * atan(1.0);
	double pi2oN = 2.0 * pi / length;
	int k, n;
	std::vector<double> X(length);
	std::vector<double> Y(length);
	std::vector<std::complex<double>> Z(length);

	f.resize(length);

	for (k = 0; k < length; k++)
	{
		X[k] = Y[k] = 0;

		for (n = 0; n < length; n++)
		{
			X[k] += x[n] * cos(pi2oN * k * n);
			Y[k] -= x[n] * sin(pi2oN * k * n);
		}

		f[k] = pi2oN * k;
		X[k] /= length;
		Y[k] /= length;
		Z[k] = { X[k], Y[k] };
	}

	return Z;
}

std::vector<double> Transform::InverseDFT(
	std::vector<double>& f,
	std::vector<std::complex<double>>& X)
{
	double imag = 0.0;
	int length = static_cast<int>(X.size());
	std::vector<double> x(length);

	for (int n = 0; n < length; n++)
	{
		imag = x[n] = 0.0;

		for (int k = 0; k < length; k++)
		{
			x[n] += X[k]._Val[0] * cos(f[k] * n)
				- X[k]._Val[1] * sin(f[k] * n);
			imag += X[k]._Val[0] * sin(f[k] * n)
				+ X[k]._Val[1] * cos(f[k] * n);
		}
	}

	return x;
}

static void FFTStep(
	std::vector<std::complex<double>>& Zinp,
	int after, int now, int before,
	std::vector<std::complex<double>>& Zout)
{
	double angle = 0.0, ratio = 0.0;
	double twoPi = 2.0 * 4.0 * atan(1.0);
	int ia = 0, ib = 0, inp = 0, j = 0;
	std::complex<double> arg = 1.0, omega = 0, value = 0;
	angle = twoPi / ((now + 1) * (after + 1));
	omega = std::complex<double>(cos(angle), -sin(angle));
	int address = 1;

	for (int i = 1; i <= now; i++)
	{
		for (int j = 1; j <= after; j++)
		{
			for (int k = 1; k <= before; k++)
			{
				address = i * j * k;

				if (address < Zout.size())
					Zout[address] = { 0.0, 0.0 };
			}
		}
	}

	address = 1;

	for (int j = 1; j <= now; j++)
	{
		for (ia = 1; ia <= after; ia++)
		{
			for (ib = 1; ib <= before; ib++)
			{
				int address = j * ia * ib;

				if (address < Zinp.size())
					value = Zinp[address];

				for (inp = now - 1; inp >= 1; inp--)
				{
					address = ia * ib * inp;

					if (address < Zinp.size())
						value = value * arg + Zinp[address];
				}

				address = ia * j * ib;

				if (address < Zout.size())
					Zout[address] = value;
			}

			arg *= omega;
		}
	}
}

void Transform::FFT(
	std::vector<std::complex<double>>& Z1,
	int& after, int& now, int& before, int& inzee,
	std::vector<std::complex<double>>& Z2)
{
	std::vector<int> prime = 
		{ 0, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43,
			47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97 };
	int next = 1, nextmx = 25;

	after = 1;
	before = (int)Z1.size();
	now = 1;

Label10:

	if (before / prime[next] * prime[next] < before)
	{
		next++;

		if (next <= nextmx)
			goto Label10;

		else
		{
			now = before;
			before = 1;
		}
	}

	else
	{
		now = prime[next];
		before /= prime[next];
	}

	if (inzee == 1)
		FFTStep(Z1, after, now, before, Z2);
	else
		FFTStep(Z2, after, now, before, Z1);

	inzee = 3 - inzee;

	if (before == 1)
		return;

	after *= now;
	goto Label10;
}

void Transform::FFT(short dir, int m,
	std::vector<double>& x, std::vector<double>& y)
{
	int n, i, i1, j, k, i2, l, l1, l2;
	double c1, c2, tx, ty, t1, t2, u1, u2, z;

	// Calculate the number of points

	n = 1;

	for (i = 0; i < m; i++)
		n *= 2;

	// Do the bit reversal

	i2 = n >> 1;
	j = 0;

	for (i = 0; i < n - 1; i++)
	{
		if (i < j)
		{
			tx = x[i];
			ty = y[i];
			x[i] = x[j];
			y[i] = y[j];
			x[j] = tx;
			y[j] = ty;
		}

		k = i2;

		while (k <= j)
		{
			j -= k;
			k >>= 1;
		}

		j += k;
	}

	// Compute the FFT

	c1 = -1.0;
	c2 = 0.0;
	l2 = 1;

	for (l = 0; l < m; l++)
	{
		l1 = l2;
		l2 <<= 1;
		u1 = 1.0;
		u2 = 0.0;

		for (j = 0; j < l1; j++)
		{
			for (i = j; i < n; i += l2)
			{
				i1 = i + l1;
				t1 = u1 * x[i1] - u2 * y[i1];
				t2 = u1 * y[i1] + u2 * x[i1];
				x[i1] = x[i] - t1;
				y[i1] = y[i] - t2;
				x[i] += t1;
				y[i] += t2;
			}

			z = u1 * c1 - u2 * c2;
			u2 = u1 * c2 + u2 * c1;
			u1 = z;
		}

		c2 = sqrt((1.0 - c1) / 2.0);

		if (dir == 1)
			c2 = -c2;

		c1 = sqrt((1.0 + c1) / 2.0);
	}

	// Scaling for forward transform

	if (dir == 1)
	{
		for (i = 0; i < n; i++)
		{
			x[i] /= n;
			y[i] /= n;
		}
	}
}

static void FFTBase(
	std::vector<std::complex<double>> a,
	std::vector<std::complex<double>> A)
{
	double pi = 4.0 * atan(1.0);
	int n = static_cast<int>(a.size());

	for (int s = 1; s <= log2(n); s++)
	{
		int m = static_cast<int>(pow(2, s));
		std::complex<double> z(0.0, 2.0 * pi / m);
		std::complex<double> omegaM = exp(z);

		for (int k = 0; k <= n - 1; k += m)
		{
			std::complex<double>  omega = { 1.0, 0.0 };
			
			for (int j = 0; j <= m / 2 - 1; j++)
			{
				std::complex<double> t = omega * A[k + j + m / 2];
				std::complex<double> u = A[k + j];
				std::complex<double> jc = { static_cast<double>(j), 0.0 };
				A[k + j] = u + jc;
				A[k + j + m / 2] = u - t;
				omega *= omegaM;
			}
		}
	}
}

static int Reverse(int k)
{
	int digits[32] = { 0 }, i = 0;

	while (k > 0)
	{
		int digit = k & 1;
		k >>= 1;
		digits[i++] = digit;
	}

	int result = digits[0];

	for (int j = 1; j < i; j++)
		result = result * 2 + digits[j];

	return result;
}

static void BitReverseCopy(
	std::vector<std::complex<double>>& a,
	std::vector<std::complex<double>>& A)
{
	int n = static_cast<int>(a.size());
	
	for (int k = 0; k <= n - 1; k++)
		A[Reverse(k)] = a[k];
}

void Transform::IterativeFFT(
	std::vector<std::complex<double>>& a,
	std::vector<std::complex<double>>& A)
{
	BitReverseCopy(a, A);

	double pi = 4.0 * atan(1.0);
	int n = static_cast<int>(a.size());

	for (int s = 1; s <= static_cast<int>(log2(n)); s++)
	{
		int m = static_cast<int>(pow(2.0, s));
		std::complex<double> z(0.0, 2.0 * pi / m);
		std::complex<double> omegaM = exp(z);
		std::complex<double> omega = { 1.0, 0.0 };

		for (int j = 0; j <= m / 2 - 1; j++)
		{
			for (int k = j; k <= n - 1; k += m)
			{
				std::complex<double> t = omega * A[k + m / 2];
				std::complex<double> u = A[k];
				A[k] = u + t;
				A[k + m / 2] = u - t;
				omega *= omegaM;
			}
		}
	}
}

std::vector<std::complex<double>> Transform::RecursiveFFT(
	std::vector<std::complex<double>>& a)
{
	int n = static_cast<int>(a.size());
	
	if (n == 1)
		return a;

	std::vector<std::complex<double>> a0;
	std::vector<std::complex<double>> a1;
	std::vector<std::complex<double>> y0;
	std::vector<std::complex<double>> y1;
	
	for (int i = 0; i <= n - 2; i++)
		a0.push_back(a[i]);
		
	for (int i = 1; i <= n - 1; i++)
		a1.push_back(a[i]);

	y0 = RecursiveFFT(a0);
	y1 = RecursiveFFT(a1);

	double pi = 4.0 * atan(1.0);
	std::complex<double> z(0.0, 2.0 * pi / n);
	std::complex<double> omegaN = exp(z);
	std::complex<double> omega(1.0, 0.0);
	std::vector<std::complex<double>> y(n, 0.0);

	for (int k = 0; k <= n / 2 - 1; k++)
	{
		y[k] = y0[k] + omega * y1[k];
		y[(long long)k + n / 2] = y0[k] - omega * y1[k];
		omega *= omegaN;
	}

	return y;
}

// CooleyTukeyConsole.cpp : This file contains the 'main' function. Program execution begins and ends there.
//

#include <algorithm>
#include <complex>
#include <iomanip>
#include <iostream>
#include <string>
#include <vector>
#include "Transform.h"

int M = 0, N = 0, Sn = 2048;

static double SimpsonsRule(
    double a, double b, int n,
    double (*f)(double))
{
    double h = (b - a) / n;
    double h2 = 2.0 * h;
    double s = 0.0;
    double t = 0.0;
    double x = a + h;

    for (int i = 1; i < n; i += 2)
    {
        s += f(x);
        x += h2;
    }

    x = a + h2;

    for (int i = 2; i < n; i += 2)
    {
        t += f(x);
        x += h2;
    }

    return h * (f(a) + 4 * s + 2 * t + f(b)) / 3.0;
}

static double f(double x)
{
    return x * x * sin(x);
}

static double AMf(double x)
{
    return f(x) * cos(M * x);
}

static double BMf(double x)
{
    return f(x) * sin(M * x);
}

static double AM(int m)
{
    M = m;
    double pi = 4.0 * atan(1.0), twoPi = pi + pi;
    double integ = SimpsonsRule(-pi, pi, Sn, AMf);
    double value = integ / twoPi;

    if (m == 0)
        return value;
    else
        return 2.0 * value;
}

static double BM(int m)
{
    M = m;
    double pi = 4.0 * atan(1.0), twoPi = pi + pi;
    double integ = SimpsonsRule(-pi, pi, Sn, BMf);
    return 2.0 * integ / twoPi;
}

static void FormatPrint(double x)
{
    if (fabs(x) < 1.0e-12)
        x = 0.0;
    std::cout << std::setw(13);
    std::cout << std::setfill(' ');
    std::cout << std::setprecision(10);
    std::cout << x << '\t';
}

static void FormatPrint(std::complex<double> z)
{
    std::cout << std::setw(13);
    std::cout << std::setfill(' ');
    std::cout << std::setprecision(10);
    std::cout << z << '\t';
}

static void GetCooleyTukeyData(
    int n0, int n1, int n2, int n3,
    std::vector<double>& a,
    std::vector<double>& b,
    std::vector<double>& p,
    std::vector<double>& x,
    std::vector<double>& fx)
{
    double pi = 4.0 * atan(1.0), twoPi = pi + pi;

    x[0] = fx[0] = 0.0;

    for (int i = 1; i <= n0; i++)
    {
        x[i] = twoPi * i / n0;
        fx[i] = f(x[i]);
    }

    a[0] = AM(0);

    for (int m = 1; m <= n0; m++)
    {
        a[m] = AM(m);
        b[m] = BM(m);
    }

    for (int n = 0; n <= n0; n++)
    {
        double asum = 0.0;

        for (int m = 1; m < n0; m++)
            asum += a[m] * cos(m * x[n]);

        double bsum = 0.0;

        for (int m = 1; m <= n0; m++)
            bsum += b[m] * sin(m * x[n]);

        p[n] = a[0] / 2.0 + asum + bsum;
    }
}

static void GetData(
    int n0,
    std::vector<double>& a,
    std::vector<double>& b,
    std::vector<double>& p,
    std::vector<double>& x,
    std::vector<double>& fx)
{
    double pi = 4.0 * atan(1.0), twoPi = pi + pi;

    for (int i = 0; i < n0; i++)
    {
        x[i] = twoPi * i / n0;
        fx[i] = f(x[i]);
    }

    a[0] = AM(0);

    for (int m = 1; m < n0; m++)
    {
        a[m] = AM(m);
        b[m] = BM(m);
    }

    for (int n = 0; n < n0; n++)
    {
        double asum = 0.0;

        for (int m = 1; m < n0; m++)
            asum += a[m] * cos(m * x[n]);

        double bsum = 0.0;

        for (int m = 1; m < n0; m++)
            bsum += b[m] * sin(m * x[n]);

        p[n] = a[0] / 2.0 + asum + bsum;
    }
}

static void TestCooleyTukey(
    int n0, int n1, int n2, int n3)
{
    int n01 = n0 + 1;
    std::vector<double> a(n01), b(n01), p(n01), x(n01), fx(n01);
    GetCooleyTukeyData(n0, n1, n2, n3, a, b, p, x, fx);
    double pi = 4.0 * atan(1.0), twoPi = pi + pi;
    int index1 = 1, inzee = 1;
    std::vector<std::complex<double>> pp(n01);
    std::vector<std::complex<double>> Z1 = { {0, 0} };
    std::vector<std::complex<double>> Z2 = { {0, 0} };

    Z1.resize(n01, 0.0);
    Z2.resize(n01, 0.0);

    for (int i = 1; i <= n0; i++)
        Z1[i] = fx[i];

    int after = 0, now = 0, before = 0;
    Transform::FFT(Z1, after, now, before, inzee, Z2);

    std::cout << "Forward Transform" << std::endl;

    for (int i = 1; index1 < n01 && i <= after; i++)
    {
        for (int j = 1; index1 < n01 && j <= now; j++)
        {
            for (int k = 1; index1 < n01 && k <= before; k++)
            {
                std::cout << std::setw(5);
                std::cout << std::setfill(' ');
                std::cout << index1 << '\t';

                FormatPrint(x[index1]);
                FormatPrint(fx[index1]);
                FormatPrint(Z1[index1]);
                FormatPrint(Z2[index1]);
                std::cout << std::endl;
                index1++;
            }
        }
    }
}

static void TestFFT(int n0)
{
    int m = static_cast<int>(log2(n0));
    std::vector<double> a(n0), b(n0), p(n0), x(n0), fx(n0);
    GetData(n0, a, b, p, x, fx);
    std::vector<double> xx(n0), yy(n0);

    for (int i = 0; i < n0; i++)
        xx[i] = fx[i];

    Transform::FFT(+1, m, xx, yy);
    std::cout << "Forward Transform" << std::endl;

    for (int i = 0; i < n0; i++)
    {
        std::cout << std::setw(5);
        std::cout << std::setfill(' ');
        std::cout << i << '\t';

        FormatPrint(x[i]);
        FormatPrint(fx[i]);
        FormatPrint(p[i]);
        FormatPrint(xx[i]);
        FormatPrint(yy[i]);

        std::cout << std::endl;

        if (i == n0 / 2)
            break;
    }

    Transform::FFT(-1, m, xx, yy);
    std::cout << "Backward Transform" << std::endl;

    for (int i = 0; i < n0; i++)
    {
        std::cout << std::setw(5);
        std::cout << std::setfill(' ');
        std::cout << i << '\t';

        FormatPrint(x[i]);
        FormatPrint(fx[i]);
        FormatPrint(p[i]);
        FormatPrint(xx[i]);
        FormatPrint(yy[i]);

        std::cout << std::endl;

        if (i == n0 / 2)
            break;
    }
}

static void TestIterativeFFT(int n0)
{
    std::vector<double> a(n0), b(n0), p(n0), x(n0), fx(n0);
    GetData(n0, a, b, p, x, fx);
    double pi = 4.0 * atan(1.0), twoPi = pi + pi;
    std::vector<std::complex<double>> AA(n0);
    std::vector<std::complex<double>> aa(n0);
    std::vector<std::complex<double>> yy(n0);

    for (int i = 0; i < n0; i++)
        aa[i] = fx[i];

    Transform::IterativeFFT(aa, AA);

    std::cout << "Forward Transform" << std::endl;

    for (int i = 0; i < n0; i++)
    {
        std::cout << std::setw(5);
        std::cout << std::setfill(' ');
        std::cout << i << '\t';

        FormatPrint(x[i]);
        FormatPrint(fx[i]);
        FormatPrint(AA[i]);

        std::cout << std::endl;
    }
}

static void TestRecursiveFFT(int n0)
{
    std::vector<double> a(n0), b(n0), p(n0), x(n0), fx(n0);
    std::vector<double> xx(n0), yy(n0);
    GetData(n0, a, b, p, x, fx);
    double pi = 4.0 * atan(1.0), twoPi = pi + pi;
    std::vector<std::complex<double>> A = { {0, 0} };
    std::vector<std::complex<double>> Y = { {0, 0} };

    A.resize(n0);

    for (int i = 0; i < n0; i++)
        A[i] = fx[i];

    Y = Transform::RecursiveFFT(A);

    std::cout << "Forward Transform" << std::endl;

    for (int i = 0; i < n0; i++)
    {
        std::cout << std::setw(5);
        std::cout << std::setfill(' ');
        std::cout << i << '\t';

        FormatPrint(x[i]);
        FormatPrint(fx[i]);
        FormatPrint(A[i]);
        FormatPrint(Y[i]);

        std::cout << std::endl;
    }
}

static void TestDFT(int n0)
{
    std::vector<double> a(n0), b(n0), p(n0), x(n0), fx(n0), ff(n0);
    GetData(n0, a, b, p, x, fx);
    std::vector<std::complex<double>> zz = Transform::DFT(fx, ff);

    std::cout << "Forward Transform" << std::endl;

    for (int i = 0; i < n0; i++)
    {
        std::cout << std::setw(5);
        std::cout << std::setfill(' ');
        std::cout << i << '\t';

        FormatPrint(x[i]);
        FormatPrint(fx[i]);
        FormatPrint(ff[i]);
        FormatPrint(zz[i]);
        std::cout << std::endl;
    }

    std::cout << "Inverse Transform" << std::endl;

    std::vector<double> inv = Transform::InverseDFT(ff, zz);

    for (int i = 0; i < n0; i++)
    {
        std::cout << std::setw(5);
        std::cout << std::setfill(' ');
        std::cout << i << '\t';

        FormatPrint(x[i]);
        FormatPrint(fx[i]);
        FormatPrint(ff[i]);
        FormatPrint(inv[i]);
        std::cout << std::endl;
    }
}

static void TestVandermondeDFT(int n0)
{
    std::vector<double> a(n0), b(n0), p(n0), x(n0), fx(n0);
    GetData(n0, a, b, p, x, fx);
    std::vector<std::complex<double>> aa(n0);
    std::vector<std::complex<double>> yy(n0);

    for (int i = 0; i < n0; i++)
        aa[i] = fx[i];

    Transform::VandermondeDFT(n0, aa, yy);

    std::cout << "Vandermonde DFT" << std::endl;

    for (int i = 0; i < n0; i++)
    {
        std::cout << std::setw(5);
        std::cout << std::setfill(' ');
        std::cout << i << '\t';

        FormatPrint(x[i]);
        FormatPrint(fx[i]);
        FormatPrint(p[i]);
        FormatPrint(yy[i]);
        std::cout << std::endl;
    }

    Transform::InverseVandermondeDFT(n0, aa, yy);

    std::cout << "Inverse Vandermonde DFT" << std::endl;

    for (int i = 0; i < n0; i++)
    {
        std::cout << std::setw(5);
        std::cout << std::setfill(' ');
        std::cout << i << '\t';

        FormatPrint(x[i]);
        FormatPrint(fx[i]);
        FormatPrint(p[i]);
        FormatPrint(aa[i]);
        std::cout << std::endl;
    }
}

static int Horner(char line[])
{
    int length = static_cast<int>(strlen(line));
    int sum = line[0] - '0';

    for (int i = 1; i < length; i++)
        sum = sum * 10 + line[i] - '0';
    
    return sum;
}

int main()
{
    char line[128];

    while (true)
    {
        std::cout << "== Menu ==" << std::endl;
        std::cout << "1 Cooley-Tukey" << std::endl;
        std::cout << "2 FFT" << std::endl;
        std::cout << "3 Iterative FFT" << std::endl;
        std::cout << "4 Recursive FFT" << std::endl;
        std::cout << "5 DFT" << std::endl;
        std::cout << "6 Vandermonde DFT" << std::endl;
        std::cout << "7 Exit" << std::endl;
        std::cout << "Option 1 - 7 = ";

        std::cin.getline(line, 128);
        int option = Horner(line);

        if (option == 7)
            break;

        if (option < 1 || option > 7)
        {
            std::cout << "Unknown Option Number" << std::endl;
            continue;
        }

        if (option == 1)
        {
            int n0 = 0, n1 = 0, n2 = 0, n3 = 0;

            std::cout << "n1 = ";
            std::cin.getline(line, 128);
            n1 = Horner(line);
            
            std::cout << "n2 = ";
            std::cin.getline(line, 128);
            n2 = Horner(line);

            std::cout << "n3 = ";
            std::cin.getline(line, 128);
            n3 = Horner(line);

            n0 = n1 * n2 * n3;
            TestCooleyTukey(n0, n1, n2, n3);
        }

        else if (option == 2)
        {
            std::cout << "n0 = ";
            std::cin.getline(line, 128);
            int n0 = Horner(line);

            if (n0 % 2 != 0)
            {
                std::cout << "n0 must be a power of 2";
                std::cout << std::endl;
                continue;
            }

            TestFFT(n0);
        }

        else if (option == 3)
        {
            std::cout << "n0 = ";
            std::cin.getline(line, 128);
            int n0 = Horner(line);

            if (n0 % 2 != 0)
            {
                std::cout << "n0 must be a power of 2";
                std::cout << std::endl;
                continue;
            }

            TestIterativeFFT(n0);
        }

        else if (option == 4)
        {
            std::cout << "n0 = ";
            std::cin.getline(line, 128);
            int n0 = Horner(line);

            if (n0 % 2 != 0)
            {
                std::cout << "n0 must be a power of 2";
                std::cout << std::endl;
                continue;
            }

            TestRecursiveFFT(n0);
        }

        else if (option == 5)
        {
            std::cout << "n0 = ";
            std::cin.getline(line, 128);
            int n0 = Horner(line);

            if (n0 % 2 != 0)
            {
                std::cout << "n0 must be a power of 2";
                std::cout << std::endl;
                continue;
            }

            TestDFT(n0);
        }

        else if (option == 6)
        {
            std::cout << "n0 = ";
            std::cin.getline(line, 128);
            int n0 = Horner(line);

            if (n0 % 2 != 0)
            {
                std::cout << "n0 must be a power of 2";
                std::cout << std::endl;
                continue;
            }

            TestVandermondeDFT(n0);
        }
    }

    return 0;
}