Category: Complex Functions

Blog Entry © Monday, June 15, 2026, by James Pate Williams, Jr. Filtering a Noisy Signal

Blog Entry Monday June 15 2026Download
#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

Blog Entry © Thursday – Saturday, June 11 – 13, 2026, by James Pate Williams, Jr. C/C++ Translation of S. D. Conte and Carl de Boor’s Cooley-Tukey Fast Fourier Transform Algorithm and Related Algorithms

Blog Entry Saturday May 13 2026Download
#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;
}

Blog Entry © Sunday, March 22, 2026, by James Pate Williams, Jr. Mueller’s Method

Blog Entry Sunday March 22 2026Download
/*
* MuellersMethod.c (c) Sunday, July 21, 2024 by
* James Pate Williams, Jr. BA, BS, MSwE, PhD
* Translated from the FORTRAN 77 source code
* found in "Elementary Numerical Analysis: An
* Algorithmic Approach" by S. D. Conte and Carl
* de Boor Originally coded in FORTRAN IV in 1982 then
* into C# in March 2015 Finished Tuesday,
* July 23, 2024 The complex division method is
* from "A Numerical Library in C for Scientists
* and Engineers" by H. T. Lau Chapter 1 
*/

#include <complex.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

//#define DEBUG

static _Lcomplex HornersMethod(_Lcomplex coeff[], _Lcomplex z, int degree) {
	int i = 0;
	_Lcomplex c = coeff[degree];

	for (i = degree; i >= 1; i--) {
		_Lcomplex product = _LCmulcc(c, z);
		c._Val[0] = product._Val[0] + coeff[i - 1]._Val[0];
		c._Val[1] = product._Val[1] + coeff[i - 1]._Val[1];
	}
#ifdef DEBUG
	_Lcomplex sum = { 0 };
	for (i = 0; i <= degree; i++) {
		_Lcomplex term =
			_LCmulcc(cpowl(z, _LCbuild(i, 0.0)), coeff[i]);
		sum._Val[0] += term._Val[0];
		sum._Val[1] += term._Val[1];
	}

	long double delta = fabsl(cabsl(c) - cabsl(sum));

	if (delta > 1.0e-12)
		exit(-5);
#endif
	return c;
}

static void comdiv(
	long double xr, long double xi,
	long double yr, long double yi,
	long double* zr, long double* zi)
{
	long double h, d;

	if (fabs(yi) < fabs(yr)) {
		if (yi == 0.0) {
			*zr = xr / yr;
			*zi = xi / yr;
		}
		else {
			h = yi / yr;
			d = h * yi + yr;
			*zr = (xr + h * xi) / d;
			*zi = (xi - h * xr) / d;
		}
	}
	else {
		h = yr / yi;
		d = h * yr + yi;
		*zr = (xr * h + xi) / d;
		*zi = (xi * h - xr) / d;
	}
}

#ifdef DEBUG
static _Lcomplex MyComplexDivide(_Lcomplex numer, _Lcomplex denom) {
	long double norm2 =
		denom._Val[0] * denom._Val[0] +
		denom._Val[1] * denom._Val[1];
	_Lcomplex result = { 0 };

	result._Val[0] = (
		numer._Val[0] * denom._Val[0] +
		numer._Val[1] * denom._Val[1]) / norm2;
	result._Val[1] = (
		numer._Val[1] * denom._Val[0] -
		numer._Val[0] * denom._Val[1]) / norm2;
	return result;
}
#endif

static _Lcomplex ComplexDivide(_Lcomplex numer, _Lcomplex denom) {
	_Lcomplex result = { 0 };

	comdiv(
		numer._Val[0], numer._Val[1],
		denom._Val[0], denom._Val[1],
		&result._Val[0], &result._Val[1]);
#ifdef DEBUG
	_Lcomplex myResult = MyComplexDivide(numer, denom);
	long double delta = fabsl(cabsl(result) - cabsl(myResult));

	if (delta > 1.0e-12)
		exit(-6);
#endif
	return result;
}

static int Deflate(
	_Lcomplex coeff[], _Lcomplex zero,
	_Lcomplex* fzero, _Lcomplex* fzrdfl,
	_Lcomplex zeros[], int i, int* count,
	int degree) {
	_Lcomplex denom = { 0 };

	(*count)++;

	*fzero = HornersMethod(coeff, zero, degree);
	*fzrdfl = *fzero;

	if (i < 1)
		return 0;

	for (int j = 1; j <= i; j++) {
		denom._Val[0] = zero._Val[0] - zeros[j - 1]._Val[0];
		denom._Val[1] = zero._Val[1] - zeros[j - 1]._Val[1];

		if (cabsl(denom) == 0.0) {
			zeros[i] = _LCmulcr(zero, 1.001);
			return 1;
		}

		else
			*fzrdfl = ComplexDivide(*fzrdfl, denom);
	}

	return 0;
}

static void Mueller(
	_Lcomplex coeff[], _Lcomplex zeros[],
	double epsilon1, double epsilon2,
	int degree, int fnreal, int maxIts, int n, int nPrev) {
	double eps1 = max(epsilon1, 1.0e-12);
	double eps2 = max(epsilon2, 1.0e-20);
	int count = 0, i = 0;
	_Lcomplex c = { 0 };
	_Lcomplex denom = { 0 };
	_Lcomplex divdf1 = { 0 };
	_Lcomplex divdf2 = { 0 };
	_Lcomplex divdf1p = { 0 };
	_Lcomplex fzero = { 0 };
	_Lcomplex fzr = { 0 };
	_Lcomplex fzdfl = { 0 };
	_Lcomplex fzrdfl = { 0 };
	_Lcomplex fzrprv = { 0 };
	_Lcomplex four = _LCbuild(4.0, 0.0);
	_Lcomplex h = { 0 };
	_Lcomplex hprev = { 0 };
	_Lcomplex sqr = { 0 };
	_Lcomplex zero = { 0 };
	_Lcomplex p5 = _LCbuild(0.5, 0.0);
	_Lcomplex zeropp5 = { 0 };
	_Lcomplex zeromp5 = { 0 };
	_Lcomplex diff = { 0 };
	_Lcomplex tadd = { 0 };
	_Lcomplex tmul = { 0 };
	_Lcomplex umul = { 0 };
	_Lcomplex vmul = { 0 };

	for (i = nPrev; i < n; i++) {
		count = 0;

	Label1:

		zero = zeros[i];
		h = p5;

		zeropp5._Val[0] = zero._Val[0] + p5._Val[0];
		zeropp5._Val[1] = zero._Val[1] + p5._Val[1];

		if (Deflate(
			coeff, zeropp5, &fzr, &divdf1p,
			zeros, i, &count, degree))
			goto Label1;

		zeromp5._Val[0] = zero._Val[0] - p5._Val[0];
		zeromp5._Val[1] = zero._Val[1] - p5._Val[1];

		if (Deflate(
			coeff, zeromp5, &fzr, &fzrprv,
			zeros, i, &count, degree))
			goto Label1;

		hprev._Val[0] = -1.0;
		hprev._Val[1] = 0.0;
		diff._Val[0] = fzrprv._Val[0] - divdf1p._Val[0];
		diff._Val[1] = fzrprv._Val[1] - divdf1p._Val[1];
		if (cabsl(hprev) == 0)
			exit(-2);
		divdf1p = ComplexDivide(diff, hprev);

		if (Deflate(
			coeff, zero, &fzr, &fzrdfl,
			zeros, i, &count, degree))
			goto Label1;

	Label2:

		diff._Val[0] = fzrdfl._Val[0] - fzrprv._Val[0];
		diff._Val[1] = fzrdfl._Val[1] - fzrprv._Val[1];
		if (cabsl(h) == 0)
			exit(-3);
		divdf1 = ComplexDivide(diff, h);
		diff._Val[0] = divdf1._Val[0] - divdf1p._Val[0];
		diff._Val[1] = divdf1._Val[1] - divdf1p._Val[1];
		tadd._Val[0] = h._Val[0] + hprev._Val[0];
		tadd._Val[1] = h._Val[1] + hprev._Val[1];
		if (cabsl(tadd) == 0)
			exit(-3);
		divdf2 = ComplexDivide(diff, tadd);
		hprev = h;
		divdf1p = divdf1;
		tmul = _LCmulcc(h, divdf2);
		c._Val[0] = divdf1._Val[0] + tmul._Val[0];
		c._Val[1] = divdf1._Val[1] + tmul._Val[1];
		tmul = _LCmulcc(c, c);
		umul = _LCmulcc(four, fzrdfl);
		vmul = _LCmulcc(umul, divdf2);
		sqr._Val[0] = tmul._Val[0] - vmul._Val[0];
		sqr._Val[1] = tmul._Val[1] - vmul._Val[1];

		if (fnreal && sqr._Val[0] < 0.0)
		{
			sqr._Val[0] = 0.0;
			sqr._Val[1] = 0.0;
		}

		sqr = csqrtl(sqr);

		if ((c._Val[0] * sqr._Val[0] + c._Val[1] * sqr._Val[1]) < 0.0) {
			denom._Val[0] = c._Val[0] - sqr._Val[0];
			denom._Val[1] = c._Val[1] - sqr._Val[1];
		}
		else {
			denom._Val[0] = c._Val[0] + sqr._Val[0];
			denom._Val[1] = c._Val[1] + sqr._Val[1];
		}
		if (cabsl(denom) <= 0.0)
		{
			denom._Val[0] = 1.0;
			denom._Val[1] = 0.0;
		}
		if (cabsl(denom) == 0)
			exit(-4);
		tmul = _LCmulcr(fzrdfl, -2.0);
		h = ComplexDivide(tmul, denom);
		fzrprv = fzrdfl;
		zero._Val[0] = zero._Val[0] + h._Val[0];
		zero._Val[1] = zero._Val[1] + h._Val[1];

		if (count > maxIts)
			goto Label4;

	Label3:

		if (Deflate(
			coeff, zero, &fzr, &fzrdfl,
			zeros, i, &count, degree))
			goto Label1;

		if (cabsl(h) < eps1 * cabsl(zero))
			goto Label4;

		if (max(cabsl(fzr), cabsl(fzdfl)) < eps2)
			goto Label4;

		if (cabsl(fzrdfl) >= 10.0 * cabsl(fzrprv)) {
			h = _LCmulcr(h, 0.5);
			zero._Val[0] = zero._Val[0] - h._Val[0];
			zero._Val[1] = zero._Val[1] - h._Val[1];
			goto Label3;
		}

		else
			goto Label2;

	Label4:

		zeros[i] = zero;
	}
}

int main(void)
{
	double epsilon1 = 1.0e-15;
	double epsilon2 = 1.0e-15;
	int degree = 0, fnreal = 0, i = 0, maxIts = 1000;
	int n = 0, nPrev = 0;
		
	while (1) {
		_Lcomplex* coeff = NULL;
		_Lcomplex* zeros = NULL;

		printf_s("Degree (0 to quit) = ");
		scanf_s("%d", &degree);

		if (degree == 0)
			break;

		n = degree;
		coeff = calloc(degree + 1, sizeof(_Lcomplex));

		if (coeff == NULL)
			exit(-1);

		zeros = calloc(n, sizeof(_Lcomplex));

		if (zeros == NULL)
			exit(-1);

		for (i = degree; i >= 0; i--) {
			printf_s("coefficient[%d].real = ", i);
			scanf_s("%Lf", &coeff[i]._Val[0]);
			printf_s("coefficient[%d].imag = ", i);
			scanf_s("%Lf", &coeff[i]._Val[1]);
		}

		Mueller(
			coeff, zeros, epsilon1,
			epsilon2, degree, fnreal,
			maxIts, n, nPrev);

		printf_s("\n");

		for (i = 0; i < degree; i++) {
			printf_s("zero[%d].real = %17.10e\t", i, zeros[i]._Val[0]);
			printf_s("zero[%d].imag = %17.10e\n", i, zeros[i]._Val[1]);
		}

		printf_s("\n");

		for (i = 0; i < degree; i++) {
			_Lcomplex func = HornersMethod(coeff, zeros[i], degree);

			printf_s("func[%d].real = %17.10e\t", i, func._Val[0]);
			printf_s("func[%d].imag = %17.10e\n", i, func._Val[1]);
		}

		printf_s("\n");

		free(coeff);
		free(zeros);
	}

	return 0;
}

Blog Entry © Sunday, November 9, 2025, by James Pate Williams, Jr. Hydrogenic Wavefunctions, Radial Probability Functions, Distribution Functions, and First Moment Integrals

Blog Entry Sunday November 09 2025Download

Blog Entry (c) Wednesday, August 13, 2025, by James Pate Williams, Jr. Exercises from an Online Textbook

Exercises Szabo OstlundDownload
#include <complex>
#include <vector>

class CmpLinearAlgebra
{
public:
    static void CmpPrintMatrix(
        int m, int n,
        std::vector<std::vector<std::complex<double>>>& Ac);
    static void CmpAddition(
        size_t m, size_t n,
        std::vector<std::vector<std::complex<double>>>& A,
        std::vector<std::vector<std::complex<double>>>& B,
        std::vector<std::vector<std::complex<double>>>& C);
    static void CmpSubtraction(
        size_t m, size_t n,
        std::vector<std::vector<std::complex<double>>>& A,
        std::vector<std::vector<std::complex<double>>>& B,
        std::vector<std::vector<std::complex<double>>>& C);
    static void CmpMultiply(
        size_t m, size_t n, size_t p,
        std::vector<std::vector<std::complex<double>>>& A,
        std::vector<std::vector<std::complex<double>>>& B,
        std::vector<std::vector<std::complex<double>>>& C);
    static void CmpAnticommutator(
        size_t n,
        std::vector<std::vector<std::complex<double>>>& A,
        std::vector<std::vector<std::complex<double>>>& B,
        std::vector<std::vector<std::complex<double>>>& C);
    static void CmpCommutator(
        size_t n,
        std::vector<std::vector<std::complex<double>>>& A,
        std::vector<std::vector<std::complex<double>>>& B,
        std::vector<std::vector<std::complex<double>>>& C);
    static void CmpAdjoint(
        size_t m, size_t n,
        std::vector<std::vector<std::complex<double>>>& Ac,
        std::vector<std::vector<std::complex<double>>>& Ad);
    static std::complex<double> CmpDeterminant(
        bool& failure, int n,
        std::vector<std::vector<std::complex<double>>>& A);
    static bool CmpGaussianElimination(
        int m, int n,
        std::vector<std::vector<std::complex<double>>>& A,
        std::vector<std::complex<double>>& b,
        std::vector<std::complex<double>>& x,
        std::vector<size_t>& pivot);
    static bool CmpGaussianFactor(
        int n, std::vector<std::vector<std::complex<double>>>& M,
        std::vector<size_t>& pivot);
    static bool CmpGaussianSolution(
        int n, std::vector<std::vector<std::complex<double>>>& M,
        std::vector<std::complex<double>>& b,
        std::vector<std::complex<double>>& x,
        std::vector<size_t>& pivot);
    static bool CmpSubstitution(
        int m, int n, std::vector<std::vector<std::complex<double>>>& M,
        std::vector<std::complex<double>>& b,
        std::vector<std::complex<double>>& x,
        std::vector<size_t>& pivot);
    static bool CmpInverse(
        int n, std::vector<std::vector<std::complex<double>>>& M,
        std::vector<std::vector<std::complex<double>>>& Mi);
    static void CmpCharPolyAndAdjoint(
        int n,
        std::vector<std::vector<std::complex<double>>>& C,
        std::vector<std::vector<std::complex<double>>>& I,
        std::vector<std::vector<std::complex<double>>>& M,
        std::vector<std::vector<std::complex<double>>>& adjoint,
        std::vector<std::complex<double>>& a);
    static void CmpMatrixKernel(
        int m, int n,
        std::vector<std::vector<std::complex<double>>>& M,
        std::vector<std::vector<std::complex<double>>>& X,
        size_t& r);
    static void CmpMatrixImage(
        int m, int n,
        std::vector<std::vector<std::complex<double>>>& M,
        std::vector<std::vector<std::complex<double>>>& N,
        std::vector<std::vector<std::complex<double>>>& X,
        int rank);
};

#include <vector>

class DblLinearAlgebra
{
public:
    static void DblPrintMatrix(
        int m, int n, std::vector<std::vector<double>>& A);
    static void DblAddition(
        size_t m, size_t n,
        std::vector<std::vector<double>>& A,
        std::vector<std::vector<double>>& B,
        std::vector<std::vector<double>>& C);
    static void DblSubtraction(
        size_t m, size_t n,
        std::vector<std::vector<double>>& A,
        std::vector<std::vector<double>>& B,
        std::vector<std::vector<double>>& C);
    static void DblMultiply(
        size_t m, size_t n, size_t p,
        std::vector<std::vector<double>>& A,
        std::vector<std::vector<double>>& B,
        std::vector<std::vector<double>>& C);
    static void DblAnticommutator(
        size_t n,
        std::vector<std::vector<double>>& A,
        std::vector<std::vector<double>>& B,
        std::vector<std::vector<double>>& C);
    static void DblCommutator(
        size_t n,
        std::vector<std::vector<double>>& A,
        std::vector<std::vector<double>>& B,
        std::vector<std::vector<double>>& C);
    static double DblDeterminant(
        bool& failure, int n,
        std::vector<std::vector<double>>& A);
    static bool DblGaussianElimination(
        int m, int n, std::vector<std::vector<double>>& A,
        std::vector<double>& b, std::vector<double>& x,
        std::vector<size_t>& pivot);
    static bool DblGaussianFactor(
        int n, std::vector<std::vector<double>>& M,
        std::vector<size_t>& pivot);
    static bool DblGaussianSolution(
        int n, std::vector<std::vector<double>>& M,
        std::vector<double>& b, std::vector<double>& x,
        std::vector<size_t>& pivot);
    static bool DblSubstitution(
        int n, std::vector<std::vector<double>>& M,
        std::vector<double>& b, std::vector<double>& x,
        std::vector<size_t>& pivot);
    static bool DblInverse(
        int n, std::vector<std::vector<double>>& M,
        std::vector<std::vector<double>>& A);
    static void DblCharPolyAndAdjoint(
        int n,
        std::vector<std::vector<double>>& C,
        std::vector<std::vector<double>>& I,
        std::vector<std::vector<double>>& M,
        std::vector<std::vector<double>>& adjoint,
        std::vector<double>& a);
    static void DblMatrixKernel(
        int m, int n,
        std::vector<std::vector<double>>& M,
        std::vector<std::vector<double>>& X,
        size_t& r);
    static void DblMatrixImage(
        int m, int n,
        std::vector<std::vector<double>>& M,
        std::vector<std::vector<double>>& N,
        std::vector<std::vector<double>>& X,
        int rank);
};

// Exercises from "Modern Quantum Chemistry An Introduction to Advanced
// Electronic Structure Theory" by Attila Szabo and Neil S. Ostlund
// https://chemistlibrary.wordpress.com/wp-content/uploads/2015/02/modern-quantum-chemistry.pdf
// Program (c) Tuesday, August 12, 2025 by James Pate Williams, Jr.

#include <complex>
#include <iomanip>
#include <iostream>
#include <vector>
#include "DblLinearAlgebra.h"
#include "CmpLinearAlgebra.h"

int main()
{
	double AArcb[3][3] = { { 2, 3, -1 }, { 4, 4, -3 }, { -2, 3, -1 } };
	double AArso[3][3] = { { 1, 1, 0 }, { 1, 2, 2 }, { 0, 2, -1 } };
	double BBrso[3][3] = { { 1, -1, 1 }, { -1, 0, 0 }, { 1, 0, 1} };
	double BBr[3][3] = { { 1, -1, 1 }, { -1 , 0, 0 }, { 1, 0, 1 } };
	double AAcr[3][3] = { { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, 9 } };
	double AAci[3][3] = { { 1, 1, 2 }, { 3, 0, 1 }, { 0, 2, 4 } };
	double BBcr[3][3] = { { 1, 0, 1 }, { 1 , 1, 0 }, { 0, 1, 1 } };
	double BBci[3][3] = { { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, 9 } };
	int m = 3, n = 3, p = 3;

	std::vector<double> br(3);
	std::vector<size_t> pivot(3);
	std::vector<std::vector<double>> Arcb(3, std::vector<double>(3));
	std::vector<std::vector<double>> Arso(3, std::vector<double>(3));
	std::vector<std::vector<double>> Brso(3, std::vector<double>(3));
	std::vector<std::vector<double>> Br(3, std::vector<double>(3));
	std::vector<std::vector<double>> Cr(3, std::vector<double>(3));
	std::vector<std::vector<double>> Ai(3, std::vector<double>(3));
	std::vector<std::vector<double>> Ari(3, std::vector<double>(3));
	std::vector<std::vector<std::complex<double>>> Ac(3,
		std::vector<std::complex<double>>(3));
	std::vector<std::vector<std::complex<double>>> Bc(3,
		std::vector<std::complex<double>>(3));
	std::vector<std::vector<std::complex<double>>> Cc(3,
		std::vector<std::complex<double>>(3));
	std::vector<std::vector<std::complex<double>>> Dc(3,
		std::vector<std::complex<double>>(3));
	std::vector<std::vector<std::complex<double>>> Ec(3,
		std::vector<std::complex<double>>(3));
	std::vector<std::vector<std::complex<double>>> Fc(3,
		std::vector<std::complex<double>>(3));
	std::vector<std::vector<std::complex<double>>> Gc(3,
		std::vector<std::complex<double>>(3));
	
	for (int i = 0; i < m; i++)
	{
		for (int j = 0; j < p; j++)
		{
			Arcb[i][j] = AArcb[i][j];
			Arso[i][j] = AArso[i][j];
			Brso[i][j] = BBrso[i][j];
			Ac[i][j]._Val[0] = AAcr[i][j];
			Ac[i][j]._Val[1] = AAci[i][j];
		}
	}

	for (int i = 0; i < p; i++)
	{
		for (int j = 0; j < n; j++)
		{
			Br[i][j] = BBr[i][j];
			Bc[i][j]._Val[0] = BBcr[i][j];
			Bc[i][j]._Val[1] = BBci[i][j];
		}
	}
	
	DblLinearAlgebra::DblMultiply(3, 3, 3, Arcb, Br, Cr);
	std::cout << "Ar * Br = Cr Conte & de Boor" << std::endl;
	DblLinearAlgebra::DblPrintMatrix(3, 3, Cr);
	std::cout << std::endl;
	CmpLinearAlgebra::CmpMultiply(3, 3, 3, Ac, Bc, Cc);
	std::cout << "Ac * Bc = Cc" << std::endl;
	CmpLinearAlgebra::CmpPrintMatrix(3, 3, Cc);
	std::cout << std::endl;
	// Exercise 1.2
	std::cout << "Exercise 1.2 page 5 Commutator" << std::endl;
	DblLinearAlgebra::DblCommutator(3, Arso, Brso, Cr);
	DblLinearAlgebra::DblPrintMatrix(3, 3, Cr);
	std::cout << std::endl;
	std::cout << "Exercise 1.2 page 5 Anticommutator" << std::endl;
	DblLinearAlgebra::DblAnticommutator(3, Arso, Brso, Cr);
	DblLinearAlgebra::DblPrintMatrix(3, 3, Cr);
	std::cout << std::endl;
	CmpLinearAlgebra::CmpAdjoint(3, 3, Cc, Dc);
	std::cout << "Exercise 1.3 page 6 Cc adjoint" << std::endl;
	CmpLinearAlgebra::CmpPrintMatrix(3, 3, Dc);
	std::cout << std::endl;
	CmpLinearAlgebra::CmpAdjoint(3, 3, Ac, Ec);
	CmpLinearAlgebra::CmpAdjoint(3, 3, Bc, Fc);
	CmpLinearAlgebra::CmpMultiply(3, 3, 3, Fc, Ec, Gc);
	std::cout << "Exercise 1.3 page 6 Bc adjoint * Ac adjoint" << std::endl;
	CmpLinearAlgebra::CmpPrintMatrix(3, 3, Gc);
	std::cout << std::endl;
	std::cout << "Ar matrix" << std::endl;
	DblLinearAlgebra::DblPrintMatrix(3, 3, Arcb);
	bool inv = DblLinearAlgebra::DblInverse(n, Arcb, Ai);
	std::cout << std::endl;
	std::cout << "Ar Conte & de Boor inverse = " << inv << std::endl;
	DblLinearAlgebra::DblPrintMatrix(3, 3, Ai);
	std::cout << std::endl;
	std::cout << "Ar * Ar inverse" << std::endl;
	DblLinearAlgebra::DblMultiply(3, 3, 3, Arcb, Ai, Ari);
	DblLinearAlgebra::DblPrintMatrix(3, 3, Ari);
	std::cout << std::endl;
	std::cout << "Ac" << std::endl;
	CmpLinearAlgebra::CmpPrintMatrix(3, 3, Ac);
	std::cout << std::endl;
	inv = CmpLinearAlgebra::CmpInverse(3, Ac, Bc);
	std::cout << "Ac inverse = " << inv << std::endl;
	CmpLinearAlgebra::CmpPrintMatrix(3, 3, Bc);
	CmpLinearAlgebra::CmpMultiply(3, 3, 3, Ac, Bc, Cc);
	std::cout << std::endl;
	std::cout << "Ac * AC inverse" << std::endl;
	CmpLinearAlgebra::CmpPrintMatrix(3, 3, Cc);
}

Blog Entry © Sunday, August 10, 2025, First-Order Perturbation Treatment of the Helium Atom by James Pate Williams, Jr.

Blog Entry Sunday August 10 2025Download

Approximation of the Ground-State Total Energy of a Beryllium Atom © Sunday, March 30 to Tuesday April 1, 2025, by James Pate Williams, Jr., BA, BS, Master of Software Engineering, PhD Computer Science

Approximation of the Ground-State Total Energy of a Beryllium AtomDownload

Blog Entry © Thursday, March 27, 2025, by James Pate Williams, Jr., BA, BS, Master of Software Engineering, PhD Lithium (Li, Z = 3) Total Ground-State Energy Numerical Experiments

Lithium Atom ExperimentsDownload

Blog Entry © Tuesday, March 25, 2025, by James Pate Williams, Jr. Hydrogen Radial Wavefunctions and Related Functions

Hydrogen Radial Wavefunctions and Related FunctionsDownload