Blog Entry © Monday, June 15, 2026, by James Pate Williams, Jr. Filtering a Noisy Signal
#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
#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 © Monday to Tuesday, June 8 – 9, 2026, by James Pate Williams, Jr., DampedNewton’s Method for a System of Equations
Blog Entry © Saturday – Tuesday, May 30 – June 2, 2026, by James Pate Williams, Jr., theMicrosoft Bing Copilot, the M365 Copilot Partial Reproduction of Figures 8 and 9 fromChapter 2 of Quantum Mechanics Third Edition by Leonard I. Schiff
class Figure
{
public:
static void ComputeFigure(
bool eight,
double& xi2, double& eta2,
std::vector<double>& radius,
std::vector<double>& xi,
std::vector<double>& eta,
std::vector<double>& vix,
std::vector<double>& viy,
std::vector<double>& energy1,
double (*f)(double),
double (*g)(double));
};
#include "Figure.h"
#include <algorithm>
#include <cmath>
#include <iterator>
#include <set>
#include <vector>
// Function to compute intersection of two sorted containers of doubles with tolerance
template <typename InputIt1, typename InputIt2, typename OutputIt>
void set_intersection_with_tolerance(InputIt1 first1, InputIt1 last1,
InputIt2 first2, InputIt2 last2,
OutputIt d_first, double tolerance) {
while (first1 != last1 && first2 != last2) {
double a = *first1;
double b = *first2;
if (std::fabs(a - b) <= tolerance) {
// Values are considered equal within tolerance
*d_first++ = /*a; // or */ (a + b) / 2.0; //if you want averaged value
++first1;
++first2;
}
else if (a < b - tolerance) {
++first1;
}
else {
++first2;
}
}
}
static double f(double xi)
{
return xi * tan(xi);
}
static double g(double xi)
{
return -xi * 1.0 / tan(xi);
}
static void Intersection(
bool eight, std::vector<double>& intersection)
{
double del = 0.0, radius2 = 0.0;
double tolerance = 0.0;
double (*h)(double);
int ilimit = 0;
std::set<double> eta, rad;
if (eight)
{
ilimit = 10000;
radius2 = 1.0;
tolerance = 1.0e-1;
h = f;
}
else
{
ilimit = 10000;
radius2 = 4.0;
tolerance = 1.0e-1;
h = g;
}
del = radius2 / ilimit;
intersection.clear();
for (size_t i = 0; i < ilimit; i++)
{
double xi0 = i * del;
if (eight && xi0 >= 0.65)
{
double et = h(xi0);
double r2 = xi0 * xi0 + et * et;
if (fabs(r2 - radius2) < tolerance)
{
eta.insert(et);
rad.insert(r2);
}
}
else if (!eight && xi0 >= 1.8125)
{
double et = h(xi0);
double r2 = xi0 * xi0 + et * et;
if (fabs(r2 - radius2) < tolerance)
{
eta.insert(et);
rad.insert(r2);
}
}
}
if (eight)
{
size_t count = 0, index = 0;
for (double val : eta)
{
if (index == eta.size() - 1)
{
intersection.push_back(val);
break;
}
count++;
index++;
}
}
else
{
int count = 24, index = 0;
for (double val : eta)
{
if (index == eta.size() - count)
{
intersection.push_back(val);
break;
}
index++;
}
}
/*set_intersection_with_tolerance(
eta.begin(), eta.end(),
rad.begin(), rad.end(),
std::back_inserter(intersection),
tolerance);*/
}
void Figure::ComputeFigure(
bool eight,
double& xi2, double& eta2,
std::vector<double>& radius,
std::vector<double>& xi,
std::vector<double>& eta,
std::vector<double>& vix,
std::vector<double>& viy,
std::vector<double>& intersection,
double (*f)(double),
double (*g)(double))
{
double xi0 = 0.0;
double xi1 = 3.5;
double eta0 = 0.0, eta1 = 0.0;
double (*h)(double);
xi.clear();
eta.clear();
if (eight)
h = f;
else
{
xi0 = 1.5;
h = g;
}
eta0 = h(xi0);
eta1 = h(xi1);
double deltaXi = (xi1 - xi0) / 10000.0;
int count = 0;
vix.clear();
viy.clear();
for (int j = 0; j <= 10000; j++)
{
double x = j * deltaXi;
double hx = h(x);
double vx = x * x + hx * hx;
if (eight && x >= xi0 && x <= xi1)
{
if (count == 0 && x >= xi0)
{
xi.push_back(x);
eta.push_back(hx);
count = 1;
}
else
{
xi.push_back(x);
eta.push_back(hx);
}
for (int k = 0; k < 2; k++)
{
double r = radius[k];
vix.push_back(xi[xi.size() - 1]);
viy.push_back(r * r);
}
}
else if (!eight && x >= xi0 && x <= xi1)
{
if (count == 0 && x >= xi0)
{
xi.push_back(x);
eta.push_back(hx);
count = 1;
}
else if (count == 1)
{
xi.push_back(x);
eta.push_back(hx);
}
if (xi.size() > 0)
{
for (int k = 0; k < 1; k++)
{
double r = radius[k];
vix.push_back(xi[xi.size() - 1]);
viy.push_back(r * r);
}
}
}
if (eight && xi[j] >= 0.0 && eta[j] >= 3.5)
break;
else if (!eight &&
xi.size() > 0 &&
xi[xi.size() - 1] > 1.5 &&
eta.size() > 0 &&
eta[eta.size() - 1] >= 3.5)
break;
}
Intersection(eight, intersection);
}
// SchiffChapter2Fig8and9.cpp : Defines the entry point for the application.
//
#include "framework.h"
#include "SchiffChapter2Fig8and9.h"
#include <time.h>
#include <algorithm>
#include <vector>
#include "Figure.h"
#define MAX_LOADSTRING 100
typedef struct tagPoint2d
{
double x, y;
} 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
WCHAR line[8192], text[8192]; // wide character buffers
WCHAR title[65536]; // window title
bool eight; // true = plot figure 8
double xi2, eta2; // energy cordinates
std::vector<Point2d> points; // plotting 2d points
std::vector<double> radius; // radius vector { 1.0, 2.0 }
std::vector<double> xi; // Greek x-coordinate
std::vector<double> eta; // Greek y-coordinate
std::vector<double> V0; // potential of the well
std::vector<double> vix; // potential x-coordinate
std::vector<double> viy; // potential y-coordinate
std::vector<double> energy1; // energy eigenvalues
// 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 DrawEtaDialog(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_SCHIFFCHAPTER2FIG8AND9, szWindowClass, MAX_LOADSTRING);
MyRegisterClass(hInstance);
// Perform application initialization:
if (!InitInstance (hInstance, nCmdShow))
{
return FALSE;
}
HACCEL hAccelTable = LoadAccelerators(hInstance, MAKEINTRESOURCE(IDC_SCHIFFCHAPTER2FIG8AND9));
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_SCHIFFCHAPTER2FIG8AND9));
wcex.hCursor = LoadCursor(nullptr, IDC_ARROW);
wcex.hbrBackground = (HBRUSH)(COLOR_WINDOW+1);
wcex.lpszMenuName = MAKEINTRESOURCEW(IDC_SCHIFFCHAPTER2FIG8AND9);
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 fx(double x)
{
return x * tan(x);
}
static double gx(double x)
{
return -x * 1.0 / tan(x);
}
//
// 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_FIGURE8:
{
eight = true;
text[0] = L'\0';
clock_t clock0 = clock();
radius.push_back(1.0);
radius.push_back(2.0);
Figure::ComputeFigure(
eight, xi2, eta2, radius, xi, eta, vix, viy, energy1, fx, gx);
clock_t clock1 = clock() - clock0;
double runtime = (double)clock1 / CLOCKS_PER_SEC;
swprintf_s(line, 8192, L"Runtime in Seconds = %lf\r\n", runtime);
wcscat_s(text, 8192, line);
points.clear();
for (size_t j = 0; j < xi.size(); j++)
{
Point2d pt = { 0 };
pt.x = xi[j];
pt.y = eta[j];
points.push_back(pt);
}
for (size_t i = 0; i < energy1.size(); i++)
{
swprintf_s(line, 8192,
L"E[%zu] = %lf\r\n", i, energy1[i]);
wcscat_s(text, 8192, line);
}
MessageBox(hWnd, text, L"Energy Information",
MB_OK | MB_ICONINFORMATION);
DialogBox(hInst, MAKEINTRESOURCE(IDD_DRAW_ETA_DIALOG), hWnd, DrawEtaDialog);
break;
}
case IDM_FIGURE9:
{
eight = false;
text[0] = L'\0';
clock_t clock0 = clock();
radius.push_back(2.0);
Figure::ComputeFigure(
eight, xi2, eta2, radius, xi, eta, vix, viy, energy1, fx, gx);
clock_t clock1 = clock() - clock0;
double runtime = (double)clock1 / CLOCKS_PER_SEC;
swprintf_s(line, 8192, L"Runtime in Seconds = %lf\r\n", runtime);
wcscat_s(text, 8192, line);
points.clear();
for (size_t j = 0; j < xi.size(); j++)
{
Point2d pt = { 0 };
pt.x = xi[j];
pt.y = eta[j];
points.push_back(pt);
}
for (size_t i = 0; i < energy1.size(); i++)
{
swprintf_s(line, 8192,
L"E[%zu] = %lf\r\n", i, energy1[i]);
wcscat_s(text, 8192, line);
}
MessageBox(hWnd, text, L"Energy Information",
MB_OK | MB_ICONINFORMATION);
DialogBox(hInst, MAKEINTRESOURCE(IDD_DRAW_ETA_DIALOG), hWnd, DrawEtaDialog);
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 FindMinMax(
std::vector<Point2d>& points,
double& xMin, double& xMax,
double& yMin, double& yMax)
{
// uses global 2D double points structure
xMin = yMin = DBL_MAX;
xMax = yMax = DBL_MIN;
for (size_t i = 0; i < points.size(); i++)
{
Point2d pt = points[i];
double x = pt.x;
double y = pt.y;
if (x < xMin)
xMin = x;
if (x > xMax)
xMax = x;
if (y < yMin)
yMin = y;
if (y > yMax)
yMax = y;
}
}
static void DrawFormattedText(HDC hdc, char text[], RECT rect)
{
// Draw the text with formatting options
DrawTextA(hdc, text, -1, &rect, DT_SINGLELINE | DT_NOCLIP);
}
static void DrawQuarterCircleArc(
HDC hdc,
float xSlope, float ySlope,
float xInter, float yInter,
float radius, bool topToRight)
{
auto mapX = [&](float x)
{
return (int)lroundf(xSlope * x + xInter);
};
auto mapY = [&](float y)
{
return (int)lroundf(ySlope * y + yInter);
};
int x1 = mapX(-radius);
int y1 = mapY(+radius);
int x2 = mapX(+radius);
int y2 = mapY(-radius);
int left = min(x1, x2);
int right = max(x1, x2);
int top = min(y1, y2);
int bottom = max(y1, y2);
int xTop = mapX(0.0f);
int yTop = mapY(radius);
int xRight = mapX(radius);
int yRight = mapY(0.0f);
if (topToRight)
{
Arc(hdc, left, top, right, bottom,
xTop, yTop, xRight, yRight);
}
else
{
Arc(hdc, left, top, right, bottom,
xRight, yRight, xTop, yTop);
}
}
INT_PTR CALLBACK DrawEtaDialog(
HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
UNREFERENCED_PARAMETER(lParam);
switch (message)
{
case WM_INITDIALOG:
SetWindowText(hDlg, title);
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, pi = 0, plm = 0, theta = 0;
double xMax = 0, xMin = 0, yMax = 0, yMin = 0;
FindMinMax(points, xMin, xMax, yMin, yMax);
float xSpan = (float)(xMax - xMin);
float ySpan = (float)(yMax - yMin);
RECT rect = { };
GetClientRect(hDlg, &rect);
float width = (float)(rect.right - rect.left + 1);
float height = (float)(rect.bottom - rect.top - 32 + 1);
float sx0 = 2.0f * width / 16.0f;
float sx1 = 14.0f * width / 16.0f;
float sy0 = 2.0f * height / 16.0f;
float sy1 = 14.0f * height / 16.0f;
float deltaX = xSpan / 8.0f;
float deltaY = ySpan / 8.0f;
float xSlope, xInter, ySlope, yInter;
xSlope = (sx1 - sx0) / xSpan;
xInter = (float)(sx0 - xSlope * xMin);
ySlope = (sy0 - sy1) / ySpan;
yInter = (float)(sy0 - ySlope * yMax);
float px = 0, py = 0, sx = 0, sy = 0;
float vx = 0, vy = 0;
PAINTSTRUCT ps;
POINT wPt = { };
HDC hdc = BeginPaint(hDlg, &ps);
int i = 0;
float x = (float)xMin;
float y = (float)yMax;
px = x;
py = y;
sx = xSlope * px + xInter;
sy = ySlope * py + yInter;
MoveToEx(hdc, (int)sx, (int)sy0, &wPt);
char buffer[128] = { };
while (i <= 8)
{
sx = xSlope * x + xInter;
wPt.x = wPt.y = 0;
MoveToEx(hdc, (int)sx, (int)sy0, &wPt);
LineTo(hdc, (int)sx, (int)sy1);
sprintf_s(buffer, "%5.4lf", x);
SIZE size = { };
GetTextExtentPoint32A(
hdc,
buffer,
(int)strlen(buffer),
&size);
RECT textRect = { };
textRect.left = (long)(sx - size.cx / 2.0f);
textRect.right = (long)(sx + size.cx / 2.0f);
textRect.top = (long)sy1;
textRect.bottom = (long)(sy1 + size.cy / 2.0f);
DrawFormattedText(hdc, buffer, textRect);
x += deltaX;
i++;
}
i = 0;
y = (float)yMin;
while (i <= 8)
{
sy = ySlope * y + yInter;
wPt.x = wPt.y = 0;
MoveToEx(hdc, (int)sx0, (int)sy, &wPt);
LineTo(hdc, (int)sx, (int)sy);
if (i != 0)
{
sprintf_s(buffer, "%5.3lf", y);
SIZE size = { };
GetTextExtentPoint32A(
hdc,
buffer,
(int)strlen(buffer),
&size);
RECT textRect = { };
textRect.left = (long)(sx0 - size.cx - size.cx / 5.0f);
textRect.right = (long)(sx0 - size.cx / 2.0f);
textRect.top = (long)(sy - size.cy / 2.0f);
textRect.bottom = (long)(sy + size.cy / 2.0f);
DrawFormattedText(hdc, buffer, textRect);
}
y += deltaY;
i++;
}
HGDIOBJ bPenNew = NULL;
HGDIOBJ hPenOld = NULL;
bPenNew = CreatePen(PS_SOLID, 2, RGB(0, 0, 255));
hPenOld = SelectObject(hdc, bPenNew);
HRGN hrgn = CreateRectRgn((int)sx0, (int)sy0, (int)sx1, (int)sy1);
// Select the clipping region into the DC
if (SelectClipRgn(hdc, hrgn) == ERROR) {
MessageBox(hDlg, L"Failed to select clip region",
L"Error", MB_ICONERROR);
return (INT_PTR)FALSE;
}
SelectClipRgn(hdc, hrgn);
px = (float)points[0].x;
py = (float)points[0].y;
sx = xSlope * px + xInter;
sy = ySlope * py + yInter;
wPt.x = wPt.y = 0;
MoveToEx(hdc, (int)sx, (int)sy, &wPt);
for (size_t j = 1; j < points.size(); j++)
{
px = (float)points[j].x;
py = (float)points[j].y;
sx = xSlope * px + xInter;
sy = ySlope * py + yInter;
LineTo(hdc, (int)sx, (int)sy);
}
float radius = 0.0f;
if (eight)
radius = 1.0f;
else
radius = 2.0f;
DrawQuarterCircleArc(
hdc, xSlope, ySlope,
xInter, yInter, radius, false);
if (eight)
{
radius = 2.0f;
DrawQuarterCircleArc(
hdc, xSlope, ySlope,
xInter, yInter, radius, false);
}
SelectObject(hdc, hPenOld);
DeleteObject(bPenNew);
return (INT_PTR)FALSE;
}
return (INT_PTR)FALSE;
}
Blog Entry © Thursday, May 28, 2026, by James Pate Williams, Jr. and Microsoft’s Copilot Solution of the Potential Equation in Rectangle using Fixed Point Iteration in Python
# NOTE:
# This implementation prioritizes clarity and correctness over optimization.
# Further performance improvements can be made if needed.
# (c) May 26, 2026 by James Pate Willims, Jr.
# I had some help from the Microsoft Copilot
# to calculate runtimes and define matrices
# Computes the potential in a rectangle
# Reference: "Boundary Value Problems
# Second Edition" by David L. Powers
# See pages 179 to 182 for the analytic
# solution of this Laplace Equation
# Stand alone application using
# Microsoft Visual Studio 2022
# Community Version
import math
import time
xi = yi = 10
u = [[0.0 for _ in range(xi + 2)] for _ in range(yi + 2)]
v = [[0.0 for _ in range(xi + 2)] for _ in range(yi + 2)]
def ComputeBoundaryValues(x, y):
if x == 0:
return 0
if x == 1:
return 0
if y == 0 or y == 1:
if x > 0.0 and x < 0.5:
return 2.0 * x
elif x >= 0.5 and x < 1.0:
return 2.0 - 2.0 * x
return 0.0
def ComputeParams(its, norm, params):
params['iterations'] = its
params['norm'] = norm
def Compute(h, k, xi, yi, maxIts, params):
# Use a simple fixed-point iteration to
# compute an approximate solution
for i in range(0, xi + 1):
for j in range(0, yi + 1):
u[i][j] = ComputeBoundaryValues(i * h, j * k)
for its in range(1, maxIts + 1):
for i in range(1, xi):
for j in range(1, yi):
u[i][j] = 0.25 * (u[i + 1][j] + u[i - 1][j] + u[i][j + 1] + u[i][j - 1]);
norm = 0
for i in range(0, xi + 1):
for j in range(0, yi + 1):
norm += math.fabs(u[i][j] * u[i][j])
norm = math.sqrt(norm)
params['iterations'] = its
params['norm'] = norm
def f(x, y):
# Analytic solution series expansion n = 1 to 100
sum = 0.0
for n in range(1, 101):
factor1 = math.sin(n * math.pi / 2.0) / (n * n)
factor2 = math.sinh(n * math.pi * y)
factor3 = math.sinh(n * math.pi * (1 - y))
factor4 = math.sin(n * math.pi * x)
term = (factor2 + factor3) / math.sinh(n * math.pi)
sum += factor1 * term * factor4
return 8.0 * sum / (math.pi * math.pi)
avgPE = 0
deltaX = 1.0 / xi
deltaY = 1.0 / yi
maxIts = 50
start_time = time.perf_counter()
for i in range(0, xi + 1):
for j in range(0, yi + 1):
v[i][j] = f(i * deltaX, j * deltaY)
minPE = +1000000000
maxPE = -1000000000
params = {}
Compute(deltaX, deltaY, xi, yi, maxIts, params)
print("Approximate\tAnalytic\tPercent Error")
for i in range(0, xi + 1):
for j in range(0, yi + 1):
if (math.fabs(u[i][j]) > 1.0e-12 and
math.fabs(v[i][j]) > 1.0e-12):
pe = 100.0 * math.fabs((v[i][j] - u[i][j]) / v[i][j])
else:
pe = 0.0
avgPE += pe
if (pe < minPE):
minPE = pe
if (pe > maxPE):
maxPE = pe
if math.fabs(pe) != 0.0:
print("{:10.8f}".format(u[i][j]), "\t", "{:10.8f}".format(v[i][j]), "\t", "{:10.8f}".format(pe))
avgPE /= (xi * yi)
end_time = time.perf_counter()
# Calculate elapsed time in milliseconds
elapsed_ms = (end_time - start_time) * 1000
print("Iterations = ", params['iterations'])
print("Frobenius Norm = ", params['norm'])
print("Minimum Percent Error = ", "{:10.8f}".format(minPE))
print("Average Percent Error = ", "{:10.8f}".format(avgPE))
print("Maximum Percent Error = ", "{:10.8f}".format(maxPE))
print("Elapsed Milliseconds = ", "{:10.8f}".format(elapsed_ms))
Approximate Analytic Percent Error
0.20000000 0.19999972 0.00013831
0.16633455 0.16663592 0.18085704
0.13739427 0.13768928 0.21425509
0.11591159 0.11605132 0.12040154
0.10292732 0.10291871 0.00836375
0.09864668 0.09854114 0.10710198
0.10305612 0.10291871 0.13351975
0.11613103 0.11605132 0.06868445
0.13763395 0.13768928 0.04018144
0.16650309 0.16663592 0.07971563
0.20000000 0.19999972 0.00013831
0.40000000 0.39999927 0.00018169
0.32813882 0.32854798 0.12453472
0.26763492 0.26776105 0.04710489
0.22369849 0.22344522 0.11334607
0.19755216 0.19705465 0.25246951
0.18899096 0.18834090 0.34515051
0.19778524 0.19705465 0.37075402
0.22409557 0.22344522 0.29105301
0.26806863 0.26776105 0.11486989
0.32844379 0.32854798 0.03171204
0.40000000 0.39999927 0.00018169
0.60000000 0.59999811 0.00031532
0.47888974 0.47875999 0.02710181
0.38176991 0.38059768 0.30799710
0.31425594 0.31267225 0.50650231
0.27518847 0.27354681 0.60013876
0.26255192 0.26082096 0.66365837
0.27549367 0.27354681 0.71170992
0.31477587 0.31267225 0.67278721
0.38233779 0.38059768 0.45720485
0.47928905 0.47875999 0.11050679
0.60000000 0.59999811 0.00031532
0.80000000 0.79999202 0.00099701
0.60602379 0.60222488 0.63081180
0.46685956 0.46199176 1.05365615
0.37704317 0.37308607 1.06063957
0.32711029 0.32392135 0.98448025
0.31121899 0.30818168 0.98555933
0.32745161 0.32392135 1.08985035
0.37762462 0.37308607 1.21648870
0.46749464 0.46199176 1.19112076
0.60647034 0.60222488 0.70496220
0.80000000 0.79999202 0.00099701
1.00000000 0.99594729 0.40692036
0.67874673 0.65811281 3.13531720
0.50319768 0.49282441 2.10486107
0.40066316 0.39470092 1.51057096
0.34574699 0.34157202 1.22228048
0.32848272 0.32468552 1.16950228
0.34608839 0.34157202 1.32223096
0.40124475 0.39470092 1.65792212
0.50383291 0.49282441 2.23375660
0.67919339 0.65811281 3.20318626
1.00000000 0.99594729 0.40692036
0.80000000 0.79999202 0.00099701
0.60615279 0.60222488 0.65223206
0.46709299 0.46199176 1.10418350
0.37734885 0.37308607 1.14257161
0.32745218 0.32392135 1.09002599
0.31156100 0.30818168 1.09653477
0.32776105 0.32392135 1.18538026
0.37787503 0.37308607 1.28360546
0.46766769 0.46199176 1.22857864
0.60655688 0.60222488 0.71933113
0.80000000 0.79999202 0.00099701
0.60000000 0.59999811 0.00031532
0.47910949 0.47875999 0.07300025
0.38216756 0.38059768 0.41247566
0.31477665 0.31267225 0.67303766
0.27577086 0.27354681 0.81304189
0.26313452 0.26082096 0.88702848
0.27602079 0.27354681 0.90440981
0.31520243 0.31267225 0.80920943
0.38263258 0.38059768 0.53465917
0.47943646 0.47875999 0.14129609
0.60000000 0.59999811 0.00031532
0.40000000 0.39999927 0.00018169
0.32837881 0.32854798 0.05149022
0.26806920 0.26776105 0.11508237
0.22426717 0.22344522 0.36785094
0.19818820 0.19705465 0.57524397
0.18962723 0.18834090 0.68297865
0.19836093 0.19705465 0.66290001
0.22456142 0.22344522 0.49953914
0.26839057 0.26776105 0.23510691
0.32860478 0.32854798 0.01728752
0.40000000 0.39999927 0.00018169
0.20000000 0.19999972 0.00013831
0.16650327 0.16663592 0.07960469
0.13769959 0.13768928 0.00748883
0.11631140 0.11605132 0.22411146
0.10337449 0.10291871 0.44285504
0.09909401 0.09854114 0.56105764
0.10346087 0.10291871 0.52678322
0.11645855 0.11605132 0.35090580
0.13786030 0.13768928 0.12420922
0.16661627 0.16663592 0.01179315
0.20000000 0.19999972 0.00013831
Iterations = 50
Frobenius Norm = 4.028216200275417
Minimum Percent Error = 0.00000000
Average Percent Error = 0.54286140
Maximum Percent Error = 3.20318626
Elapsed Milliseconds = 36.23520000
Press any key to continue . . .
Blog Entry © Wednesday, May 27, 2026, by James Pate Williams, Jr. and Microsoft’s Copilot Grade School Arithmetic
#pragma once
#include <stdint.h>
/* Algorithm due to Microsft's Coilot
function udiv_restoring(N, D, n) :
R = 0
Q = 0
negD = (~D + 1)
for i from n - 1 down to 0
{
R = (R << 1) | ((N >> i) & 1)
T = R + negD
if MSB(T) == 0:
R = T
Q = Q | (1 << i)
return (Q, R)
*/
class Arithmetic
{
public:
static bool udiv_restoring(
uint32_t numer,
uint32_t denom,
uint32_t& quo,
uint32_t& rem,
int n);
static bool umul_shift_add(
uint32_t a,
uint32_t b,
uint64_t& product,
int n);
};
#include <cstdint>
#include "Arithmetic.h"
static inline uint32_t mask_n(int bits) {
return (bits >= 32) ? 0xFFFFFFFFu : ((1u << bits) - 1u);
}
static inline uint32_t msb(uint32_t x, int bits) {
// returns top bit of a 'bits'-wide value
return (x >> (bits - 1)) & 1u;
}
bool Arithmetic::udiv_restoring(
uint32_t numer,
uint32_t denom,
uint32_t& quo,
uint32_t& rem,
int n)
{
if (denom == 0 || n <= 0 || n > 32) return false;
if (numer == 0)
{
quo = rem = 0;
return true;
}
quo = 0;
rem = 0;
if (n == 32) {
uint64_t R = 0;
uint64_t D = (uint64_t)denom;
uint64_t maskW = (1ull << 33) - 1ull; // 33-bit mask
uint64_t negD = ((~D) + 1ull) & maskW; // 33-bit two's complement
for (int i = n - 1; i >= 0; --i) {
R = ((R << 1) | ((numer >> i) & 1u)) & maskW;
uint64_t T = (R + negD) & maskW; // R - D
// Sign bit is bit 32 (the 33rd bit)
if (((T >> 32) & 1ull) == 0ull) {
R = T;
quo |= (1u << i);
}
}
rem = (uint32_t)(R & 0xFFFFFFFFu);
return true;
}
// n < 32 case: we can keep everything in uint32_t using (n+1) bits
uint32_t maskN = mask_n(n);
uint32_t maskW = mask_n(n + 1);
uint32_t N = numer & maskN;
uint32_t D = denom & maskN;
// Two's complement of D in (n+1) bits
uint32_t Dw = D; // placed in low bits of (n+1)-wide register
uint32_t negD = ((~Dw) + 1u) & maskW;
uint32_t R = 0;
for (int i = n - 1; i >= 0; --i) {
R = ((R << 1) | ((N >> i) & 1u)) & maskW;
uint32_t T = (R + negD) & maskW; // trial subtract: R - D (in w bits)
if (msb(T, n + 1) == 0) { // non-negative in (n+1) bits
R = T;
quo |= (1u << i);
}
}
rem = R & maskN; // remainder fits in n bits
return true;
}
bool Arithmetic::umul_shift_add(
uint32_t a,
uint32_t b,
uint64_t& product,
int n)
{
if (n <= 0 || n > 32) return false;
uint64_t A = a; // promote to avoid overflow
uint32_t B = b;
product = 0;
for (int i = 0; i < n; ++i) {
if (B & 1u) {
product += A;
}
A <<= 1;
B >>= 1;
}
return true;
}
#include <chrono>
#include <cstdint>
#include <iostream>
#include <limits>
#include <random>
#include <string>
#include "Arithmetic.h"
namespace {
constexpr int TESTS_PER_N = 200000;
uint32_t make_mask(int n) {
return (n == 32) ? 0xFFFFFFFFu : ((1u << n) - 1u);
}
void clear_bad_input() {
std::cin.clear();
std::cin.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
}
template <typename TrialFn>
double run_suite(const char* label, TrialFn trial, bool verbose) {
std::mt19937 rng(12345); // deterministic
auto t0 = std::chrono::high_resolution_clock::now();
for (int n = 1; n <= 32; ++n) {
const uint32_t mask = make_mask(n);
for (int i = 0; i < TESTS_PER_N; ++i) {
if (!trial(rng, mask, n)) {
std::cout << label << ": FAILED (n=" << n << ", i=" << i << ")\n";
return -1.0;
}
}
if (verbose) {
std::cout << "n=" << n << " passed\n";
}
}
auto t1 = std::chrono::high_resolution_clock::now();
double secs = std::chrono::duration<double>(t1 - t0).count();
std::cout << label << " runtime = " << secs << " sec\n";
return secs;
}
bool trial_division(std::mt19937& rng, uint32_t mask, int n) {
const uint32_t numer = rng() & mask;
const uint32_t denom = (rng() & mask) | 1u; // ensure non-zero
uint32_t q = 0, r = 0;
if (!Arithmetic::udiv_restoring(numer, denom, q, r, n)) {
std::cout << "Failure numer=" << numer << " denom=" << denom << "\n";
return false;
}
const uint32_t q2 = numer / denom;
const uint32_t r2 = numer % denom;
if (q != q2 || r != r2) {
std::cout << "Mismatch n=" << n
<< " numer=" << numer
<< " denom=" << denom
<< " got q=" << q << " r=" << r
<< " expected q=" << q2 << " r=" << r2 << "\n";
return false;
}
return true;
}
bool trial_multiplication(std::mt19937& rng, uint32_t mask, int n) {
const uint32_t a = rng() & mask;
const uint32_t b = rng() & mask;
uint64_t prod = 0;
if (!Arithmetic::umul_shift_add(a, b, prod, n)) {
std::cout << "Failure a=" << a << " b=" << b << "\n";
return false;
}
const uint64_t expected = static_cast<uint64_t>(a) * static_cast<uint64_t>(b);
if (prod != expected) {
std::cout << "Mismatch n=" << n
<< " a=" << a
<< " b=" << b
<< " got=" << prod
<< " expected=" << expected << "\n";
return false;
}
return true;
}
} // namespace
int main() {
bool verbose = true;
while (true) {
std::cout << "\nArithmetic Lab\n";
std::cout << "1. Test Division (restoring)\n";
std::cout << "2. Test Multiplication (shift-add)\n";
std::cout << "3. Run ALL tests\n";
std::cout << "4. Toggle verbose (currently " << (verbose ? "ON" : "OFF") << ")\n";
std::cout << "5. Exit\n";
std::cout << "Choice: ";
int choice = 0;
if (!(std::cin >> choice)) {
clear_bad_input();
std::cout << "Invalid input. Please enter a number.\n";
continue;
}
if (choice == 1) {
run_suite("Division test", trial_division, verbose);
}
else if (choice == 2) {
run_suite("Multiplication test", trial_multiplication, verbose);
}
else if (choice == 3) {
const double d = run_suite("Division test", trial_division, verbose);
if (d >= 0.0) run_suite("Multiplication test", trial_multiplication, verbose);
}
else if (choice == 4) {
verbose = !verbose;
}
else if (choice == 5) {
return 0;
}
else {
std::cout << "Invalid choice.\n";
}
}
}
Blog Entry (c) Tuesday, May 26, 2026, by James Pate Williams, Jr. and Microsoft’s Copilot Hydrogen-Like Atomic Radial Wave Functions
Included a downloadable PDF and Microsoft Excel Workbook.
Numerical Explorations 