1- // DCT-IV / half-shift (BDCT ) benchmark and validation.
1+ // BODFT (odd-frequency DFT ) benchmark and validation.
22//
3- // Validates the paired radix-4 half-shift kernel and the DCT-IV built on it,
4- // in both double and single precision, without any long double (there is no
5- // extended precision on the NEON target). Correctness is checked three ways:
3+ // Validates the iterative paired radix-4 BODFT kernel in both double and single
4+ // precision, without any long double (there is no extended precision on the NEON
5+ // target). Correctness is checked three ways:
66//
7- // 1. analytic exact: half-shift of a unit impulse has |H[k]| == 1 with a
8- // closed-form phase, and DCT-IV of a single basis vector is a clean spike ;
7+ // 1. analytic exact: the BODFT of a unit impulse has |H[k]| == 1 with a
8+ // closed-form phase;
99// 2. independent cross-check against a direct same-precision DFT;
1010// 3. exact inverse roundtrip.
1111//
1212// Build (from repo root):
13- // c++ -O3 -std=c++17 -Iinclude examples/dctiv_benchmark .cpp -lm -o build/examples/dctiv_benchmark
13+ // c++ -O3 -std=c++17 -Iinclude examples/bodft_benchmark .cpp -lm -o build/examples/bodft_benchmark
1414
15- #include " ../src/detail/bdct_kernel .hpp"
15+ #include " ../src/detail/bodft_kernel .hpp"
1616
1717#include < chrono>
1818#include < cmath>
@@ -33,48 +33,29 @@ std::vector<RT> make_input(int n, std::mt19937_64& rng) {
3333 return x;
3434}
3535
36- // --- correctness for one precision -----------------------------------------
37-
38- template <class HS , class DCT , class CT , class RT >
36+ template <class Plan , class CT , class RT >
3937void run_correctness (int N, std::mt19937_64& rng, const char * tag) {
40- HS hs (N);
41- DCT dct (N);
38+ Plan p (N);
4239 const int bins = N / 2 ;
4340
44- // (1a) half-shift of a unit impulse at n0: H[k] = exp(-2*pi*i*(k+1/2)*n0/N).
41+ // (1) BODFT of a unit impulse at n0: H[k] = exp(-2*pi*i*(k+1/2)*n0/N).
4542 std::vector<RT > imp (static_cast <std::size_t >(N), static_cast <RT >(0 ));
4643 const int n0 = 1 ;
4744 imp[static_cast <std::size_t >(n0)] = static_cast <RT >(1 );
4845 std::vector<CT > Himp (static_cast <std::size_t >(bins));
49- hs .forward (imp.data (), Himp.data ());
46+ p .forward (imp.data (), Himp.data ());
5047 double imp_err = 0.0 ;
5148 for (int k = 0 ; k < bins; ++k) {
52- const double ang = -2.0 * PI * (static_cast <double >(k) + 0.5 ) *
53- static_cast <double >(n0) / static_cast <double >(N);
49+ const double ang = -2.0 * PI * (k + 0.5 ) * n0 / N;
5450 const double dr = static_cast <double >(Himp[static_cast <std::size_t >(k)].re ) - std::cos (ang);
5551 const double di = static_cast <double >(Himp[static_cast <std::size_t >(k)].im ) - std::sin (ang);
5652 imp_err = std::max (imp_err, std::sqrt (dr * dr + di * di));
5753 }
5854
59- // (1b) DCT-IV of basis vector m: x[n]=cos(pi*(m+1/2)*(n+1/2)/N) -> (N/2)*delta_km.
60- const int m = (N >= 8 ) ? 3 : 1 ;
61- std::vector<RT > basis (static_cast <std::size_t >(N));
62- for (int n = 0 ; n < N; ++n) {
63- basis[static_cast <std::size_t >(n)] = static_cast <RT >(
64- std::cos (PI * (m + 0.5 ) * (n + 0.5 ) / N));
65- }
66- std::vector<RT > c (static_cast <std::size_t >(N));
67- dct.forward (basis.data (), c.data ());
68- double dct_err = 0.0 ;
69- for (int k = 0 ; k < N; ++k) {
70- const double expect = (k == m) ? 0.5 * N : 0.0 ;
71- dct_err = std::max (dct_err, std::fabs (static_cast <double >(c[static_cast <std::size_t >(k)]) - expect));
72- }
73-
7455 // (2) random signal vs direct same-precision DFT, and (3) roundtrip.
7556 std::vector<RT > x = make_input<RT >(N, rng);
7657 std::vector<CT > H (static_cast <std::size_t >(bins));
77- hs .forward (x.data (), H.data ());
58+ p .forward (x.data (), H.data ());
7859 double xcheck = 0.0 , ref_mag = 0.0 ;
7960 for (int k = 0 ; k < bins; ++k) {
8061 double sr = 0.0 , si = 0.0 ;
@@ -90,29 +71,24 @@ void run_correctness(int N, std::mt19937_64& rng, const char* tag) {
9071 }
9172
9273 std::vector<RT > xr (static_cast <std::size_t >(N));
93- hs .inverse (H.data (), xr.data ());
74+ p .inverse (H.data (), xr.data ());
9475 double round_err = 0.0 ;
9576 for (int n = 0 ; n < N; ++n) {
9677 round_err = std::max (round_err,
9778 std::fabs (static_cast <double >(xr[static_cast <std::size_t >(n)] -
9879 x[static_cast <std::size_t >(n)])));
9980 }
10081
101- std::printf (" %-6s N %7d hs_impulse %.3e dctiv_basis %.3e "
102- " hs_vs_dft(rel) %.3e roundtrip %.3e\n "
103- tag, N, imp_err, dct_err,
104- ref_mag > 0 ? xcheck / ref_mag : 0.0 , round_err);
82+ std::printf (" %-6s N %7d impulse %.3e vs_dft(rel) %.3e roundtrip %.3e\n "
83+ tag, N, imp_err, ref_mag > 0 ? xcheck / ref_mag : 0.0 , round_err);
10584}
10685
107- // --- timing ----------------------------------------------------------------
108-
109- template <class HS , class DCT , class CT , class RT >
86+ template <class Plan , class CT , class RT >
11087void run_timing (int N, std::mt19937_64& rng, const char * tag) {
11188 std::vector<RT > x = make_input<RT >(N, rng);
112- HS hs (N);
113- DCT dct (N);
89+ Plan p (N);
11490 std::vector<CT > H (static_cast <std::size_t >(N / 2 ));
115- std::vector<RT > c (static_cast <std::size_t >(N));
91+ std::vector<RT > xr (static_cast <std::size_t >(N));
11692
11793 long iters = 1 ;
11894 {
@@ -122,19 +98,19 @@ void run_timing(int N, std::mt19937_64& rng, const char* tag) {
12298 }
12399
124100 using clock = std::chrono::steady_clock;
125- hs .forward (x.data (), H.data ());
126- dct. forward (x .data (), c .data ());
101+ p .forward (x.data (), H.data ());
102+ p. inverse (H .data (), xr .data ());
127103
128104 auto t0 = clock::now ();
129- for (long i = 0 ; i < iters; ++i) hs .forward (x.data (), H.data ());
105+ for (long i = 0 ; i < iters; ++i) p .forward (x.data (), H.data ());
130106 auto t1 = clock::now ();
131- for (long i = 0 ; i < iters; ++i) dct. forward (x .data (), c .data ());
107+ for (long i = 0 ; i < iters; ++i) p. inverse (H .data (), xr .data ());
132108 auto t2 = clock::now ();
133109
134- const double hs_ns = std::chrono::duration<double , std::nano>(t1 - t0).count () / iters;
135- const double dct_ns = std::chrono::duration<double , std::nano>(t2 - t1).count () / iters;
136- std::printf (" %-6s N %7d halfshift %10.1f ns dctiv %10.1f ns (%ld iters)\n "
137- tag, N, hs_ns, dct_ns , iters);
110+ const double fwd_ns = std::chrono::duration<double , std::nano>(t1 - t0).count () / iters;
111+ const double inv_ns = std::chrono::duration<double , std::nano>(t2 - t1).count () / iters;
112+ std::printf (" %-6s N %7d forward %10.1f ns inverse %10.1f ns (%ld iters)\n "
113+ tag, N, fwd_ns, inv_ns , iters);
138114}
139115
140116} // namespace
@@ -145,20 +121,20 @@ int main(int argc, char** argv) {
145121 std::printf (" backend: %s\n " bruun::simd_backend_name ());
146122
147123 std::printf (" \n == correctness (double) ==\n "
148- for (int N = 4 ; N <= 8192 ; N <<= 1 )
149- run_correctness<bdct::half_shift_plan, bdct::dctiv_plan, bdct ::complex_t , double >(N, rng, " f64"
124+ for (int N = 2 ; N <= 8192 ; N <<= 1 )
125+ run_correctness<bodft::plan, bodft ::complex_t , double >(N, rng, " f64"
150126
151127 std::printf (" \n == correctness (float) ==\n "
152- for (int N = 4 ; N <= 8192 ; N <<= 1 )
153- run_correctness<bdct::half_shift_plan_f32, bdct::dctiv_plan_f32, bdct ::complex_f32_t , float >(N, rng, " f32"
128+ for (int N = 2 ; N <= 8192 ; N <<= 1 )
129+ run_correctness<bodft::plan_f32, bodft ::complex_f32_t , float >(N, rng, " f32"
154130
155131 std::printf (" \n == timing (double) ==\n "
156132 for (int N = 64 ; N <= (1 << 20 ); N <<= 1 )
157- run_timing<bdct::half_shift_plan, bdct::dctiv_plan, bdct ::complex_t , double >(N, rng, " f64"
133+ run_timing<bodft::plan, bodft ::complex_t , double >(N, rng, " f64"
158134
159135 std::printf (" \n == timing (float) ==\n "
160136 for (int N = 64 ; N <= (1 << 20 ); N <<= 1 )
161- run_timing<bdct::half_shift_plan_f32, bdct::dctiv_plan_f32, bdct ::complex_f32_t , float >(N, rng, " f32"
137+ run_timing<bodft::plan_f32, bodft ::complex_f32_t , float >(N, rng, " f32"
162138
163139 (void )argc;
164140 (void )argv;
0 commit comments