From ed93c83a7ac559a6598b76f9e82b3ca6bc9e6996 Mon Sep 17 00:00:00 2001 From: Mitchell Thompkins Date: Wed, 24 Jun 2026 18:20:31 +0000 Subject: [PATCH 01/11] Evaluate Faddeev-LeVerrier char_poly via consteig, no accuracy gain --- .../profiling/results/accuracy_gcc_15.2.0.csv | 80 +++++++------- .../profiling/results/accuracy_gcc_15.2.0.png | 4 +- .../results/compile_times_gcc_15.2.0.csv | 102 +++++++++--------- .../results/compile_times_gcc_15.2.0.png | 4 +- docs/profiling/results/runtime_gcc_15.2.0.csv | 78 +++++++------- docs/profiling/results/runtime_gcc_15.2.0.png | 4 +- include/constfilt/discretize.hpp | 30 +----- .../vendor/consteig/matrix/operations.hpp | 30 ++++++ 8 files changed, 170 insertions(+), 162 deletions(-) diff --git a/docs/profiling/results/accuracy_gcc_15.2.0.csv b/docs/profiling/results/accuracy_gcc_15.2.0.csv index dcfad8b..0e826b3 100644 --- a/docs/profiling/results/accuracy_gcc_15.2.0.csv +++ b/docs/profiling/results/accuracy_gcc_15.2.0.csv @@ -4,16 +4,16 @@ library,filter_type,order,method,max_b_err,max_a_err,max_step_err constfilt,butterworth,1,zoh,0.000000e+00,0.000000e+00,0.000000e+00 constfilt,butterworth,2,zoh,3.608225e-16,2.220446e-16,1.110223e-15 -constfilt,butterworth,3,zoh,9.020562e-16,1.776357e-15,7.771561e-15 -constfilt,butterworth,4,zoh,8.971990e-15,8.881784e-16,1.587619e-14 -constfilt,butterworth,5,zoh,6.945833e-15,3.552714e-15,1.809664e-14 -constfilt,butterworth,6,zoh,3.284352e-14,1.598721e-14,6.572520e-14 -constfilt,butterworth,7,zoh,1.100200e-11,4.796163e-14,1.254094e-12 -constfilt,butterworth,8,zoh,2.955275e-08,5.062617e-13,2.134253e-09 -constfilt,butterworth,9,zoh,7.577398e-05,8.526513e-13,3.692378e-06 -constfilt,butterworth,10,zoh,2.065719e-01,2.030731e-11,6.461003e-03 -constfilt,butterworth,11,zoh,2.543305e+02,3.566485e-08,4.685698e+00 -constfilt,butterworth,12,zoh,1.799631e+03,3.915801e-10,2.955113e+02 +constfilt,butterworth,3,zoh,9.020562e-16,1.332268e-15,1.088019e-14 +constfilt,butterworth,4,zoh,8.971990e-15,1.110223e-15,1.321165e-14 +constfilt,butterworth,5,zoh,6.944098e-15,7.549517e-15,1.543210e-14 +constfilt,butterworth,6,zoh,3.286173e-14,9.769963e-15,5.095924e-14 +constfilt,butterworth,7,zoh,1.100198e-11,2.131628e-14,1.254094e-12 +constfilt,butterworth,8,zoh,2.955275e-08,4.849454e-13,2.134253e-09 +constfilt,butterworth,9,zoh,7.577398e-05,8.100187e-13,3.692378e-06 +constfilt,butterworth,10,zoh,2.065719e-01,2.042100e-11,6.461003e-03 +constfilt,butterworth,11,zoh,2.543305e+02,3.731445e-07,4.685698e+00 +constfilt,butterworth,12,zoh,1.799631e+03,3.914806e-10,2.955113e+02 constfilt,butterworth,1,matchedz,5.551115e-17,0.000000e+00,2.220446e-16 constfilt,butterworth,2,matchedz,2.775558e-17,2.220446e-16,1.110223e-15 constfilt,butterworth,3,matchedz,1.387779e-17,1.776357e-15,1.998401e-15 @@ -27,28 +27,28 @@ constfilt,butterworth,10,matchedz,1.523304e-16,5.579892e-11,1.363021e-11 constfilt,butterworth,11,matchedz,2.168404e-19,1.847411e-13,4.028200e-11 constfilt,butterworth,12,matchedz,9.662952e-17,3.253149e-10,1.482581e-10 constfilt,butterworth,1,prewarp,8.326673e-17,2.220446e-16,2.220446e-16 -constfilt,butterworth,2,prewarp,1.249001e-16,0.000000e+00,2.331468e-15 -constfilt,butterworth,3,prewarp,1.283695e-16,2.442491e-15,1.776357e-15 -constfilt,butterworth,4,prewarp,1.630640e-16,3.108624e-15,1.432188e-14 -constfilt,butterworth,5,prewarp,1.647987e-16,8.881784e-16,2.930989e-14 -constfilt,butterworth,6,prewarp,1.231654e-16,1.154632e-14,1.159073e-13 -constfilt,butterworth,7,prewarp,2.493665e-16,2.664535e-14,1.717515e-13 -constfilt,butterworth,8,prewarp,5.490400e-16,5.684342e-14,5.651035e-13 -constfilt,butterworth,9,prewarp,4.832289e-16,2.842171e-14,4.035883e-12 -constfilt,butterworth,10,prewarp,3.041296e-15,2.486900e-14,6.625700e-12 -constfilt,butterworth,11,prewarp,1.426024e-15,3.197442e-13,1.841038e-11 -constfilt,butterworth,12,prewarp,1.350920e-14,1.563194e-13,5.809664e-11 +constfilt,butterworth,2,prewarp,1.249001e-16,5.551115e-17,1.998401e-15 +constfilt,butterworth,3,prewarp,5.204170e-17,4.440892e-16,3.330669e-15 +constfilt,butterworth,4,prewarp,1.734723e-16,2.220446e-15,1.598721e-14 +constfilt,butterworth,5,prewarp,1.682682e-16,2.220446e-15,2.131628e-14 +constfilt,butterworth,6,prewarp,1.214306e-16,1.065814e-14,3.508305e-14 +constfilt,butterworth,7,prewarp,2.474149e-16,7.993606e-15,2.481348e-13 +constfilt,butterworth,8,prewarp,5.555452e-16,3.019807e-14,7.693846e-13 +constfilt,butterworth,9,prewarp,4.776995e-16,2.131628e-14,4.931056e-12 +constfilt,butterworth,10,prewarp,3.020696e-15,3.552714e-14,4.900191e-12 +constfilt,butterworth,11,prewarp,1.400681e-15,6.394885e-14,3.339040e-11 +constfilt,butterworth,12,prewarp,1.350855e-14,1.705303e-13,3.085243e-11 constfilt,elliptic,2,zoh,2.395861e-13,2.842171e-14,1.409099e-13 -constfilt,elliptic,3,zoh,5.400819e-13,9.858780e-14,3.640976e-13 -constfilt,elliptic,4,zoh,5.393082e-13,1.074696e-13,8.994472e-13 -constfilt,elliptic,5,zoh,1.066647e-12,4.938272e-13,1.474598e-12 -constfilt,elliptic,6,zoh,1.760619e-12,1.481482e-12,1.819322e-12 -constfilt,elliptic,7,zoh,1.055545e-11,1.796252e-11,9.916956e-12 -constfilt,elliptic,8,zoh,2.154665e-09,2.641301e-09,7.933327e-10 -constfilt,elliptic,9,zoh,3.582365e-07,1.304602e-07,9.765907e-07 -constfilt,elliptic,10,zoh,3.461312e-03,3.213908e-06,6.243202e-03 -constfilt,elliptic,11,zoh,5.384745e-01,4.802802e-05,7.544774e-02 -constfilt,elliptic,12,zoh,5.474468e+02,5.003185e-04,2.311049e+01 +constfilt,elliptic,3,zoh,5.399570e-13,9.481305e-14,3.640976e-13 +constfilt,elliptic,4,zoh,5.392874e-13,1.083578e-13,8.994472e-13 +constfilt,elliptic,5,zoh,1.066758e-12,4.858336e-13,1.475042e-12 +constfilt,elliptic,6,zoh,1.760855e-12,1.453060e-12,1.816436e-12 +constfilt,elliptic,7,zoh,1.055522e-11,1.797673e-11,9.918288e-12 +constfilt,elliptic,8,zoh,2.154666e-09,2.641208e-09,7.933643e-10 +constfilt,elliptic,9,zoh,3.582365e-07,1.304601e-07,9.761215e-07 +constfilt,elliptic,10,zoh,3.461312e-03,3.213907e-06,6.241612e-03 +constfilt,elliptic,11,zoh,5.384745e-01,4.802802e-05,7.545471e-02 +constfilt,elliptic,12,zoh,5.474468e+02,5.003185e-04,2.311052e+01 constfilt,elliptic,2,matchedz,3.854481e-11,2.786660e-14,1.928402e-11 constfilt,elliptic,3,matchedz,5.143178e-13,9.592327e-14,3.533007e-13 constfilt,elliptic,4,matchedz,1.324132e-13,1.101341e-13,8.861800e-13 @@ -61,15 +61,15 @@ constfilt,elliptic,10,matchedz,7.442555e-07,3.213346e-06,2.409790e-07 constfilt,elliptic,11,matchedz,2.360767e-05,4.802825e-05,1.950004e-06 constfilt,elliptic,12,matchedz,1.156955e-04,5.002874e-04,1.090689e-05 constfilt,elliptic,2,prewarp,1.637024e-13,3.286260e-14,9.213463e-14 -constfilt,elliptic,3,prewarp,1.596015e-13,8.082424e-14,3.510803e-13 -constfilt,elliptic,4,prewarp,3.788428e-13,1.079137e-13,8.840706e-13 -constfilt,elliptic,5,prewarp,2.870135e-13,4.884981e-13,1.437517e-12 -constfilt,elliptic,6,prewarp,1.340761e-12,1.419309e-12,1.797784e-12 -constfilt,elliptic,7,prewarp,2.618128e-12,1.753975e-11,9.648726e-12 -constfilt,elliptic,8,prewarp,1.244403e-09,2.586319e-09,7.183503e-10 -constfilt,elliptic,9,prewarp,1.482448e-08,1.278831e-07,1.805560e-08 -constfilt,elliptic,10,prewarp,1.376367e-06,3.130827e-06,2.329530e-07 -constfilt,elliptic,11,prewarp,4.570721e-06,4.696437e-05,1.897978e-06 +constfilt,elliptic,3,prewarp,1.595737e-13,8.104628e-14,3.510803e-13 +constfilt,elliptic,4,prewarp,3.788705e-13,1.092459e-13,8.841816e-13 +constfilt,elliptic,5,prewarp,2.869857e-13,4.902745e-13,1.438516e-12 +constfilt,elliptic,6,prewarp,1.340678e-12,1.428191e-12,1.799783e-12 +constfilt,elliptic,7,prewarp,2.617934e-12,1.758949e-11,9.648282e-12 +constfilt,elliptic,8,prewarp,1.244404e-09,2.586262e-09,7.183454e-10 +constfilt,elliptic,9,prewarp,1.482448e-08,1.278829e-07,1.805562e-08 +constfilt,elliptic,10,prewarp,1.376367e-06,3.130828e-06,2.329530e-07 +constfilt,elliptic,11,prewarp,4.570721e-06,4.696437e-05,1.897979e-06 constfilt,elliptic,12,prewarp,2.111470e-04,4.866247e-04,1.055330e-05 iir1,butterworth,1,prewarp,n/a,n/a,4.440892e-16 iir1,butterworth,2,prewarp,n/a,n/a,8.881784e-16 diff --git a/docs/profiling/results/accuracy_gcc_15.2.0.png b/docs/profiling/results/accuracy_gcc_15.2.0.png index 996c68d..f63dbb9 100644 --- a/docs/profiling/results/accuracy_gcc_15.2.0.png +++ b/docs/profiling/results/accuracy_gcc_15.2.0.png @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:3063f67ee8c7ca693c9bc690f2bf9e9e3a35e86b8025f1938552ae6ef2271908 -size 169290 +oid sha256:89d68110f966dd40c7c502e1c38033ce69e54d7985522432ef624ad7d8e9af2a +size 169078 diff --git a/docs/profiling/results/compile_times_gcc_15.2.0.csv b/docs/profiling/results/compile_times_gcc_15.2.0.csv index d7f20e2..183c566 100644 --- a/docs/profiling/results/compile_times_gcc_15.2.0.csv +++ b/docs/profiling/results/compile_times_gcc_15.2.0.csv @@ -2,54 +2,54 @@ # compiler: g++ (Alpine 15.2.0) 15.2.0 # os: Alpine Linux filter_type,order,method,compile_time_sec,max_rss_kb,exit_code -butterworth,10,matchedz,0.37,73000,0 -butterworth,10,tustin,0.45,83112,0 -butterworth,10,zoh,1.07,158552,0 -butterworth,11,matchedz,0.50,88108,0 -butterworth,11,tustin,0.62,102420,0 -butterworth,11,zoh,1.42,203688,0 -butterworth,12,matchedz,0.63,105460,0 -butterworth,12,tustin,0.81,126272,0 -butterworth,12,zoh,1.93,264624,0 -butterworth,1,matchedz,0.07,36764,0 -butterworth,1,tustin,0.06,36776,0 -butterworth,1,zoh,0.07,37024,0 -butterworth,2,matchedz,0.07,37096,0 -butterworth,2,tustin,0.06,37268,0 -butterworth,2,zoh,0.07,38084,0 -butterworth,4,matchedz,0.09,39528,0 -butterworth,4,tustin,0.09,40116,0 -butterworth,4,zoh,0.14,44948,0 -butterworth,6,matchedz,0.12,44468,0 -butterworth,6,tustin,0.14,45600,0 -butterworth,6,zoh,0.27,60956,0 -butterworth,8,matchedz,0.21,54360,0 -butterworth,8,tustin,0.25,58428,0 -butterworth,8,zoh,0.54,94996,0 -butterworth,9,matchedz,0.28,62652,0 -butterworth,9,tustin,0.34,69284,0 -butterworth,9,zoh,0.77,123156,0 -elliptic,10,matchedz,0.72,112964,0 -elliptic,10,tustin,0.53,89548,0 -elliptic,10,zoh,1.13,163684,0 -elliptic,11,matchedz,0.95,140320,0 -elliptic,11,tustin,0.69,109080,0 -elliptic,11,zoh,1.52,210084,0 -elliptic,12,matchedz,1.30,181428,0 -elliptic,12,tustin,0.92,134184,0 -elliptic,12,zoh,2.03,271788,0 -elliptic,2,matchedz,0.11,40780,0 -elliptic,2,tustin,0.10,40608,0 -elliptic,2,zoh,0.11,41888,0 -elliptic,4,matchedz,0.14,45328,0 -elliptic,4,tustin,0.13,43936,0 -elliptic,4,zoh,0.18,49076,0 -elliptic,6,matchedz,0.22,53920,0 -elliptic,6,tustin,0.20,50064,0 -elliptic,6,zoh,0.31,64840,0 -elliptic,8,matchedz,0.40,74308,0 -elliptic,8,tustin,0.31,63908,0 -elliptic,8,zoh,0.61,100160,0 -elliptic,9,matchedz,0.53,89588,0 -elliptic,9,tustin,0.40,74688,0 -elliptic,9,zoh,0.84,128436,0 +butterworth,10,matchedz,0.37,72924,0 +butterworth,10,tustin,0.26,60356,0 +butterworth,10,zoh,0.89,139956,0 +butterworth,11,matchedz,0.49,88088,0 +butterworth,11,tustin,0.34,70304,0 +butterworth,11,zoh,1.21,177720,0 +butterworth,12,matchedz,0.64,105620,0 +butterworth,12,tustin,0.45,83064,0 +butterworth,12,zoh,1.59,225348,0 +butterworth,1,matchedz,0.07,36864,0 +butterworth,1,tustin,0.06,36184,0 +butterworth,1,zoh,0.07,37032,0 +butterworth,2,matchedz,0.07,37020,0 +butterworth,2,tustin,0.06,36620,0 +butterworth,2,zoh,0.08,38152,0 +butterworth,4,matchedz,0.09,39224,0 +butterworth,4,tustin,0.07,37548,0 +butterworth,4,zoh,0.12,44124,0 +butterworth,6,matchedz,0.13,44296,0 +butterworth,6,tustin,0.10,40444,0 +butterworth,6,zoh,0.25,58156,0 +butterworth,8,matchedz,0.21,54336,0 +butterworth,8,tustin,0.15,47184,0 +butterworth,8,zoh,0.48,86952,0 +butterworth,9,matchedz,0.29,62732,0 +butterworth,9,tustin,0.19,52700,0 +butterworth,9,zoh,0.67,109980,0 +elliptic,10,matchedz,0.73,112956,0 +elliptic,10,tustin,0.33,66488,0 +elliptic,10,zoh,0.98,145432,0 +elliptic,11,matchedz,0.93,140208,0 +elliptic,11,tustin,0.41,75984,0 +elliptic,11,zoh,1.29,184024,0 +elliptic,12,matchedz,1.30,181384,0 +elliptic,12,tustin,0.54,89480,0 +elliptic,12,zoh,1.68,232768,0 +elliptic,2,matchedz,0.10,40604,0 +elliptic,2,tustin,0.09,40108,0 +elliptic,2,zoh,0.11,42152,0 +elliptic,4,matchedz,0.15,45480,0 +elliptic,4,tustin,0.12,41884,0 +elliptic,4,zoh,0.17,48260,0 +elliptic,6,matchedz,0.22,53660,0 +elliptic,6,tustin,0.15,45224,0 +elliptic,6,zoh,0.29,62044,0 +elliptic,8,matchedz,0.41,74268,0 +elliptic,8,tustin,0.21,52560,0 +elliptic,8,zoh,0.53,91808,0 +elliptic,9,matchedz,0.52,89608,0 +elliptic,9,tustin,0.26,58540,0 +elliptic,9,zoh,0.73,115224,0 diff --git a/docs/profiling/results/compile_times_gcc_15.2.0.png b/docs/profiling/results/compile_times_gcc_15.2.0.png index 0a1b52f..d2dce27 100644 --- a/docs/profiling/results/compile_times_gcc_15.2.0.png +++ b/docs/profiling/results/compile_times_gcc_15.2.0.png @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:4a4f746ddc0d9c0f7f9380a8c7fe02c44288304df30b974042213c7456a1af73 -size 150965 +oid sha256:b2d38127103dd59828b20e490ebe7f6c1b1f2ca00e8e3be08d90c51f02bf28fc +size 147973 diff --git a/docs/profiling/results/runtime_gcc_15.2.0.csv b/docs/profiling/results/runtime_gcc_15.2.0.csv index 5b7cc99..35a0a76 100644 --- a/docs/profiling/results/runtime_gcc_15.2.0.csv +++ b/docs/profiling/results/runtime_gcc_15.2.0.csv @@ -2,42 +2,42 @@ # compiler: g++ (Alpine 15.2.0) 15.2.0 # os: Alpine Linux library,filter_type,order,method,ns_per_sample,msa_per_s,dc_gain -constfilt,Butterworth,1,tustin,6.415,155.9,1.00000000 -constfilt,Butterworth,2,tustin,10.876,91.9,1.00000000 -constfilt,Butterworth,4,tustin,18.954,52.8,1.00000000 -constfilt,Butterworth,6,tustin,26.164,38.2,1.00000000 -constfilt,Butterworth,8,tustin,34.550,28.9,1.00000000 -constfilt,Butterworth,1,zoh,6.420,155.8,1.00000000 -constfilt,Butterworth,2,zoh,10.826,92.4,1.00000000 -constfilt,Butterworth,4,zoh,18.535,54.0,1.00000000 -constfilt,Butterworth,6,zoh,26.066,38.4,1.00000000 -constfilt,Butterworth,8,zoh,34.415,29.1,1.00000000 -constfilt,Butterworth,1,matchedz,6.423,155.7,1.00000000 -constfilt,Butterworth,2,matchedz,10.839,92.3,1.00000000 -constfilt,Butterworth,4,matchedz,18.858,53.0,1.00000000 -constfilt,Butterworth,6,matchedz,26.130,38.3,1.00000000 -constfilt,Butterworth,8,matchedz,34.357,29.1,1.00000000 -constfilt,Elliptic,2,tustin,10.844,92.2,0.94406088 -constfilt,Elliptic,4,tustin,18.373,54.4,0.94406088 -constfilt,Elliptic,6,tustin,26.217,38.1,0.94406088 -constfilt,Elliptic,8,tustin,34.485,29.0,0.94406088 -constfilt,Elliptic,2,zoh,10.828,92.4,0.94406088 -constfilt,Elliptic,4,zoh,18.359,54.5,0.94406088 -constfilt,Elliptic,6,zoh,26.095,38.3,0.94406088 -constfilt,Elliptic,8,zoh,34.249,29.2,0.94406088 -constfilt,Elliptic,2,matchedz,11.091,90.2,0.94406088 -constfilt,Elliptic,4,matchedz,18.904,52.9,0.94406088 -constfilt,Elliptic,6,matchedz,26.024,38.4,0.94406088 -constfilt,Elliptic,8,matchedz,34.507,29.0,0.94406088 -iir1,Butterworth,2,runtime,10.732,93.2,1.00000000 -iir1,Butterworth,4,runtime,18.155,55.1,1.00000000 -iir1,Butterworth,6,runtime,30.092,33.2,1.00000000 -iir1,Butterworth,8,runtime,40.322,24.8,1.00000000 -kfr,Butterworth,2,runtime+simd,607.261,1.6,1.00000000 -kfr,Butterworth,4,runtime+simd,402.769,2.5,1.00000000 -kfr,Butterworth,6,runtime+simd,307.412,3.3,1.00000000 -kfr,Butterworth,8,runtime+simd,307.261,3.3,1.00000000 -kfr,Elliptic,2,runtime+simd,605.280,1.7,0.94406088 -kfr,Elliptic,4,runtime+simd,401.691,2.5,0.94406088 -kfr,Elliptic,6,runtime+simd,306.657,3.3,0.94406088 -kfr,Elliptic,8,runtime+simd,308.208,3.2,0.94406088 +constfilt,Butterworth,1,tustin,6.418,155.8,1.00000000 +constfilt,Butterworth,2,tustin,10.855,92.1,1.00000000 +constfilt,Butterworth,4,tustin,18.741,53.4,1.00000000 +constfilt,Butterworth,6,tustin,26.035,38.4,1.00000000 +constfilt,Butterworth,8,tustin,34.363,29.1,1.00000000 +constfilt,Butterworth,1,zoh,6.486,154.2,1.00000000 +constfilt,Butterworth,2,zoh,10.817,92.4,1.00000000 +constfilt,Butterworth,4,zoh,18.597,53.8,1.00000000 +constfilt,Butterworth,6,zoh,25.921,38.6,1.00000000 +constfilt,Butterworth,8,zoh,34.208,29.2,1.00000000 +constfilt,Butterworth,1,matchedz,6.421,155.7,1.00000000 +constfilt,Butterworth,2,matchedz,10.851,92.2,1.00000000 +constfilt,Butterworth,4,matchedz,18.754,53.3,1.00000000 +constfilt,Butterworth,6,matchedz,26.121,38.3,1.00000000 +constfilt,Butterworth,8,matchedz,34.793,28.7,1.00000000 +constfilt,Elliptic,2,tustin,10.829,92.3,0.94406088 +constfilt,Elliptic,4,tustin,18.308,54.6,0.94406088 +constfilt,Elliptic,6,tustin,26.122,38.3,0.94406088 +constfilt,Elliptic,8,tustin,34.242,29.2,0.94406088 +constfilt,Elliptic,2,zoh,10.817,92.5,0.94406088 +constfilt,Elliptic,4,zoh,18.313,54.6,0.94406088 +constfilt,Elliptic,6,zoh,26.092,38.3,0.94406088 +constfilt,Elliptic,8,zoh,34.251,29.2,0.94406088 +constfilt,Elliptic,2,matchedz,11.076,90.3,0.94406088 +constfilt,Elliptic,4,matchedz,18.853,53.0,0.94406088 +constfilt,Elliptic,6,matchedz,25.973,38.5,0.94406088 +constfilt,Elliptic,8,matchedz,34.539,29.0,0.94406088 +iir1,Butterworth,2,runtime,10.830,92.3,1.00000000 +iir1,Butterworth,4,runtime,18.134,55.1,1.00000000 +iir1,Butterworth,6,runtime,29.181,34.3,1.00000000 +iir1,Butterworth,8,runtime,40.300,24.8,1.00000000 +kfr,Butterworth,2,runtime+simd,588.824,1.7,1.00000000 +kfr,Butterworth,4,runtime+simd,381.513,2.6,1.00000000 +kfr,Butterworth,6,runtime+simd,297.600,3.4,1.00000000 +kfr,Butterworth,8,runtime+simd,296.401,3.4,1.00000000 +kfr,Elliptic,2,runtime+simd,589.859,1.7,0.94406088 +kfr,Elliptic,4,runtime+simd,382.360,2.6,0.94406088 +kfr,Elliptic,6,runtime+simd,296.054,3.4,0.94406088 +kfr,Elliptic,8,runtime+simd,297.250,3.4,0.94406088 diff --git a/docs/profiling/results/runtime_gcc_15.2.0.png b/docs/profiling/results/runtime_gcc_15.2.0.png index 931a63a..3627dcb 100644 --- a/docs/profiling/results/runtime_gcc_15.2.0.png +++ b/docs/profiling/results/runtime_gcc_15.2.0.png @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:e23c6d526144320c1c2baf0928589539fa06e92946d213ac38e5050900538506 -size 113395 +oid sha256:b734988bf539406d81d6488b72cfc27df649ccf461564721f215d1b6f12eb21c +size 112455 diff --git a/include/constfilt/discretize.hpp b/include/constfilt/discretize.hpp index 28d2d5c..acd3d96 100644 --- a/include/constfilt/discretize.hpp +++ b/include/constfilt/discretize.hpp @@ -204,37 +204,15 @@ constexpr StateSpace zoh_discretize(const StateSpace &sys_c, T Ts, // Fills monic characteristic polynomial of Ad: // [1, c_1, c_2, ..., c_N] (N+1 coefficients) -// Uses consteig::eigenvalues to obtain lam_1..lam_N, then builds -// (z - lam_1)(z - lam_2)...(z - lam_N) in complex arithmetic. -// Real parts are extracted at the end (imaginary parts cancel for real A). +// Uses the Faddeev-LeVerrier algorithm via consteig::char_poly: +// iterative matrix multiplications and traces, no eigenvalues. template constexpr void char_poly(const consteig::Matrix &Ad, T (&coeffs)[N + 1u]) { - using Complex = consteig::Complex; - - const auto evals = consteig::eigenvalues(Ad); // Matrix - - // p[0..k] holds the monic polynomial of degree k after k iterations. - Complex p[N + 1u]{}; - p[0] = Complex{static_cast(1), static_cast(0)}; - - for (consteig::Size k = 0; k < N; ++k) - { - const Complex lam = evals(k, 0); - // Multiply degree-k poly by (z - lam), working high-to-low in-place. - p[k + 1u] = Complex{static_cast(0), static_cast(0)} - lam * p[k]; - for (consteig::Size i = k; i > 0u; --i) - { - p[i] = p[i] - lam * p[i - 1u]; - } - // p[0] is unchanged (stays 1) - } - + const auto res = consteig::char_poly(Ad); for (consteig::Size i = 0; i <= N; ++i) - { - coeffs[i] = p[i].real; - } + coeffs[i] = res(i, 0u); } // Markov numerator diff --git a/include/constfilt/vendor/consteig/matrix/operations.hpp b/include/constfilt/vendor/consteig/matrix/operations.hpp index 27a6c3e..345e4eb 100644 --- a/include/constfilt/vendor/consteig/matrix/operations.hpp +++ b/include/constfilt/vendor/consteig/matrix/operations.hpp @@ -379,6 +379,36 @@ constexpr T trace(const Matrix &mat) return result; } +/// @brief Monic characteristic polynomial via the Faddeev-LeVerrier algorithm. +/// +/// Computes det(lam*I - A) = lam^N + c_1*lam^(N-1) + ... + c_N using only +/// matrix multiplications and traces; no eigenvalues, no complex arithmetic. +/// Susceptible to catastrophic cancellation for near-repeated eigenvalues. +/// +/// @tparam T Scalar type. +/// @tparam N Matrix dimension. +/// @param A Square NxN matrix. +/// @return Column vector of N+1 coefficients in descending power order, +/// with result(0,0) = 1 (monic leading term). +template +constexpr Matrix char_poly(const Matrix &A) +{ + Matrix coeffs{}; + coeffs(0u, 0u) = static_cast(1); + + Matrix M = eye(); + + for (Size k = 1u; k <= N; ++k) + { + const Matrix B = A * M; + const T ck = -trace(B) / static_cast(k); + coeffs(k, 0u) = ck; + M = B + ck * eye(); + } + + return coeffs; +} + /// @brief Element-wise approximate equality within an absolute tolerance. /// /// Returns `true` if every element satisfies From 9c44604522b661503adb672071a505b2a5850a37 Mon Sep 17 00:00:00 2001 From: Mitchell Thompkins Date: Thu, 25 Jun 2026 04:07:12 +0000 Subject: [PATCH 02/11] Revert "Evaluate Faddeev-LeVerrier char_poly via consteig, no accuracy gain" This reverts commit ed93c83a7ac559a6598b76f9e82b3ca6bc9e6996. --- .../profiling/results/accuracy_gcc_15.2.0.csv | 80 +++++++------- .../profiling/results/accuracy_gcc_15.2.0.png | 4 +- .../results/compile_times_gcc_15.2.0.csv | 102 +++++++++--------- .../results/compile_times_gcc_15.2.0.png | 4 +- docs/profiling/results/runtime_gcc_15.2.0.csv | 78 +++++++------- docs/profiling/results/runtime_gcc_15.2.0.png | 4 +- include/constfilt/discretize.hpp | 30 +++++- .../vendor/consteig/matrix/operations.hpp | 30 ------ 8 files changed, 162 insertions(+), 170 deletions(-) diff --git a/docs/profiling/results/accuracy_gcc_15.2.0.csv b/docs/profiling/results/accuracy_gcc_15.2.0.csv index 0e826b3..dcfad8b 100644 --- a/docs/profiling/results/accuracy_gcc_15.2.0.csv +++ b/docs/profiling/results/accuracy_gcc_15.2.0.csv @@ -4,16 +4,16 @@ library,filter_type,order,method,max_b_err,max_a_err,max_step_err constfilt,butterworth,1,zoh,0.000000e+00,0.000000e+00,0.000000e+00 constfilt,butterworth,2,zoh,3.608225e-16,2.220446e-16,1.110223e-15 -constfilt,butterworth,3,zoh,9.020562e-16,1.332268e-15,1.088019e-14 -constfilt,butterworth,4,zoh,8.971990e-15,1.110223e-15,1.321165e-14 -constfilt,butterworth,5,zoh,6.944098e-15,7.549517e-15,1.543210e-14 -constfilt,butterworth,6,zoh,3.286173e-14,9.769963e-15,5.095924e-14 -constfilt,butterworth,7,zoh,1.100198e-11,2.131628e-14,1.254094e-12 -constfilt,butterworth,8,zoh,2.955275e-08,4.849454e-13,2.134253e-09 -constfilt,butterworth,9,zoh,7.577398e-05,8.100187e-13,3.692378e-06 -constfilt,butterworth,10,zoh,2.065719e-01,2.042100e-11,6.461003e-03 -constfilt,butterworth,11,zoh,2.543305e+02,3.731445e-07,4.685698e+00 -constfilt,butterworth,12,zoh,1.799631e+03,3.914806e-10,2.955113e+02 +constfilt,butterworth,3,zoh,9.020562e-16,1.776357e-15,7.771561e-15 +constfilt,butterworth,4,zoh,8.971990e-15,8.881784e-16,1.587619e-14 +constfilt,butterworth,5,zoh,6.945833e-15,3.552714e-15,1.809664e-14 +constfilt,butterworth,6,zoh,3.284352e-14,1.598721e-14,6.572520e-14 +constfilt,butterworth,7,zoh,1.100200e-11,4.796163e-14,1.254094e-12 +constfilt,butterworth,8,zoh,2.955275e-08,5.062617e-13,2.134253e-09 +constfilt,butterworth,9,zoh,7.577398e-05,8.526513e-13,3.692378e-06 +constfilt,butterworth,10,zoh,2.065719e-01,2.030731e-11,6.461003e-03 +constfilt,butterworth,11,zoh,2.543305e+02,3.566485e-08,4.685698e+00 +constfilt,butterworth,12,zoh,1.799631e+03,3.915801e-10,2.955113e+02 constfilt,butterworth,1,matchedz,5.551115e-17,0.000000e+00,2.220446e-16 constfilt,butterworth,2,matchedz,2.775558e-17,2.220446e-16,1.110223e-15 constfilt,butterworth,3,matchedz,1.387779e-17,1.776357e-15,1.998401e-15 @@ -27,28 +27,28 @@ constfilt,butterworth,10,matchedz,1.523304e-16,5.579892e-11,1.363021e-11 constfilt,butterworth,11,matchedz,2.168404e-19,1.847411e-13,4.028200e-11 constfilt,butterworth,12,matchedz,9.662952e-17,3.253149e-10,1.482581e-10 constfilt,butterworth,1,prewarp,8.326673e-17,2.220446e-16,2.220446e-16 -constfilt,butterworth,2,prewarp,1.249001e-16,5.551115e-17,1.998401e-15 -constfilt,butterworth,3,prewarp,5.204170e-17,4.440892e-16,3.330669e-15 -constfilt,butterworth,4,prewarp,1.734723e-16,2.220446e-15,1.598721e-14 -constfilt,butterworth,5,prewarp,1.682682e-16,2.220446e-15,2.131628e-14 -constfilt,butterworth,6,prewarp,1.214306e-16,1.065814e-14,3.508305e-14 -constfilt,butterworth,7,prewarp,2.474149e-16,7.993606e-15,2.481348e-13 -constfilt,butterworth,8,prewarp,5.555452e-16,3.019807e-14,7.693846e-13 -constfilt,butterworth,9,prewarp,4.776995e-16,2.131628e-14,4.931056e-12 -constfilt,butterworth,10,prewarp,3.020696e-15,3.552714e-14,4.900191e-12 -constfilt,butterworth,11,prewarp,1.400681e-15,6.394885e-14,3.339040e-11 -constfilt,butterworth,12,prewarp,1.350855e-14,1.705303e-13,3.085243e-11 +constfilt,butterworth,2,prewarp,1.249001e-16,0.000000e+00,2.331468e-15 +constfilt,butterworth,3,prewarp,1.283695e-16,2.442491e-15,1.776357e-15 +constfilt,butterworth,4,prewarp,1.630640e-16,3.108624e-15,1.432188e-14 +constfilt,butterworth,5,prewarp,1.647987e-16,8.881784e-16,2.930989e-14 +constfilt,butterworth,6,prewarp,1.231654e-16,1.154632e-14,1.159073e-13 +constfilt,butterworth,7,prewarp,2.493665e-16,2.664535e-14,1.717515e-13 +constfilt,butterworth,8,prewarp,5.490400e-16,5.684342e-14,5.651035e-13 +constfilt,butterworth,9,prewarp,4.832289e-16,2.842171e-14,4.035883e-12 +constfilt,butterworth,10,prewarp,3.041296e-15,2.486900e-14,6.625700e-12 +constfilt,butterworth,11,prewarp,1.426024e-15,3.197442e-13,1.841038e-11 +constfilt,butterworth,12,prewarp,1.350920e-14,1.563194e-13,5.809664e-11 constfilt,elliptic,2,zoh,2.395861e-13,2.842171e-14,1.409099e-13 -constfilt,elliptic,3,zoh,5.399570e-13,9.481305e-14,3.640976e-13 -constfilt,elliptic,4,zoh,5.392874e-13,1.083578e-13,8.994472e-13 -constfilt,elliptic,5,zoh,1.066758e-12,4.858336e-13,1.475042e-12 -constfilt,elliptic,6,zoh,1.760855e-12,1.453060e-12,1.816436e-12 -constfilt,elliptic,7,zoh,1.055522e-11,1.797673e-11,9.918288e-12 -constfilt,elliptic,8,zoh,2.154666e-09,2.641208e-09,7.933643e-10 -constfilt,elliptic,9,zoh,3.582365e-07,1.304601e-07,9.761215e-07 -constfilt,elliptic,10,zoh,3.461312e-03,3.213907e-06,6.241612e-03 -constfilt,elliptic,11,zoh,5.384745e-01,4.802802e-05,7.545471e-02 -constfilt,elliptic,12,zoh,5.474468e+02,5.003185e-04,2.311052e+01 +constfilt,elliptic,3,zoh,5.400819e-13,9.858780e-14,3.640976e-13 +constfilt,elliptic,4,zoh,5.393082e-13,1.074696e-13,8.994472e-13 +constfilt,elliptic,5,zoh,1.066647e-12,4.938272e-13,1.474598e-12 +constfilt,elliptic,6,zoh,1.760619e-12,1.481482e-12,1.819322e-12 +constfilt,elliptic,7,zoh,1.055545e-11,1.796252e-11,9.916956e-12 +constfilt,elliptic,8,zoh,2.154665e-09,2.641301e-09,7.933327e-10 +constfilt,elliptic,9,zoh,3.582365e-07,1.304602e-07,9.765907e-07 +constfilt,elliptic,10,zoh,3.461312e-03,3.213908e-06,6.243202e-03 +constfilt,elliptic,11,zoh,5.384745e-01,4.802802e-05,7.544774e-02 +constfilt,elliptic,12,zoh,5.474468e+02,5.003185e-04,2.311049e+01 constfilt,elliptic,2,matchedz,3.854481e-11,2.786660e-14,1.928402e-11 constfilt,elliptic,3,matchedz,5.143178e-13,9.592327e-14,3.533007e-13 constfilt,elliptic,4,matchedz,1.324132e-13,1.101341e-13,8.861800e-13 @@ -61,15 +61,15 @@ constfilt,elliptic,10,matchedz,7.442555e-07,3.213346e-06,2.409790e-07 constfilt,elliptic,11,matchedz,2.360767e-05,4.802825e-05,1.950004e-06 constfilt,elliptic,12,matchedz,1.156955e-04,5.002874e-04,1.090689e-05 constfilt,elliptic,2,prewarp,1.637024e-13,3.286260e-14,9.213463e-14 -constfilt,elliptic,3,prewarp,1.595737e-13,8.104628e-14,3.510803e-13 -constfilt,elliptic,4,prewarp,3.788705e-13,1.092459e-13,8.841816e-13 -constfilt,elliptic,5,prewarp,2.869857e-13,4.902745e-13,1.438516e-12 -constfilt,elliptic,6,prewarp,1.340678e-12,1.428191e-12,1.799783e-12 -constfilt,elliptic,7,prewarp,2.617934e-12,1.758949e-11,9.648282e-12 -constfilt,elliptic,8,prewarp,1.244404e-09,2.586262e-09,7.183454e-10 -constfilt,elliptic,9,prewarp,1.482448e-08,1.278829e-07,1.805562e-08 -constfilt,elliptic,10,prewarp,1.376367e-06,3.130828e-06,2.329530e-07 -constfilt,elliptic,11,prewarp,4.570721e-06,4.696437e-05,1.897979e-06 +constfilt,elliptic,3,prewarp,1.596015e-13,8.082424e-14,3.510803e-13 +constfilt,elliptic,4,prewarp,3.788428e-13,1.079137e-13,8.840706e-13 +constfilt,elliptic,5,prewarp,2.870135e-13,4.884981e-13,1.437517e-12 +constfilt,elliptic,6,prewarp,1.340761e-12,1.419309e-12,1.797784e-12 +constfilt,elliptic,7,prewarp,2.618128e-12,1.753975e-11,9.648726e-12 +constfilt,elliptic,8,prewarp,1.244403e-09,2.586319e-09,7.183503e-10 +constfilt,elliptic,9,prewarp,1.482448e-08,1.278831e-07,1.805560e-08 +constfilt,elliptic,10,prewarp,1.376367e-06,3.130827e-06,2.329530e-07 +constfilt,elliptic,11,prewarp,4.570721e-06,4.696437e-05,1.897978e-06 constfilt,elliptic,12,prewarp,2.111470e-04,4.866247e-04,1.055330e-05 iir1,butterworth,1,prewarp,n/a,n/a,4.440892e-16 iir1,butterworth,2,prewarp,n/a,n/a,8.881784e-16 diff --git a/docs/profiling/results/accuracy_gcc_15.2.0.png b/docs/profiling/results/accuracy_gcc_15.2.0.png index f63dbb9..996c68d 100644 --- a/docs/profiling/results/accuracy_gcc_15.2.0.png +++ b/docs/profiling/results/accuracy_gcc_15.2.0.png @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:89d68110f966dd40c7c502e1c38033ce69e54d7985522432ef624ad7d8e9af2a -size 169078 +oid sha256:3063f67ee8c7ca693c9bc690f2bf9e9e3a35e86b8025f1938552ae6ef2271908 +size 169290 diff --git a/docs/profiling/results/compile_times_gcc_15.2.0.csv b/docs/profiling/results/compile_times_gcc_15.2.0.csv index 183c566..d7f20e2 100644 --- a/docs/profiling/results/compile_times_gcc_15.2.0.csv +++ b/docs/profiling/results/compile_times_gcc_15.2.0.csv @@ -2,54 +2,54 @@ # compiler: g++ (Alpine 15.2.0) 15.2.0 # os: Alpine Linux filter_type,order,method,compile_time_sec,max_rss_kb,exit_code -butterworth,10,matchedz,0.37,72924,0 -butterworth,10,tustin,0.26,60356,0 -butterworth,10,zoh,0.89,139956,0 -butterworth,11,matchedz,0.49,88088,0 -butterworth,11,tustin,0.34,70304,0 -butterworth,11,zoh,1.21,177720,0 -butterworth,12,matchedz,0.64,105620,0 -butterworth,12,tustin,0.45,83064,0 -butterworth,12,zoh,1.59,225348,0 -butterworth,1,matchedz,0.07,36864,0 -butterworth,1,tustin,0.06,36184,0 -butterworth,1,zoh,0.07,37032,0 -butterworth,2,matchedz,0.07,37020,0 -butterworth,2,tustin,0.06,36620,0 -butterworth,2,zoh,0.08,38152,0 -butterworth,4,matchedz,0.09,39224,0 -butterworth,4,tustin,0.07,37548,0 -butterworth,4,zoh,0.12,44124,0 -butterworth,6,matchedz,0.13,44296,0 -butterworth,6,tustin,0.10,40444,0 -butterworth,6,zoh,0.25,58156,0 -butterworth,8,matchedz,0.21,54336,0 -butterworth,8,tustin,0.15,47184,0 -butterworth,8,zoh,0.48,86952,0 -butterworth,9,matchedz,0.29,62732,0 -butterworth,9,tustin,0.19,52700,0 -butterworth,9,zoh,0.67,109980,0 -elliptic,10,matchedz,0.73,112956,0 -elliptic,10,tustin,0.33,66488,0 -elliptic,10,zoh,0.98,145432,0 -elliptic,11,matchedz,0.93,140208,0 -elliptic,11,tustin,0.41,75984,0 -elliptic,11,zoh,1.29,184024,0 -elliptic,12,matchedz,1.30,181384,0 -elliptic,12,tustin,0.54,89480,0 -elliptic,12,zoh,1.68,232768,0 -elliptic,2,matchedz,0.10,40604,0 -elliptic,2,tustin,0.09,40108,0 -elliptic,2,zoh,0.11,42152,0 -elliptic,4,matchedz,0.15,45480,0 -elliptic,4,tustin,0.12,41884,0 -elliptic,4,zoh,0.17,48260,0 -elliptic,6,matchedz,0.22,53660,0 -elliptic,6,tustin,0.15,45224,0 -elliptic,6,zoh,0.29,62044,0 -elliptic,8,matchedz,0.41,74268,0 -elliptic,8,tustin,0.21,52560,0 -elliptic,8,zoh,0.53,91808,0 -elliptic,9,matchedz,0.52,89608,0 -elliptic,9,tustin,0.26,58540,0 -elliptic,9,zoh,0.73,115224,0 +butterworth,10,matchedz,0.37,73000,0 +butterworth,10,tustin,0.45,83112,0 +butterworth,10,zoh,1.07,158552,0 +butterworth,11,matchedz,0.50,88108,0 +butterworth,11,tustin,0.62,102420,0 +butterworth,11,zoh,1.42,203688,0 +butterworth,12,matchedz,0.63,105460,0 +butterworth,12,tustin,0.81,126272,0 +butterworth,12,zoh,1.93,264624,0 +butterworth,1,matchedz,0.07,36764,0 +butterworth,1,tustin,0.06,36776,0 +butterworth,1,zoh,0.07,37024,0 +butterworth,2,matchedz,0.07,37096,0 +butterworth,2,tustin,0.06,37268,0 +butterworth,2,zoh,0.07,38084,0 +butterworth,4,matchedz,0.09,39528,0 +butterworth,4,tustin,0.09,40116,0 +butterworth,4,zoh,0.14,44948,0 +butterworth,6,matchedz,0.12,44468,0 +butterworth,6,tustin,0.14,45600,0 +butterworth,6,zoh,0.27,60956,0 +butterworth,8,matchedz,0.21,54360,0 +butterworth,8,tustin,0.25,58428,0 +butterworth,8,zoh,0.54,94996,0 +butterworth,9,matchedz,0.28,62652,0 +butterworth,9,tustin,0.34,69284,0 +butterworth,9,zoh,0.77,123156,0 +elliptic,10,matchedz,0.72,112964,0 +elliptic,10,tustin,0.53,89548,0 +elliptic,10,zoh,1.13,163684,0 +elliptic,11,matchedz,0.95,140320,0 +elliptic,11,tustin,0.69,109080,0 +elliptic,11,zoh,1.52,210084,0 +elliptic,12,matchedz,1.30,181428,0 +elliptic,12,tustin,0.92,134184,0 +elliptic,12,zoh,2.03,271788,0 +elliptic,2,matchedz,0.11,40780,0 +elliptic,2,tustin,0.10,40608,0 +elliptic,2,zoh,0.11,41888,0 +elliptic,4,matchedz,0.14,45328,0 +elliptic,4,tustin,0.13,43936,0 +elliptic,4,zoh,0.18,49076,0 +elliptic,6,matchedz,0.22,53920,0 +elliptic,6,tustin,0.20,50064,0 +elliptic,6,zoh,0.31,64840,0 +elliptic,8,matchedz,0.40,74308,0 +elliptic,8,tustin,0.31,63908,0 +elliptic,8,zoh,0.61,100160,0 +elliptic,9,matchedz,0.53,89588,0 +elliptic,9,tustin,0.40,74688,0 +elliptic,9,zoh,0.84,128436,0 diff --git a/docs/profiling/results/compile_times_gcc_15.2.0.png b/docs/profiling/results/compile_times_gcc_15.2.0.png index d2dce27..0a1b52f 100644 --- a/docs/profiling/results/compile_times_gcc_15.2.0.png +++ b/docs/profiling/results/compile_times_gcc_15.2.0.png @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:b2d38127103dd59828b20e490ebe7f6c1b1f2ca00e8e3be08d90c51f02bf28fc -size 147973 +oid sha256:4a4f746ddc0d9c0f7f9380a8c7fe02c44288304df30b974042213c7456a1af73 +size 150965 diff --git a/docs/profiling/results/runtime_gcc_15.2.0.csv b/docs/profiling/results/runtime_gcc_15.2.0.csv index 35a0a76..5b7cc99 100644 --- a/docs/profiling/results/runtime_gcc_15.2.0.csv +++ b/docs/profiling/results/runtime_gcc_15.2.0.csv @@ -2,42 +2,42 @@ # compiler: g++ (Alpine 15.2.0) 15.2.0 # os: Alpine Linux library,filter_type,order,method,ns_per_sample,msa_per_s,dc_gain -constfilt,Butterworth,1,tustin,6.418,155.8,1.00000000 -constfilt,Butterworth,2,tustin,10.855,92.1,1.00000000 -constfilt,Butterworth,4,tustin,18.741,53.4,1.00000000 -constfilt,Butterworth,6,tustin,26.035,38.4,1.00000000 -constfilt,Butterworth,8,tustin,34.363,29.1,1.00000000 -constfilt,Butterworth,1,zoh,6.486,154.2,1.00000000 -constfilt,Butterworth,2,zoh,10.817,92.4,1.00000000 -constfilt,Butterworth,4,zoh,18.597,53.8,1.00000000 -constfilt,Butterworth,6,zoh,25.921,38.6,1.00000000 -constfilt,Butterworth,8,zoh,34.208,29.2,1.00000000 -constfilt,Butterworth,1,matchedz,6.421,155.7,1.00000000 -constfilt,Butterworth,2,matchedz,10.851,92.2,1.00000000 -constfilt,Butterworth,4,matchedz,18.754,53.3,1.00000000 -constfilt,Butterworth,6,matchedz,26.121,38.3,1.00000000 -constfilt,Butterworth,8,matchedz,34.793,28.7,1.00000000 -constfilt,Elliptic,2,tustin,10.829,92.3,0.94406088 -constfilt,Elliptic,4,tustin,18.308,54.6,0.94406088 -constfilt,Elliptic,6,tustin,26.122,38.3,0.94406088 -constfilt,Elliptic,8,tustin,34.242,29.2,0.94406088 -constfilt,Elliptic,2,zoh,10.817,92.5,0.94406088 -constfilt,Elliptic,4,zoh,18.313,54.6,0.94406088 -constfilt,Elliptic,6,zoh,26.092,38.3,0.94406088 -constfilt,Elliptic,8,zoh,34.251,29.2,0.94406088 -constfilt,Elliptic,2,matchedz,11.076,90.3,0.94406088 -constfilt,Elliptic,4,matchedz,18.853,53.0,0.94406088 -constfilt,Elliptic,6,matchedz,25.973,38.5,0.94406088 -constfilt,Elliptic,8,matchedz,34.539,29.0,0.94406088 -iir1,Butterworth,2,runtime,10.830,92.3,1.00000000 -iir1,Butterworth,4,runtime,18.134,55.1,1.00000000 -iir1,Butterworth,6,runtime,29.181,34.3,1.00000000 -iir1,Butterworth,8,runtime,40.300,24.8,1.00000000 -kfr,Butterworth,2,runtime+simd,588.824,1.7,1.00000000 -kfr,Butterworth,4,runtime+simd,381.513,2.6,1.00000000 -kfr,Butterworth,6,runtime+simd,297.600,3.4,1.00000000 -kfr,Butterworth,8,runtime+simd,296.401,3.4,1.00000000 -kfr,Elliptic,2,runtime+simd,589.859,1.7,0.94406088 -kfr,Elliptic,4,runtime+simd,382.360,2.6,0.94406088 -kfr,Elliptic,6,runtime+simd,296.054,3.4,0.94406088 -kfr,Elliptic,8,runtime+simd,297.250,3.4,0.94406088 +constfilt,Butterworth,1,tustin,6.415,155.9,1.00000000 +constfilt,Butterworth,2,tustin,10.876,91.9,1.00000000 +constfilt,Butterworth,4,tustin,18.954,52.8,1.00000000 +constfilt,Butterworth,6,tustin,26.164,38.2,1.00000000 +constfilt,Butterworth,8,tustin,34.550,28.9,1.00000000 +constfilt,Butterworth,1,zoh,6.420,155.8,1.00000000 +constfilt,Butterworth,2,zoh,10.826,92.4,1.00000000 +constfilt,Butterworth,4,zoh,18.535,54.0,1.00000000 +constfilt,Butterworth,6,zoh,26.066,38.4,1.00000000 +constfilt,Butterworth,8,zoh,34.415,29.1,1.00000000 +constfilt,Butterworth,1,matchedz,6.423,155.7,1.00000000 +constfilt,Butterworth,2,matchedz,10.839,92.3,1.00000000 +constfilt,Butterworth,4,matchedz,18.858,53.0,1.00000000 +constfilt,Butterworth,6,matchedz,26.130,38.3,1.00000000 +constfilt,Butterworth,8,matchedz,34.357,29.1,1.00000000 +constfilt,Elliptic,2,tustin,10.844,92.2,0.94406088 +constfilt,Elliptic,4,tustin,18.373,54.4,0.94406088 +constfilt,Elliptic,6,tustin,26.217,38.1,0.94406088 +constfilt,Elliptic,8,tustin,34.485,29.0,0.94406088 +constfilt,Elliptic,2,zoh,10.828,92.4,0.94406088 +constfilt,Elliptic,4,zoh,18.359,54.5,0.94406088 +constfilt,Elliptic,6,zoh,26.095,38.3,0.94406088 +constfilt,Elliptic,8,zoh,34.249,29.2,0.94406088 +constfilt,Elliptic,2,matchedz,11.091,90.2,0.94406088 +constfilt,Elliptic,4,matchedz,18.904,52.9,0.94406088 +constfilt,Elliptic,6,matchedz,26.024,38.4,0.94406088 +constfilt,Elliptic,8,matchedz,34.507,29.0,0.94406088 +iir1,Butterworth,2,runtime,10.732,93.2,1.00000000 +iir1,Butterworth,4,runtime,18.155,55.1,1.00000000 +iir1,Butterworth,6,runtime,30.092,33.2,1.00000000 +iir1,Butterworth,8,runtime,40.322,24.8,1.00000000 +kfr,Butterworth,2,runtime+simd,607.261,1.6,1.00000000 +kfr,Butterworth,4,runtime+simd,402.769,2.5,1.00000000 +kfr,Butterworth,6,runtime+simd,307.412,3.3,1.00000000 +kfr,Butterworth,8,runtime+simd,307.261,3.3,1.00000000 +kfr,Elliptic,2,runtime+simd,605.280,1.7,0.94406088 +kfr,Elliptic,4,runtime+simd,401.691,2.5,0.94406088 +kfr,Elliptic,6,runtime+simd,306.657,3.3,0.94406088 +kfr,Elliptic,8,runtime+simd,308.208,3.2,0.94406088 diff --git a/docs/profiling/results/runtime_gcc_15.2.0.png b/docs/profiling/results/runtime_gcc_15.2.0.png index 3627dcb..931a63a 100644 --- a/docs/profiling/results/runtime_gcc_15.2.0.png +++ b/docs/profiling/results/runtime_gcc_15.2.0.png @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:b734988bf539406d81d6488b72cfc27df649ccf461564721f215d1b6f12eb21c -size 112455 +oid sha256:e23c6d526144320c1c2baf0928589539fa06e92946d213ac38e5050900538506 +size 113395 diff --git a/include/constfilt/discretize.hpp b/include/constfilt/discretize.hpp index acd3d96..28d2d5c 100644 --- a/include/constfilt/discretize.hpp +++ b/include/constfilt/discretize.hpp @@ -204,15 +204,37 @@ constexpr StateSpace zoh_discretize(const StateSpace &sys_c, T Ts, // Fills monic characteristic polynomial of Ad: // [1, c_1, c_2, ..., c_N] (N+1 coefficients) -// Uses the Faddeev-LeVerrier algorithm via consteig::char_poly: -// iterative matrix multiplications and traces, no eigenvalues. +// Uses consteig::eigenvalues to obtain lam_1..lam_N, then builds +// (z - lam_1)(z - lam_2)...(z - lam_N) in complex arithmetic. +// Real parts are extracted at the end (imaginary parts cancel for real A). template constexpr void char_poly(const consteig::Matrix &Ad, T (&coeffs)[N + 1u]) { - const auto res = consteig::char_poly(Ad); + using Complex = consteig::Complex; + + const auto evals = consteig::eigenvalues(Ad); // Matrix + + // p[0..k] holds the monic polynomial of degree k after k iterations. + Complex p[N + 1u]{}; + p[0] = Complex{static_cast(1), static_cast(0)}; + + for (consteig::Size k = 0; k < N; ++k) + { + const Complex lam = evals(k, 0); + // Multiply degree-k poly by (z - lam), working high-to-low in-place. + p[k + 1u] = Complex{static_cast(0), static_cast(0)} - lam * p[k]; + for (consteig::Size i = k; i > 0u; --i) + { + p[i] = p[i] - lam * p[i - 1u]; + } + // p[0] is unchanged (stays 1) + } + for (consteig::Size i = 0; i <= N; ++i) - coeffs[i] = res(i, 0u); + { + coeffs[i] = p[i].real; + } } // Markov numerator diff --git a/include/constfilt/vendor/consteig/matrix/operations.hpp b/include/constfilt/vendor/consteig/matrix/operations.hpp index 345e4eb..27a6c3e 100644 --- a/include/constfilt/vendor/consteig/matrix/operations.hpp +++ b/include/constfilt/vendor/consteig/matrix/operations.hpp @@ -379,36 +379,6 @@ constexpr T trace(const Matrix &mat) return result; } -/// @brief Monic characteristic polynomial via the Faddeev-LeVerrier algorithm. -/// -/// Computes det(lam*I - A) = lam^N + c_1*lam^(N-1) + ... + c_N using only -/// matrix multiplications and traces; no eigenvalues, no complex arithmetic. -/// Susceptible to catastrophic cancellation for near-repeated eigenvalues. -/// -/// @tparam T Scalar type. -/// @tparam N Matrix dimension. -/// @param A Square NxN matrix. -/// @return Column vector of N+1 coefficients in descending power order, -/// with result(0,0) = 1 (monic leading term). -template -constexpr Matrix char_poly(const Matrix &A) -{ - Matrix coeffs{}; - coeffs(0u, 0u) = static_cast(1); - - Matrix M = eye(); - - for (Size k = 1u; k <= N; ++k) - { - const Matrix B = A * M; - const T ck = -trace(B) / static_cast(k); - coeffs(k, 0u) = ck; - M = B + ck * eye(); - } - - return coeffs; -} - /// @brief Element-wise approximate equality within an absolute tolerance. /// /// Returns `true` if every element satisfies From 0e873f20a91bc1a87df168925a932b2f3be497df Mon Sep 17 00:00:00 2001 From: Mitchell Thompkins Date: Thu, 25 Jun 2026 05:41:08 +0000 Subject: [PATCH 03/11] Thread analytic poles/zeros through ZOH and Matched-Z discretization Butterworth and Elliptic now carry their analytically-computed poles and zeros as a FactoredTF through to the discretization step. ZOH skips the QR eigenvalue search on the CCF matrix; Matched-Z skips both companion matrix eigendecompositions entirely. Tustin is unchanged. Adds FactoredTF struct, matrix_exp_with_evals, zoh_discretize_with_evals, matched_z_discretize_factored, and discretize_with_factored dispatch overloads to discretize.hpp. --- include/constfilt/analog_filter.hpp | 18 +++ include/constfilt/butterworth.hpp | 39 ++++- include/constfilt/discretize.hpp | 242 ++++++++++++++++++++++++++++ include/constfilt/elliptic.hpp | 197 ++++++++++++++++------ 4 files changed, 450 insertions(+), 46 deletions(-) diff --git a/include/constfilt/analog_filter.hpp b/include/constfilt/analog_filter.hpp index 11b19e9..204895e 100644 --- a/include/constfilt/analog_filter.hpp +++ b/include/constfilt/analog_filter.hpp @@ -79,6 +79,15 @@ class AnalogFilter : public Filter { } + constexpr AnalogFilter(TransferFunction continuous_tf, + const FactoredTF &ftf, T sample_rate_hz, + BoundMethod method_tag) + : AnalogFilter(checked_discretize_factored(continuous_tf.b, + continuous_tf.a, ftf, + sample_rate_hz, method_tag)) + { + } + private: constexpr explicit AnalogFilter( TransferFunction digital_tf) @@ -120,6 +129,15 @@ class AnalogFilter : public Filter return analog_to_digital( b_c, a_c, static_cast(1) / sample_rate_hz, method_tag); } + + static constexpr TransferFunction + checked_discretize_factored(const T (&b_c)[N + 1u], const T (&a_c)[N + 1u], + const FactoredTF &ftf, T sample_rate_hz, + BoundMethod method_tag) + { + return discretize_with_factored( + b_c, a_c, ftf, static_cast(1) / sample_rate_hz, method_tag); + } }; } // namespace constfilt diff --git a/include/constfilt/butterworth.hpp b/include/constfilt/butterworth.hpp index b50d4e5..88304c4 100644 --- a/include/constfilt/butterworth.hpp +++ b/include/constfilt/butterworth.hpp @@ -34,7 +34,8 @@ class Butterworth // Construct from filter specification; all math is constexpr. constexpr Butterworth(T cutoff_hz, T sample_rate_hz) : AnalogFilter( - compute_continuous_tf(cutoff_hz), sample_rate_hz, + compute_continuous_tf(cutoff_hz), + compute_factored_tf(cutoff_hz, FilterType{}), sample_rate_hz, make_tustin_tag(cutoff_hz, BoundMethod{})) { } @@ -91,6 +92,42 @@ class Butterworth } } + // LP: poles at wc*exp(j*theta_k), no finite zeros, gain = wc^N. + static constexpr FactoredTF compute_factored_tf(T cutoff_hz, LowPass) + { + const T wc = static_cast(2) * static_cast(GCEM_PI) * cutoff_hz; + FactoredTF ftf{}; + ftf.nz = 0; + ftf.gain = gcem::pow(wc, static_cast(N)); + for (consteig::Size k = 1u; k <= N; ++k) + { + const T theta = static_cast(GCEM_PI) * + static_cast(2u * k + N - 1u) / + static_cast(2u * N); + ftf.poles[k - 1u] = {wc * gcem::cos(theta), wc * gcem::sin(theta)}; + } + return ftf; + } + + // HP: poles at wc*exp(-j*theta_k) (magnitude wc), N zeros at s=0, gain=1. + static constexpr FactoredTF compute_factored_tf(T cutoff_hz, HighPass) + { + using Complex = consteig::Complex; + const T wc = static_cast(2) * static_cast(GCEM_PI) * cutoff_hz; + FactoredTF ftf{}; + ftf.nz = N; + ftf.gain = static_cast(1); + for (consteig::Size k = 1u; k <= N; ++k) + { + const T theta = static_cast(GCEM_PI) * + static_cast(2u * k + N - 1u) / + static_cast(2u * N); + ftf.poles[k - 1u] = {wc * gcem::cos(theta), -wc * gcem::sin(theta)}; + ftf.zeros[k - 1u] = Complex{static_cast(0), static_cast(0)}; + } + return ftf; + } + // Normalized Butterworth denominator coefficients (wc=1, monic). // Fills result in descending power order: [1, p[N-1], ..., p[0]] // diff --git a/include/constfilt/discretize.hpp b/include/constfilt/discretize.hpp index 28d2d5c..31a59df 100644 --- a/include/constfilt/discretize.hpp +++ b/include/constfilt/discretize.hpp @@ -103,6 +103,20 @@ struct TransferFunction T a[NA]{}; }; +// An analog transfer function in factored (pole/zero/gain) form. +// +// H(s) = gain * prod_i(s - zeros[i]) / prod_j(s - poles[j]) +// +// Only the first `nz` entries of zeros[] are meaningful; the rest are +// zero-initialized and must not be read by callers. +template struct FactoredTF +{ + consteig::Complex poles[N]{}; // all N analog poles + consteig::Complex zeros[N]{}; // first nz finite analog zeros + consteig::Size nz{0}; // number of populated zeros entries + T gain{}; // k_c = b[d_b] / a[0] +}; + // Matrix exponential // matrix_exp(A) via eigendecomposition: @@ -164,6 +178,49 @@ constexpr consteig::Matrix matrix_exp( return result; } +// matrix_exp variant that accepts precomputed eigenvalues of A. +// Identical to matrix_exp except consteig::eigenvalues(A) is skipped; +// the caller supplies the eigenvalues directly via `evals`. +// For ZOH use: pass (Ts*Ac, scaled_evals) where scaled_evals(i,0) = Ts*pole_i. +template +constexpr consteig::Matrix matrix_exp_with_evals( + const consteig::Matrix &A, + const consteig::Matrix, N, 1> &evals) +{ + using Complex = consteig::Complex; + using ComplexMat_NN = consteig::Matrix; + using ComplexMat_N1 = consteig::Matrix; + + const auto V = consteig::eigenvectors(A, evals); + + const auto lu_V = consteig::lu(V); + ComplexMat_NN V_inv{}; + for (consteig::Size col = 0; col < N; ++col) + { + ComplexMat_N1 e_col{}; + e_col(col, 0) = Complex{static_cast(1), static_cast(0)}; + auto col_vec = consteig::lu_solve(lu_V, e_col); + for (consteig::Size row = 0; row < N; ++row) + V_inv(row, col) = col_vec(row, 0); + } + + ComplexMat_NN result_c{}; + for (consteig::Size i = 0; i < N; ++i) + { + Complex exp_lambda = consteig::exp(evals(i, 0)); + for (consteig::Size r = 0; r < N; ++r) + for (consteig::Size c = 0; c < N; ++c) + result_c(r, c) = + result_c(r, c) + exp_lambda * V(r, i) * V_inv(i, c); + } + + consteig::Matrix result{}; + for (consteig::Size r = 0; r < N; ++r) + for (consteig::Size c = 0; c < N; ++c) + result(r, c) = result_c(r, c).real; + return result; +} + // ZOH discretization // ZOH: Ad = matrix_exp(Ac*Ts), Bd = Ac^{-1} * (Ad - I) * Bc @@ -200,6 +257,38 @@ constexpr StateSpace zoh_discretize(const StateSpace &sys_c, T Ts, return sys_d; } +// ZOH discretization using analytically known poles from a FactoredTF. +// Avoids the QR eigenvalue search inside matrix_exp by supplying +// the scaled analog poles (eigenvalues of Ts*Ac) directly. +template +constexpr StateSpace zoh_discretize_with_evals( + const StateSpace &sys_c, T Ts, const FactoredTF &ftf) +{ + using Complex = consteig::Complex; + + const consteig::Matrix TsAc = Ts * sys_c.A; + + consteig::Matrix scaled_evals{}; + for (consteig::Size i = 0; i < N; ++i) + scaled_evals(i, 0) = + Complex{Ts * ftf.poles[i].real, Ts * ftf.poles[i].imag}; + + const consteig::Matrix Ad = + matrix_exp_with_evals(TsAc, scaled_evals); + + const consteig::Matrix AdmI = Ad - consteig::eye(); + const consteig::Matrix rhs = AdmI * sys_c.B; + const auto lu_Ac = consteig::lu(sys_c.A); + const auto Bd = consteig::lu_solve(lu_Ac, rhs); + + StateSpace sys_d{}; + sys_d.A = Ad; + sys_d.B = Bd; + sys_d.C = sys_c.C; + sys_d.D = sys_c.D; + return sys_d; +} + // Characteristic polynomial // Fills monic characteristic polynomial of Ad: @@ -576,6 +665,130 @@ constexpr TransferFunction matched_z_discretize_tf( return tf; } +// Matched-Z discretization from a FactoredTF (precomputed poles/zeros). +// Bypasses both companion-matrix eigendecompositions in +// matched_z_discretize_tf. Uses ftf.poles, ftf.zeros[0..ftf.nz-1], and ftf.gain +// directly. +template +constexpr TransferFunction matched_z_discretize_factored( + const FactoredTF &ftf, T Ts) +{ + using Complex = consteig::Complex; + + const consteig::Size nz = ftf.nz; + const T k_c = ftf.gain; + + // Map poles to z-domain; build monic denominator polynomial. + Complex p_d_vals[N]{}; + Complex pole_poly[N + 1u]{}; + pole_poly[0] = Complex{static_cast(1), static_cast(0)}; + for (consteig::Size k = 0; k < N; ++k) + { + p_d_vals[k] = + consteig::exp(ftf.poles[k] * Complex{Ts, static_cast(0)}); + const Complex zk = p_d_vals[k]; + pole_poly[k + 1u] = + Complex{static_cast(0), static_cast(0)} - zk * pole_poly[k]; + for (consteig::Size i = k; i > 0u; --i) + pole_poly[i] = pole_poly[i] - zk * pole_poly[i - 1u]; + } + + // Map finite zeros to z-domain; build monic numerator polynomial. + // No companion matrix or spurious-eigenvalue filtering needed. + Complex z_d_finite[N]{}; + Complex zero_poly[N + 1u]{}; + zero_poly[0] = Complex{static_cast(1), static_cast(0)}; + for (consteig::Size k = 0; k < nz; ++k) + { + z_d_finite[k] = + consteig::exp(ftf.zeros[k] * Complex{Ts, static_cast(0)}); + const Complex zk = z_d_finite[k]; + zero_poly[k + 1u] = + Complex{static_cast(0), static_cast(0)} - zk * zero_poly[k]; + for (consteig::Size i = k; i > 0u; --i) + zero_poly[i] = zero_poly[i] - zk * zero_poly[i - 1u]; + } + + // Pad with zeros at z = -1 to reach numerator degree N-1. + const consteig::Size n_extra = (nz + 1u < N) ? (N - nz - 1u) : 0u; + for (consteig::Size e = 0; e < n_extra; ++e) + { + const consteig::Size cur_deg = nz + e; + zero_poly[cur_deg + 1u] = zero_poly[cur_deg + 1u] + zero_poly[cur_deg]; + for (consteig::Size i = cur_deg; i > 0u; --i) + zero_poly[i] = zero_poly[i] + zero_poly[i - 1u]; + } + const consteig::Size num_deg = nz + n_extra; + + // Find matching frequency w_c (avoid collision with poles/zeros). + const T tol = gcem::sqrt(consteig::epsilon()); + T w_c = static_cast(0); + for (consteig::Size attempt = 0; attempt < 1000u; ++attempt) + { + bool collision = false; + for (consteig::Size i = 0; i < N && !collision; ++i) + { + const T dr = w_c - ftf.poles[i].real; + const T di = ftf.poles[i].imag; + if (dr * dr + di * di < tol * tol) + collision = true; + } + for (consteig::Size i = 0; i < nz && !collision; ++i) + { + const T dr = w_c - ftf.zeros[i].real; + const T di = ftf.zeros[i].imag; + if (dr * dr + di * di < tol * tol) + collision = true; + } + if (!collision) + break; + w_c += static_cast(0.1) / Ts; + } + + // Compute discrete gain k_d matching H_d(w_d) = H_c(w_c). + const Complex w_c_cx{w_c, static_cast(0)}; + const Complex w_d_cx = + consteig::exp(w_c_cx * Complex{Ts, static_cast(0)}); + + Complex num_c_cx{static_cast(1), static_cast(0)}; + for (consteig::Size i = 0; i < nz; ++i) + num_c_cx = num_c_cx * (w_c_cx - ftf.zeros[i]); + + Complex den_c_cx{static_cast(1), static_cast(0)}; + for (consteig::Size i = 0; i < N; ++i) + den_c_cx = den_c_cx * (w_c_cx - ftf.poles[i]); + + Complex num_d_cx{static_cast(1), static_cast(0)}; + for (consteig::Size i = 0; i < N; ++i) + num_d_cx = num_d_cx * (w_d_cx - p_d_vals[i]); + + Complex den_d_cx{static_cast(1), static_cast(0)}; + for (consteig::Size i = 0; i < nz; ++i) + den_d_cx = den_d_cx * (w_d_cx - z_d_finite[i]); + for (consteig::Size e = 0; e < n_extra; ++e) + den_d_cx = den_d_cx * + (w_d_cx - Complex{static_cast(-1), static_cast(0)}); + + const Complex gain_num = + Complex{k_c, static_cast(0)} * num_c_cx * num_d_cx; + const Complex gain_den = den_c_cx * den_d_cx; + const T gain_den_sq = + gain_den.real * gain_den.real + gain_den.imag * gain_den.imag; + const T k_d = + (gain_num.real * gain_den.real + gain_num.imag * gain_den.imag) / + gain_den_sq; + + // Assemble output TF. + TransferFunction tf{}; + for (consteig::Size i = 0; i <= N; ++i) + tf.a[i] = pole_poly[i].real; + const consteig::Size pad = N - num_deg; + for (consteig::Size i = 0; i <= num_deg; ++i) + tf.b[pad + i] = k_d * zero_poly[i].real; + + return tf; +} + // Tustin (bilinear) discretization // // Parameterized by alpha = 2/Ts (standard) or wc/tan(wc*Ts/2) (prewarped). @@ -699,6 +912,35 @@ constexpr TransferFunction analog_to_digital( return ss_to_tf(tustin_discretize(sys_c, Ts, tag)); } +// discretize_with_factored: tag-dispatched TF discretization using FactoredTF. +// ZOH and MatchedZ use the factored path; Tustin falls back to the polynomial +// path. + +template +constexpr TransferFunction discretize_with_factored( + const T (&b_c)[N + 1u], const T (&a_c)[N + 1u], const FactoredTF &ftf, + T Ts, ZOH) +{ + return ss_to_tf( + zoh_discretize_with_evals(tf_to_ss(b_c, a_c), Ts, ftf)); +} + +template +constexpr TransferFunction discretize_with_factored( + const T (&b_c)[N + 1u], const T (&a_c)[N + 1u], const FactoredTF &ftf, + T Ts, MatchedZ) +{ + return matched_z_discretize_factored(ftf, Ts); +} + +template +constexpr TransferFunction discretize_with_factored( + const T (&b_c)[N + 1u], const T (&a_c)[N + 1u], + const FactoredTF & /*ftf*/, T Ts, M method_tag) +{ + return analog_to_digital(b_c, a_c, Ts, method_tag); +} + // analog_to_digital (TF overloads) template diff --git a/include/constfilt/elliptic.hpp b/include/constfilt/elliptic.hpp index 7b6f03e..434e8dd 100644 --- a/include/constfilt/elliptic.hpp +++ b/include/constfilt/elliptic.hpp @@ -43,6 +43,8 @@ class Elliptic T sample_rate_hz) : AnalogFilter( compute_continuous_tf(cutoff_hz, ripple_db, attenuation_db), + compute_factored_tf(cutoff_hz, ripple_db, attenuation_db, + FilterType{}), sample_rate_hz, make_tustin_tag(cutoff_hz, BoundMethod{})) { } @@ -306,13 +308,51 @@ class Elliptic Complex{static_cast(0), static_cast(0)} - root * poly[0]; } + // Fills poles[] and zeros[] with the normalized (wc=1) prototype poles and + // zeros derived from the elliptic machinery (steps 3-4 of elliptic_tf). + // Poles: M conjugate pairs + optional real pole for odd N. + // Zeros: M conjugate pairs on the imaginary axis (+/-j*omega_z). + static constexpr void compute_prototype_poles_zeros( + T q, T sig0, T k, Complex (&poles)[N], consteig::Size &pole_cnt, + Complex (&zeros)[N], consteig::Size &zero_cnt) + { + const T ws = static_cast(1) / k; + const T sqrt_ws = gcem::sqrt(ws); + const T w = gcem::sqrt((static_cast(1) + k * sig0 * sig0) * + (static_cast(1) + sig0 * sig0 / k)); + + pole_cnt = 0u; + zero_cnt = 0u; + + for (consteig::Size ii = 1u; ii <= M; ++ii) + { + const T wi = compute_wi(ii, q); + const T Vi = gcem::sqrt((static_cast(1) - k * wi * wi) * + (static_cast(1) - wi * wi / k)); + + const T omega_z = sqrt_ws / wi; + zeros[zero_cnt++] = Complex{static_cast(0), omega_z}; + zeros[zero_cnt++] = Complex{static_cast(0), -omega_z}; + + const T denom = static_cast(1) + sig0 * sig0 * wi * wi; + const T p_re = sqrt_ws * (-sig0 * Vi) / denom; + const T p_im = sqrt_ws * (wi * w) / denom; + poles[pole_cnt++] = Complex{p_re, p_im}; + poles[pole_cnt++] = Complex{p_re, -p_im}; + } + + if (N % 2u == 1u) + poles[pole_cnt++] = Complex{-sig0 * sqrt_ws, static_cast(0)}; + } + // Low-pass elliptic transfer function (ncauer theta-function algorithm). // // Steps: // 1. Compute nome q from k1 via modular identity: q = q1^(1/N). // 2. Recover design modulus k from q via theta functions. // 3. Pole-shift sig0 via theta series. - // 4. Zero positions wi via theta series. + // 4. Zero positions wi via theta series (via + // compute_prototype_poles_zeros). // 5. Build s-domain polynomials from poles/zeros, scale by sqrt(ws). // 6. Gain normalization: H(0)=1 (odd N), H(0)=Gp (even N). // 7. Scale for passband cutoff wc. @@ -346,21 +386,13 @@ class Elliptic // k) via theta series, where l = acoth(10^(Rp/20)) / N. const T sig0 = compute_sig0(ripple_db, q); - // Step 4: zero and pole positions for each conjugate pair ii=1..M. - // ws = 1/k normalized stopband edge (>1) - // wi = sn(mu*K(k)/N, k) via theta series (zero/pole spacing in - // elliptic frequency - // space) - // Vi = cn(mu*K(k)/N, k) * dn(mu*K(k)/N, k) - // w = sqrt((1 + k*sig0^2)(1 + sig0^2/k)) - // - // zeros: +/-j * sqrt(ws) / wi (on the imaginary axis) - // poles: sqrt(ws) * (-sig0*Vi +/- j*wi*w) / (1 + sig0^2*wi^2) - // real pole (odd N only): -sig0 * sqrt(ws) - const T ws = static_cast(1) / k; - const T sqrt_ws = gcem::sqrt(ws); - const T w = gcem::sqrt((static_cast(1) + k * sig0 * sig0) * - (static_cast(1) + sig0 * sig0 / k)); + // Step 4: prototype poles and zeros. + Complex poles_proto[N]{}; + Complex zeros_proto[N]{}; + consteig::Size pole_cnt = 0u; + consteig::Size zero_cnt = 0u; + compute_prototype_poles_zeros(q, sig0, k, poles_proto, pole_cnt, + zeros_proto, zero_cnt); Complex poly_a[N + 1u]{}; Complex poly_b[N + 1u]{}; @@ -368,36 +400,12 @@ class Elliptic poly_b[0] = Complex{static_cast(1), static_cast(0)}; consteig::Size deg_a = 0u; - consteig::Size deg_b = 0u; - - for (consteig::Size ii = 1u; ii <= M; ++ii) - { - const T wi = compute_wi(ii, q); - const T Vi = gcem::sqrt((static_cast(1) - k * wi * wi) * - (static_cast(1) - wi * wi / k)); - - const T omega_z = sqrt_ws / wi; - ++deg_b; - poly_mul_root(poly_b, deg_b, Complex{static_cast(0), omega_z}); - ++deg_b; - poly_mul_root(poly_b, deg_b, Complex{static_cast(0), -omega_z}); + for (consteig::Size j = 0u; j < pole_cnt; ++j) + poly_mul_root(poly_a, ++deg_a, poles_proto[j]); - const T denom = static_cast(1) + sig0 * sig0 * wi * wi; - const T p_re = sqrt_ws * (-sig0 * Vi) / denom; - const T p_im = sqrt_ws * (wi * w) / denom; - - ++deg_a; - poly_mul_root(poly_a, deg_a, Complex{p_re, p_im}); - ++deg_a; - poly_mul_root(poly_a, deg_a, Complex{p_re, -p_im}); - } - - if (N % 2u == 1u) - { - ++deg_a; - poly_mul_root(poly_a, deg_a, - Complex{-sig0 * sqrt_ws, static_cast(0)}); - } + consteig::Size deg_b = 0u; + for (consteig::Size j = 0u; j < zero_cnt; ++j) + poly_mul_root(poly_b, ++deg_b, zeros_proto[j]); // Step 5: convert ascending -> descending; take real parts (imag ~ 0). for (consteig::Size i = 0u; i <= N; ++i) @@ -449,6 +457,105 @@ class Elliptic b[j] = b_lp[N - j] * sc; } } + + // LP FactoredTF: prototype poles/zeros scaled by wc; gain from polynomial. + static constexpr FactoredTF compute_factored_tf(T cutoff_hz, + T ripple_db, + T attenuation_db, + LowPass) + { + const T wc = static_cast(2) * static_cast(GCEM_PI) * cutoff_hz; + const T ep = gcem::sqrt(from_db10(ripple_db) - static_cast(1)); + const T es = gcem::sqrt(from_db10(attenuation_db) - static_cast(1)); + const T k1 = ep / es; + const T q1 = compute_nome(k1); + const T q = gcem::exp(gcem::log(q1) / static_cast(N)); + const T k = modulus_from_nome(q); + const T sig0 = compute_sig0(ripple_db, q); + + Complex poles_proto[N]{}; + Complex zeros_proto[N]{}; + consteig::Size pole_cnt = 0u; + consteig::Size zero_cnt = 0u; + compute_prototype_poles_zeros(q, sig0, k, poles_proto, pole_cnt, + zeros_proto, zero_cnt); + + FactoredTF ftf{}; + ftf.nz = zero_cnt; + for (consteig::Size i = 0u; i < pole_cnt; ++i) + ftf.poles[i] = + Complex{wc * poles_proto[i].real, wc * poles_proto[i].imag}; + for (consteig::Size i = 0u; i < zero_cnt; ++i) + ftf.zeros[i] = + Complex{wc * zeros_proto[i].real, wc * zeros_proto[i].imag}; + + // Gain from the polynomial TF (captures the normalization from steps + // 6-7). + T b_tmp[N + 1u]{}; + T a_tmp[N + 1u]{}; + elliptic_tf(wc, ripple_db, attenuation_db, b_tmp, a_tmp, LowPass{}); + consteig::Size d_b = 0u; + while (d_b <= N && b_tmp[d_b] == static_cast(0)) + ++d_b; + ftf.gain = (d_b > N) ? static_cast(0) : b_tmp[d_b] / a_tmp[0]; + + return ftf; + } + + // HP FactoredTF: derive from LP prototype via LP-to-HP transform (s -> + // wc/s). LP pole p_lp -> HP pole wc/p_lp; LP zero j*omega_z -> HP zero + // -j*wc/omega_z. For odd N: one extra zero at s=0 (LP strictly proper -> HP + // has zero at origin). + static constexpr FactoredTF compute_factored_tf(T cutoff_hz, + T ripple_db, + T attenuation_db, + HighPass) + { + const T wc = static_cast(2) * static_cast(GCEM_PI) * cutoff_hz; + + // Normalized LP prototype at wc=1 (cutoff_hz = 1/(2*pi)). + const T norm_cutoff = + static_cast(1) / (static_cast(2) * static_cast(GCEM_PI)); + const FactoredTF lp = compute_factored_tf( + norm_cutoff, ripple_db, attenuation_db, LowPass{}); + + FactoredTF ftf{}; + + // HP poles: wc / lp_pole (complex division) + for (consteig::Size i = 0u; i < N; ++i) + { + const Complex &p = lp.poles[i]; + const T denom_sq = p.real * p.real + p.imag * p.imag; + ftf.poles[i] = + Complex{wc * p.real / denom_sq, -wc * p.imag / denom_sq}; + } + + // HP zeros: wc / lp_zero (LP zeros are pure imaginary: {0, + // +/-omega_z}) + consteig::Size hp_nz = 0u; + for (consteig::Size i = 0u; i < lp.nz; ++i) + { + const Complex &z = lp.zeros[i]; + const T denom_sq = z.real * z.real + z.imag * z.imag; + ftf.zeros[hp_nz++] = + Complex{wc * z.real / denom_sq, -wc * z.imag / denom_sq}; + } + // For odd N: LP->HP adds a zero at s=0 (from the strictly-proper LP). + if (N % 2u == 1u) + ftf.zeros[hp_nz++] = Complex{static_cast(0), static_cast(0)}; + ftf.nz = hp_nz; + + // Gain from the HP polynomial TF. + T b_tmp[N + 1u]{}; + T a_tmp[N + 1u]{}; + elliptic_tf(wc, ripple_db, attenuation_db, b_tmp, a_tmp, HighPass{}); + consteig::Size d_b = 0u; + while (d_b <= N && b_tmp[d_b] == static_cast(0)) + ++d_b; + ftf.gain = (d_b > N) ? static_cast(0) : b_tmp[d_b] / a_tmp[0]; + + return ftf; + } }; } // namespace constfilt From d4f503f6e0a7e6769b7b974a1be73b0672d93479 Mon Sep 17 00:00:00 2001 From: Mitchell Thompkins Date: Thu, 25 Jun 2026 19:18:56 +0000 Subject: [PATCH 04/11] Use Vandermonde eigenvectors for ZOH CCF matrix exponential CCF eigenvectors are Vandermonde columns, so V is built directly from the analytic poles rather than via inverse iteration on the ill-conditioned CCF matrix. Fixes Butterworth ZOH accuracy to N=12. --- .../profiling/results/accuracy_gcc_15.2.0.csv | 82 +++++++------- .../profiling/results/accuracy_gcc_15.2.0.png | 4 +- .../results/compile_times_gcc_15.2.0.csv | 102 +++++++++--------- .../results/compile_times_gcc_15.2.0.png | 4 +- docs/profiling/results/runtime_gcc_15.2.0.csv | 78 +++++++------- docs/profiling/results/runtime_gcc_15.2.0.png | 4 +- include/constfilt/analog_filter.hpp | 5 + include/constfilt/discretize.hpp | 54 +++++----- 8 files changed, 172 insertions(+), 161 deletions(-) diff --git a/docs/profiling/results/accuracy_gcc_15.2.0.csv b/docs/profiling/results/accuracy_gcc_15.2.0.csv index dcfad8b..2da23ae 100644 --- a/docs/profiling/results/accuracy_gcc_15.2.0.csv +++ b/docs/profiling/results/accuracy_gcc_15.2.0.csv @@ -3,29 +3,29 @@ # os: Alpine Linux library,filter_type,order,method,max_b_err,max_a_err,max_step_err constfilt,butterworth,1,zoh,0.000000e+00,0.000000e+00,0.000000e+00 -constfilt,butterworth,2,zoh,3.608225e-16,2.220446e-16,1.110223e-15 -constfilt,butterworth,3,zoh,9.020562e-16,1.776357e-15,7.771561e-15 -constfilt,butterworth,4,zoh,8.971990e-15,8.881784e-16,1.587619e-14 -constfilt,butterworth,5,zoh,6.945833e-15,3.552714e-15,1.809664e-14 -constfilt,butterworth,6,zoh,3.284352e-14,1.598721e-14,6.572520e-14 -constfilt,butterworth,7,zoh,1.100200e-11,4.796163e-14,1.254094e-12 -constfilt,butterworth,8,zoh,2.955275e-08,5.062617e-13,2.134253e-09 -constfilt,butterworth,9,zoh,7.577398e-05,8.526513e-13,3.692378e-06 -constfilt,butterworth,10,zoh,2.065719e-01,2.030731e-11,6.461003e-03 -constfilt,butterworth,11,zoh,2.543305e+02,3.566485e-08,4.685698e+00 -constfilt,butterworth,12,zoh,1.799631e+03,3.915801e-10,2.955113e+02 +constfilt,butterworth,2,zoh,4.163336e-17,2.220446e-16,1.221245e-15 +constfilt,butterworth,3,zoh,9.714451e-17,1.776357e-15,8.548717e-15 +constfilt,butterworth,4,zoh,1.110223e-15,5.329071e-15,1.487699e-14 +constfilt,butterworth,5,zoh,7.826855e-16,2.220446e-15,1.565414e-14 +constfilt,butterworth,6,zoh,1.077333e-13,8.260059e-14,1.327827e-13 +constfilt,butterworth,7,zoh,4.350617e-13,2.184919e-13,4.956036e-13 +constfilt,butterworth,8,zoh,1.781413e-12,7.887024e-13,1.174533e-12 +constfilt,butterworth,9,zoh,2.683072e-11,2.824407e-12,1.341810e-11 +constfilt,butterworth,10,zoh,2.428169e-10,4.609291e-11,8.981571e-11 +constfilt,butterworth,11,zoh,1.938709e-10,6.137313e-10,8.004775e-11 +constfilt,butterworth,12,zoh,5.023309e-09,3.361052e-09,9.453389e-10 constfilt,butterworth,1,matchedz,5.551115e-17,0.000000e+00,2.220446e-16 -constfilt,butterworth,2,matchedz,2.775558e-17,2.220446e-16,1.110223e-15 -constfilt,butterworth,3,matchedz,1.387779e-17,1.776357e-15,1.998401e-15 -constfilt,butterworth,4,matchedz,3.122502e-17,1.776357e-15,9.103829e-15 -constfilt,butterworth,5,matchedz,4.163336e-17,4.884981e-15,1.709743e-14 -constfilt,butterworth,6,matchedz,2.168404e-17,4.884981e-14,1.030287e-13 -constfilt,butterworth,7,matchedz,1.170938e-17,9.414691e-14,3.456124e-13 -constfilt,butterworth,8,matchedz,6.505213e-19,2.664535e-14,9.876544e-13 -constfilt,butterworth,9,matchedz,6.179952e-18,6.927792e-13,1.509015e-12 -constfilt,butterworth,10,matchedz,1.523304e-16,5.579892e-11,1.363021e-11 -constfilt,butterworth,11,matchedz,2.168404e-19,1.847411e-13,4.028200e-11 -constfilt,butterworth,12,matchedz,9.662952e-17,3.253149e-10,1.482581e-10 +constfilt,butterworth,2,matchedz,0.000000e+00,2.220446e-16,1.332268e-15 +constfilt,butterworth,3,matchedz,1.387779e-17,4.440892e-16,1.998401e-15 +constfilt,butterworth,4,matchedz,3.816392e-17,1.776357e-15,4.440892e-15 +constfilt,butterworth,5,matchedz,1.040834e-17,3.552714e-15,1.487699e-14 +constfilt,butterworth,6,matchedz,3.469447e-18,1.776357e-15,5.384582e-14 +constfilt,butterworth,7,matchedz,2.602085e-18,7.105427e-15,3.296252e-13 +constfilt,butterworth,8,matchedz,2.168404e-19,5.329071e-15,9.885426e-13 +constfilt,butterworth,9,matchedz,7.589415e-19,4.440892e-15,1.743494e-12 +constfilt,butterworth,10,matchedz,1.626303e-19,7.105427e-15,7.760681e-12 +constfilt,butterworth,11,matchedz,5.421011e-20,1.421085e-14,4.277068e-11 +constfilt,butterworth,12,matchedz,1.490778e-19,1.847411e-13,9.285439e-11 constfilt,butterworth,1,prewarp,8.326673e-17,2.220446e-16,2.220446e-16 constfilt,butterworth,2,prewarp,1.249001e-16,0.000000e+00,2.331468e-15 constfilt,butterworth,3,prewarp,1.283695e-16,2.442491e-15,1.776357e-15 @@ -38,26 +38,26 @@ constfilt,butterworth,9,prewarp,4.832289e-16,2.842171e-14,4.035883e-12 constfilt,butterworth,10,prewarp,3.041296e-15,2.486900e-14,6.625700e-12 constfilt,butterworth,11,prewarp,1.426024e-15,3.197442e-13,1.841038e-11 constfilt,butterworth,12,prewarp,1.350920e-14,1.563194e-13,5.809664e-11 -constfilt,elliptic,2,zoh,2.395861e-13,2.842171e-14,1.409099e-13 -constfilt,elliptic,3,zoh,5.400819e-13,9.858780e-14,3.640976e-13 -constfilt,elliptic,4,zoh,5.393082e-13,1.074696e-13,8.994472e-13 -constfilt,elliptic,5,zoh,1.066647e-12,4.938272e-13,1.474598e-12 -constfilt,elliptic,6,zoh,1.760619e-12,1.481482e-12,1.819322e-12 -constfilt,elliptic,7,zoh,1.055545e-11,1.796252e-11,9.916956e-12 -constfilt,elliptic,8,zoh,2.154665e-09,2.641301e-09,7.933327e-10 -constfilt,elliptic,9,zoh,3.582365e-07,1.304602e-07,9.765907e-07 -constfilt,elliptic,10,zoh,3.461312e-03,3.213908e-06,6.243202e-03 -constfilt,elliptic,11,zoh,5.384745e-01,4.802802e-05,7.544774e-02 -constfilt,elliptic,12,zoh,5.474468e+02,5.003185e-04,2.311049e+01 +constfilt,elliptic,2,zoh,2.395861e-13,2.831069e-14,1.409099e-13 +constfilt,elliptic,3,zoh,5.401790e-13,9.658940e-14,3.639866e-13 +constfilt,elliptic,4,zoh,5.369767e-13,1.052491e-13,9.001133e-13 +constfilt,elliptic,5,zoh,1.077680e-12,4.902745e-13,1.481482e-12 +constfilt,elliptic,6,zoh,1.853295e-12,1.490363e-12,1.784017e-12 +constfilt,elliptic,7,zoh,1.118322e-11,1.780975e-11,1.006395e-11 +constfilt,elliptic,8,zoh,1.626981e-09,2.650339e-09,7.395398e-10 +constfilt,elliptic,9,zoh,7.139073e-08,1.304869e-07,1.850400e-08 +constfilt,elliptic,10,zoh,1.794838e-06,3.207371e-06,2.383539e-07 +constfilt,elliptic,11,zoh,2.479973e-05,4.805415e-05,1.939256e-06 +constfilt,elliptic,12,zoh,2.788489e-04,4.995743e-04,3.710091e-03 constfilt,elliptic,2,matchedz,3.854481e-11,2.786660e-14,1.928402e-11 -constfilt,elliptic,3,matchedz,5.143178e-13,9.592327e-14,3.533007e-13 -constfilt,elliptic,4,matchedz,1.324132e-13,1.101341e-13,8.861800e-13 -constfilt,elliptic,5,matchedz,1.020281e-12,4.991563e-13,1.489919e-12 -constfilt,elliptic,6,matchedz,6.720215e-13,1.428191e-12,1.833089e-12 -constfilt,elliptic,7,matchedz,1.035810e-11,1.789147e-11,1.008205e-11 -constfilt,elliptic,8,matchedz,6.505636e-10,2.641826e-09,7.457389e-10 -constfilt,elliptic,9,matchedz,6.787615e-08,1.304656e-07,1.862001e-08 -constfilt,elliptic,10,matchedz,7.442555e-07,3.213346e-06,2.409790e-07 +constfilt,elliptic,3,matchedz,5.143108e-13,9.614531e-14,3.533285e-13 +constfilt,elliptic,4,matchedz,1.324114e-13,1.088019e-13,8.872902e-13 +constfilt,elliptic,5,matchedz,1.020448e-12,4.813927e-13,1.496137e-12 +constfilt,elliptic,6,matchedz,6.720492e-13,1.449507e-12,1.837530e-12 +constfilt,elliptic,7,matchedz,1.035799e-11,1.789147e-11,1.007983e-11 +constfilt,elliptic,8,matchedz,6.505633e-10,2.641826e-09,7.456755e-10 +constfilt,elliptic,9,matchedz,6.787615e-08,1.304655e-07,1.861996e-08 +constfilt,elliptic,10,matchedz,7.442555e-07,3.213341e-06,2.409786e-07 constfilt,elliptic,11,matchedz,2.360767e-05,4.802825e-05,1.950004e-06 constfilt,elliptic,12,matchedz,1.156955e-04,5.002874e-04,1.090689e-05 constfilt,elliptic,2,prewarp,1.637024e-13,3.286260e-14,9.213463e-14 diff --git a/docs/profiling/results/accuracy_gcc_15.2.0.png b/docs/profiling/results/accuracy_gcc_15.2.0.png index 996c68d..55c616b 100644 --- a/docs/profiling/results/accuracy_gcc_15.2.0.png +++ b/docs/profiling/results/accuracy_gcc_15.2.0.png @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:3063f67ee8c7ca693c9bc690f2bf9e9e3a35e86b8025f1938552ae6ef2271908 -size 169290 +oid sha256:cfa95043df54bef12e6843007eb4561f5f086a51eb90afc46c72febefc5d4868 +size 168736 diff --git a/docs/profiling/results/compile_times_gcc_15.2.0.csv b/docs/profiling/results/compile_times_gcc_15.2.0.csv index d7f20e2..c924dfb 100644 --- a/docs/profiling/results/compile_times_gcc_15.2.0.csv +++ b/docs/profiling/results/compile_times_gcc_15.2.0.csv @@ -2,54 +2,54 @@ # compiler: g++ (Alpine 15.2.0) 15.2.0 # os: Alpine Linux filter_type,order,method,compile_time_sec,max_rss_kb,exit_code -butterworth,10,matchedz,0.37,73000,0 -butterworth,10,tustin,0.45,83112,0 -butterworth,10,zoh,1.07,158552,0 -butterworth,11,matchedz,0.50,88108,0 -butterworth,11,tustin,0.62,102420,0 -butterworth,11,zoh,1.42,203688,0 -butterworth,12,matchedz,0.63,105460,0 -butterworth,12,tustin,0.81,126272,0 -butterworth,12,zoh,1.93,264624,0 -butterworth,1,matchedz,0.07,36764,0 -butterworth,1,tustin,0.06,36776,0 -butterworth,1,zoh,0.07,37024,0 -butterworth,2,matchedz,0.07,37096,0 -butterworth,2,tustin,0.06,37268,0 -butterworth,2,zoh,0.07,38084,0 -butterworth,4,matchedz,0.09,39528,0 -butterworth,4,tustin,0.09,40116,0 -butterworth,4,zoh,0.14,44948,0 -butterworth,6,matchedz,0.12,44468,0 -butterworth,6,tustin,0.14,45600,0 -butterworth,6,zoh,0.27,60956,0 -butterworth,8,matchedz,0.21,54360,0 -butterworth,8,tustin,0.25,58428,0 -butterworth,8,zoh,0.54,94996,0 -butterworth,9,matchedz,0.28,62652,0 -butterworth,9,tustin,0.34,69284,0 -butterworth,9,zoh,0.77,123156,0 -elliptic,10,matchedz,0.72,112964,0 -elliptic,10,tustin,0.53,89548,0 -elliptic,10,zoh,1.13,163684,0 -elliptic,11,matchedz,0.95,140320,0 -elliptic,11,tustin,0.69,109080,0 -elliptic,11,zoh,1.52,210084,0 -elliptic,12,matchedz,1.30,181428,0 -elliptic,12,tustin,0.92,134184,0 -elliptic,12,zoh,2.03,271788,0 -elliptic,2,matchedz,0.11,40780,0 -elliptic,2,tustin,0.10,40608,0 -elliptic,2,zoh,0.11,41888,0 -elliptic,4,matchedz,0.14,45328,0 -elliptic,4,tustin,0.13,43936,0 -elliptic,4,zoh,0.18,49076,0 -elliptic,6,matchedz,0.22,53920,0 -elliptic,6,tustin,0.20,50064,0 -elliptic,6,zoh,0.31,64840,0 -elliptic,8,matchedz,0.40,74308,0 -elliptic,8,tustin,0.31,63908,0 -elliptic,8,zoh,0.61,100160,0 -elliptic,9,matchedz,0.53,89588,0 -elliptic,9,tustin,0.40,74688,0 -elliptic,9,zoh,0.84,128436,0 +butterworth,10,matchedz,0.07,38624,0 +butterworth,10,tustin,0.46,84536,0 +butterworth,10,zoh,0.55,96428,0 +butterworth,11,matchedz,0.07,38560,0 +butterworth,11,tustin,0.63,103800,0 +butterworth,11,zoh,0.74,119380,0 +butterworth,12,matchedz,0.07,38796,0 +butterworth,12,tustin,0.83,127672,0 +butterworth,12,zoh,0.96,146664,0 +butterworth,1,matchedz,0.06,37368,0 +butterworth,1,tustin,0.07,38072,0 +butterworth,1,zoh,0.07,38176,0 +butterworth,2,matchedz,0.07,37676,0 +butterworth,2,tustin,0.07,38636,0 +butterworth,2,zoh,0.08,39132,0 +butterworth,4,matchedz,0.06,37848,0 +butterworth,4,tustin,0.09,41480,0 +butterworth,4,zoh,0.11,42648,0 +butterworth,6,matchedz,0.07,38144,0 +butterworth,6,tustin,0.14,47152,0 +butterworth,6,zoh,0.16,50104,0 +butterworth,8,matchedz,0.07,38308,0 +butterworth,8,tustin,0.26,59960,0 +butterworth,8,zoh,0.30,66636,0 +butterworth,9,matchedz,0.07,38380,0 +butterworth,9,tustin,0.34,70816,0 +butterworth,9,zoh,0.42,79868,0 +elliptic,10,matchedz,0.27,56120,0 +elliptic,10,tustin,0.66,102124,0 +elliptic,10,zoh,0.77,114252,0 +elliptic,11,matchedz,0.27,54692,0 +elliptic,11,tustin,0.82,120504,0 +elliptic,11,zoh,0.94,135696,0 +elliptic,12,matchedz,0.30,57064,0 +elliptic,12,tustin,1.04,147284,0 +elliptic,12,zoh,1.23,166888,0 +elliptic,2,matchedz,0.16,46156,0 +elliptic,2,tustin,0.17,47204,0 +elliptic,2,zoh,0.17,47716,0 +elliptic,4,matchedz,0.20,49112,0 +elliptic,4,tustin,0.22,52360,0 +elliptic,4,zoh,0.23,53788,0 +elliptic,6,matchedz,0.22,50672,0 +elliptic,6,tustin,0.29,59368,0 +elliptic,6,zoh,0.32,62320,0 +elliptic,8,matchedz,0.24,52772,0 +elliptic,8,tustin,0.44,74612,0 +elliptic,8,zoh,0.48,81096,0 +elliptic,9,matchedz,0.26,54120,0 +elliptic,9,tustin,0.53,86104,0 +elliptic,9,zoh,0.60,95208,0 diff --git a/docs/profiling/results/compile_times_gcc_15.2.0.png b/docs/profiling/results/compile_times_gcc_15.2.0.png index 0a1b52f..9ccca11 100644 --- a/docs/profiling/results/compile_times_gcc_15.2.0.png +++ b/docs/profiling/results/compile_times_gcc_15.2.0.png @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:4a4f746ddc0d9c0f7f9380a8c7fe02c44288304df30b974042213c7456a1af73 -size 150965 +oid sha256:51da0826e6f341c9889e1209104a8a9affde570aa3dd9c642befc9f92a071dbc +size 138263 diff --git a/docs/profiling/results/runtime_gcc_15.2.0.csv b/docs/profiling/results/runtime_gcc_15.2.0.csv index 5b7cc99..e82b683 100644 --- a/docs/profiling/results/runtime_gcc_15.2.0.csv +++ b/docs/profiling/results/runtime_gcc_15.2.0.csv @@ -2,42 +2,42 @@ # compiler: g++ (Alpine 15.2.0) 15.2.0 # os: Alpine Linux library,filter_type,order,method,ns_per_sample,msa_per_s,dc_gain -constfilt,Butterworth,1,tustin,6.415,155.9,1.00000000 -constfilt,Butterworth,2,tustin,10.876,91.9,1.00000000 -constfilt,Butterworth,4,tustin,18.954,52.8,1.00000000 -constfilt,Butterworth,6,tustin,26.164,38.2,1.00000000 -constfilt,Butterworth,8,tustin,34.550,28.9,1.00000000 -constfilt,Butterworth,1,zoh,6.420,155.8,1.00000000 -constfilt,Butterworth,2,zoh,10.826,92.4,1.00000000 -constfilt,Butterworth,4,zoh,18.535,54.0,1.00000000 -constfilt,Butterworth,6,zoh,26.066,38.4,1.00000000 -constfilt,Butterworth,8,zoh,34.415,29.1,1.00000000 -constfilt,Butterworth,1,matchedz,6.423,155.7,1.00000000 -constfilt,Butterworth,2,matchedz,10.839,92.3,1.00000000 -constfilt,Butterworth,4,matchedz,18.858,53.0,1.00000000 -constfilt,Butterworth,6,matchedz,26.130,38.3,1.00000000 -constfilt,Butterworth,8,matchedz,34.357,29.1,1.00000000 -constfilt,Elliptic,2,tustin,10.844,92.2,0.94406088 -constfilt,Elliptic,4,tustin,18.373,54.4,0.94406088 -constfilt,Elliptic,6,tustin,26.217,38.1,0.94406088 -constfilt,Elliptic,8,tustin,34.485,29.0,0.94406088 -constfilt,Elliptic,2,zoh,10.828,92.4,0.94406088 -constfilt,Elliptic,4,zoh,18.359,54.5,0.94406088 -constfilt,Elliptic,6,zoh,26.095,38.3,0.94406088 -constfilt,Elliptic,8,zoh,34.249,29.2,0.94406088 -constfilt,Elliptic,2,matchedz,11.091,90.2,0.94406088 -constfilt,Elliptic,4,matchedz,18.904,52.9,0.94406088 -constfilt,Elliptic,6,matchedz,26.024,38.4,0.94406088 -constfilt,Elliptic,8,matchedz,34.507,29.0,0.94406088 -iir1,Butterworth,2,runtime,10.732,93.2,1.00000000 -iir1,Butterworth,4,runtime,18.155,55.1,1.00000000 -iir1,Butterworth,6,runtime,30.092,33.2,1.00000000 -iir1,Butterworth,8,runtime,40.322,24.8,1.00000000 -kfr,Butterworth,2,runtime+simd,607.261,1.6,1.00000000 -kfr,Butterworth,4,runtime+simd,402.769,2.5,1.00000000 -kfr,Butterworth,6,runtime+simd,307.412,3.3,1.00000000 -kfr,Butterworth,8,runtime+simd,307.261,3.3,1.00000000 -kfr,Elliptic,2,runtime+simd,605.280,1.7,0.94406088 -kfr,Elliptic,4,runtime+simd,401.691,2.5,0.94406088 -kfr,Elliptic,6,runtime+simd,306.657,3.3,0.94406088 -kfr,Elliptic,8,runtime+simd,308.208,3.2,0.94406088 +constfilt,Butterworth,1,tustin,6.411,156.0,1.00000000 +constfilt,Butterworth,2,tustin,10.844,92.2,1.00000000 +constfilt,Butterworth,4,tustin,18.902,52.9,1.00000000 +constfilt,Butterworth,6,tustin,25.953,38.5,1.00000000 +constfilt,Butterworth,8,tustin,34.284,29.2,1.00000000 +constfilt,Butterworth,1,zoh,6.412,156.0,1.00000000 +constfilt,Butterworth,2,zoh,10.842,92.2,1.00000000 +constfilt,Butterworth,4,zoh,18.534,54.0,1.00000000 +constfilt,Butterworth,6,zoh,25.968,38.5,1.00000000 +constfilt,Butterworth,8,zoh,34.259,29.2,1.00000000 +constfilt,Butterworth,1,matchedz,6.419,155.8,1.00000000 +constfilt,Butterworth,2,matchedz,10.846,92.2,1.00000000 +constfilt,Butterworth,4,matchedz,18.803,53.2,1.00000000 +constfilt,Butterworth,6,matchedz,26.181,38.2,1.00000000 +constfilt,Butterworth,8,matchedz,34.334,29.1,1.00000000 +constfilt,Elliptic,2,tustin,10.835,92.3,0.94406088 +constfilt,Elliptic,4,tustin,18.386,54.4,0.94406088 +constfilt,Elliptic,6,tustin,25.975,38.5,0.94406088 +constfilt,Elliptic,8,tustin,34.344,29.1,0.94406088 +constfilt,Elliptic,2,zoh,10.848,92.2,0.94406088 +constfilt,Elliptic,4,zoh,18.370,54.4,0.94406088 +constfilt,Elliptic,6,zoh,25.969,38.5,0.94406088 +constfilt,Elliptic,8,zoh,34.320,29.1,0.94406088 +constfilt,Elliptic,2,matchedz,11.065,90.4,0.94406088 +constfilt,Elliptic,4,matchedz,18.840,53.1,0.94406088 +constfilt,Elliptic,6,matchedz,26.006,38.5,0.94406088 +constfilt,Elliptic,8,matchedz,34.280,29.2,0.94406088 +iir1,Butterworth,2,runtime,10.633,94.1,1.00000000 +iir1,Butterworth,4,runtime,18.076,55.3,1.00000000 +iir1,Butterworth,6,runtime,28.678,34.9,1.00000000 +iir1,Butterworth,8,runtime,39.906,25.1,1.00000000 +kfr,Butterworth,2,runtime+simd,581.479,1.7,1.00000000 +kfr,Butterworth,4,runtime+simd,392.699,2.5,1.00000000 +kfr,Butterworth,6,runtime+simd,308.324,3.2,1.00000000 +kfr,Butterworth,8,runtime+simd,308.872,3.2,1.00000000 +kfr,Elliptic,2,runtime+simd,581.987,1.7,0.94406088 +kfr,Elliptic,4,runtime+simd,390.947,2.6,0.94406088 +kfr,Elliptic,6,runtime+simd,308.693,3.2,0.94406088 +kfr,Elliptic,8,runtime+simd,307.146,3.3,0.94406088 diff --git a/docs/profiling/results/runtime_gcc_15.2.0.png b/docs/profiling/results/runtime_gcc_15.2.0.png index 931a63a..4a90f62 100644 --- a/docs/profiling/results/runtime_gcc_15.2.0.png +++ b/docs/profiling/results/runtime_gcc_15.2.0.png @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:e23c6d526144320c1c2baf0928589539fa06e92946d213ac38e5050900538506 -size 113395 +oid sha256:897a23f8a1483607f16882fd33041d0585b73c88ff2993661cffff2a07ec50ee +size 111057 diff --git a/include/constfilt/analog_filter.hpp b/include/constfilt/analog_filter.hpp index 204895e..c4b7452 100644 --- a/include/constfilt/analog_filter.hpp +++ b/include/constfilt/analog_filter.hpp @@ -37,6 +37,11 @@ namespace constfilt // AnalogFilter(continuous_tf, sample_rate_hz, method_tag) // method_tag - method instance; for TustinPW supply // constfilt::prewarp(warp_hz) +// +// AnalogFilter(continuous_tf, ftf, sample_rate_hz, method_tag) +// ftf - factored (pole/zero/gain) form of continuous_tf; +// used by Butterworth and Elliptic to bypass roundtrip +// eigendecompositions in ZOH and MatchedZ discretization template class AnalogFilter : public Filter diff --git a/include/constfilt/discretize.hpp b/include/constfilt/discretize.hpp index 31a59df..9cbf1ad 100644 --- a/include/constfilt/discretize.hpp +++ b/include/constfilt/discretize.hpp @@ -178,20 +178,35 @@ constexpr consteig::Matrix matrix_exp( return result; } -// matrix_exp variant that accepts precomputed eigenvalues of A. -// Identical to matrix_exp except consteig::eigenvalues(A) is skipped; -// the caller supplies the eigenvalues directly via `evals`. -// For ZOH use: pass (Ts*Ac, scaled_evals) where scaled_evals(i,0) = Ts*pole_i. +// matrix_exp for controllable-canonical-form (CCF) matrices given analytic poles. +// +// The CCF matrix A_c has eigenvectors v_i = [1, p_i, p_i^2, ..., p_i^(N-1)]^T +// where p_i are the analog poles (eigenvalues of A_c). This Vandermonde +// structure lets us build V analytically from the poles, bypassing +// consteig::eigenvectors (inverse iteration on the ill-conditioned CCF matrix). +// +// exp(Ts*Ac) = V * diag(exp(Ts*p_i)) * V^{-1} +// V[r][i] = p_i^r (Vandermonde, built from UNSCALED poles) +// exp factor = exp(Ts * p_i) (spectral scaling) template -constexpr consteig::Matrix matrix_exp_with_evals( - const consteig::Matrix &A, - const consteig::Matrix, N, 1> &evals) +constexpr consteig::Matrix matrix_exp_ccf( + const FactoredTF &ftf, T Ts) { using Complex = consteig::Complex; using ComplexMat_NN = consteig::Matrix; using ComplexMat_N1 = consteig::Matrix; - const auto V = consteig::eigenvectors(A, evals); + // Build Vandermonde V: V[r][i] = poles[i]^r (unscaled analog poles) + ComplexMat_NN V{}; + for (consteig::Size i = 0; i < N; ++i) + { + Complex power{static_cast(1), static_cast(0)}; + for (consteig::Size r = 0; r < N; ++r) + { + V(r, i) = power; + power = power * ftf.poles[i]; + } + } const auto lu_V = consteig::lu(V); ComplexMat_NN V_inv{}; @@ -207,7 +222,8 @@ constexpr consteig::Matrix matrix_exp_with_evals( ComplexMat_NN result_c{}; for (consteig::Size i = 0; i < N; ++i) { - Complex exp_lambda = consteig::exp(evals(i, 0)); + const Complex exp_lambda = consteig::exp( + Complex{Ts * ftf.poles[i].real, Ts * ftf.poles[i].imag}); for (consteig::Size r = 0; r < N; ++r) for (consteig::Size c = 0; c < N; ++c) result_c(r, c) = @@ -258,23 +274,13 @@ constexpr StateSpace zoh_discretize(const StateSpace &sys_c, T Ts, } // ZOH discretization using analytically known poles from a FactoredTF. -// Avoids the QR eigenvalue search inside matrix_exp by supplying -// the scaled analog poles (eigenvalues of Ts*Ac) directly. +// Bypasses both the QR eigenvalue search and eigenvector inverse iteration +// by exploiting the Vandermonde structure of the CCF eigenvectors. template constexpr StateSpace zoh_discretize_with_evals( const StateSpace &sys_c, T Ts, const FactoredTF &ftf) { - using Complex = consteig::Complex; - - const consteig::Matrix TsAc = Ts * sys_c.A; - - consteig::Matrix scaled_evals{}; - for (consteig::Size i = 0; i < N; ++i) - scaled_evals(i, 0) = - Complex{Ts * ftf.poles[i].real, Ts * ftf.poles[i].imag}; - - const consteig::Matrix Ad = - matrix_exp_with_evals(TsAc, scaled_evals); + const consteig::Matrix Ad = matrix_exp_ccf(ftf, Ts); const consteig::Matrix AdmI = Ad - consteig::eye(); const consteig::Matrix rhs = AdmI * sys_c.B; @@ -927,8 +933,8 @@ constexpr TransferFunction discretize_with_factored( template constexpr TransferFunction discretize_with_factored( - const T (&b_c)[N + 1u], const T (&a_c)[N + 1u], const FactoredTF &ftf, - T Ts, MatchedZ) + const T (& /*b_c*/)[N + 1u], const T (& /*a_c*/)[N + 1u], + const FactoredTF &ftf, T Ts, MatchedZ) { return matched_z_discretize_factored(ftf, Ts); } From 8c118685854aeb313003a9aeaca5c417e7ac7ff2 Mon Sep 17 00:00:00 2001 From: Mitchell Thompkins Date: Thu, 25 Jun 2026 19:31:39 +0000 Subject: [PATCH 05/11] restore math for elliptic filter, no need to have changed that. Add braces --- include/constfilt/discretize.hpp | 43 +++++++++++++++++-- include/constfilt/elliptic.hpp | 71 +++++++++++++++++++++++++------- 2 files changed, 97 insertions(+), 17 deletions(-) diff --git a/include/constfilt/discretize.hpp b/include/constfilt/discretize.hpp index 9cbf1ad..2cd18ec 100644 --- a/include/constfilt/discretize.hpp +++ b/include/constfilt/discretize.hpp @@ -178,7 +178,8 @@ constexpr consteig::Matrix matrix_exp( return result; } -// matrix_exp for controllable-canonical-form (CCF) matrices given analytic poles. +// matrix_exp for controllable-canonical-form (CCF) matrices given analytic +// poles. // // The CCF matrix A_c has eigenvectors v_i = [1, p_i, p_i^2, ..., p_i^(N-1)]^T // where p_i are the analog poles (eigenvalues of A_c). This Vandermonde @@ -189,8 +190,8 @@ constexpr consteig::Matrix matrix_exp( // V[r][i] = p_i^r (Vandermonde, built from UNSCALED poles) // exp factor = exp(Ts * p_i) (spectral scaling) template -constexpr consteig::Matrix matrix_exp_ccf( - const FactoredTF &ftf, T Ts) +constexpr consteig::Matrix matrix_exp_ccf(const FactoredTF &ftf, + T Ts) { using Complex = consteig::Complex; using ComplexMat_NN = consteig::Matrix; @@ -216,7 +217,9 @@ constexpr consteig::Matrix matrix_exp_ccf( e_col(col, 0) = Complex{static_cast(1), static_cast(0)}; auto col_vec = consteig::lu_solve(lu_V, e_col); for (consteig::Size row = 0; row < N; ++row) + { V_inv(row, col) = col_vec(row, 0); + } } ComplexMat_NN result_c{}; @@ -225,15 +228,23 @@ constexpr consteig::Matrix matrix_exp_ccf( const Complex exp_lambda = consteig::exp( Complex{Ts * ftf.poles[i].real, Ts * ftf.poles[i].imag}); for (consteig::Size r = 0; r < N; ++r) + { for (consteig::Size c = 0; c < N; ++c) + { result_c(r, c) = result_c(r, c) + exp_lambda * V(r, i) * V_inv(i, c); + } + } } consteig::Matrix result{}; for (consteig::Size r = 0; r < N; ++r) + { for (consteig::Size c = 0; c < N; ++c) + { result(r, c) = result_c(r, c).real; + } + } return result; } @@ -696,7 +707,9 @@ constexpr TransferFunction matched_z_discretize_factored( pole_poly[k + 1u] = Complex{static_cast(0), static_cast(0)} - zk * pole_poly[k]; for (consteig::Size i = k; i > 0u; --i) + { pole_poly[i] = pole_poly[i] - zk * pole_poly[i - 1u]; + } } // Map finite zeros to z-domain; build monic numerator polynomial. @@ -712,7 +725,9 @@ constexpr TransferFunction matched_z_discretize_factored( zero_poly[k + 1u] = Complex{static_cast(0), static_cast(0)} - zk * zero_poly[k]; for (consteig::Size i = k; i > 0u; --i) + { zero_poly[i] = zero_poly[i] - zk * zero_poly[i - 1u]; + } } // Pad with zeros at z = -1 to reach numerator degree N-1. @@ -722,7 +737,9 @@ constexpr TransferFunction matched_z_discretize_factored( const consteig::Size cur_deg = nz + e; zero_poly[cur_deg + 1u] = zero_poly[cur_deg + 1u] + zero_poly[cur_deg]; for (consteig::Size i = cur_deg; i > 0u; --i) + { zero_poly[i] = zero_poly[i] + zero_poly[i - 1u]; + } } const consteig::Size num_deg = nz + n_extra; @@ -737,17 +754,23 @@ constexpr TransferFunction matched_z_discretize_factored( const T dr = w_c - ftf.poles[i].real; const T di = ftf.poles[i].imag; if (dr * dr + di * di < tol * tol) + { collision = true; + } } for (consteig::Size i = 0; i < nz && !collision; ++i) { const T dr = w_c - ftf.zeros[i].real; const T di = ftf.zeros[i].imag; if (dr * dr + di * di < tol * tol) + { collision = true; + } } if (!collision) + { break; + } w_c += static_cast(0.1) / Ts; } @@ -758,22 +781,32 @@ constexpr TransferFunction matched_z_discretize_factored( Complex num_c_cx{static_cast(1), static_cast(0)}; for (consteig::Size i = 0; i < nz; ++i) + { num_c_cx = num_c_cx * (w_c_cx - ftf.zeros[i]); + } Complex den_c_cx{static_cast(1), static_cast(0)}; for (consteig::Size i = 0; i < N; ++i) + { den_c_cx = den_c_cx * (w_c_cx - ftf.poles[i]); + } Complex num_d_cx{static_cast(1), static_cast(0)}; for (consteig::Size i = 0; i < N; ++i) + { num_d_cx = num_d_cx * (w_d_cx - p_d_vals[i]); + } Complex den_d_cx{static_cast(1), static_cast(0)}; for (consteig::Size i = 0; i < nz; ++i) + { den_d_cx = den_d_cx * (w_d_cx - z_d_finite[i]); + } for (consteig::Size e = 0; e < n_extra; ++e) + { den_d_cx = den_d_cx * (w_d_cx - Complex{static_cast(-1), static_cast(0)}); + } const Complex gain_num = Complex{k_c, static_cast(0)} * num_c_cx * num_d_cx; @@ -787,10 +820,14 @@ constexpr TransferFunction matched_z_discretize_factored( // Assemble output TF. TransferFunction tf{}; for (consteig::Size i = 0; i <= N; ++i) + { tf.a[i] = pole_poly[i].real; + } const consteig::Size pad = N - num_deg; for (consteig::Size i = 0; i <= num_deg; ++i) + { tf.b[pad + i] = k_d * zero_poly[i].real; + } return tf; } diff --git a/include/constfilt/elliptic.hpp b/include/constfilt/elliptic.hpp index 434e8dd..0568365 100644 --- a/include/constfilt/elliptic.hpp +++ b/include/constfilt/elliptic.hpp @@ -342,7 +342,9 @@ class Elliptic } if (N % 2u == 1u) + { poles[pole_cnt++] = Complex{-sig0 * sqrt_ws, static_cast(0)}; + } } // Low-pass elliptic transfer function (ncauer theta-function algorithm). @@ -351,8 +353,7 @@ class Elliptic // 1. Compute nome q from k1 via modular identity: q = q1^(1/N). // 2. Recover design modulus k from q via theta functions. // 3. Pole-shift sig0 via theta series. - // 4. Zero positions wi via theta series (via - // compute_prototype_poles_zeros). + // 4. Zero positions wi via theta series. // 5. Build s-domain polynomials from poles/zeros, scale by sqrt(ws). // 6. Gain normalization: H(0)=1 (odd N), H(0)=Gp (even N). // 7. Scale for passband cutoff wc. @@ -386,13 +387,21 @@ class Elliptic // k) via theta series, where l = acoth(10^(Rp/20)) / N. const T sig0 = compute_sig0(ripple_db, q); - // Step 4: prototype poles and zeros. - Complex poles_proto[N]{}; - Complex zeros_proto[N]{}; - consteig::Size pole_cnt = 0u; - consteig::Size zero_cnt = 0u; - compute_prototype_poles_zeros(q, sig0, k, poles_proto, pole_cnt, - zeros_proto, zero_cnt); + // Step 4: zero and pole positions for each conjugate pair ii=1..M. + // ws = 1/k normalized stopband edge (>1) + // wi = sn(mu*K(k)/N, k) via theta series (zero/pole spacing in + // elliptic frequency + // space) + // Vi = cn(mu*K(k)/N, k) * dn(mu*K(k)/N, k) + // w = sqrt((1 + k*sig0^2)(1 + sig0^2/k)) + // + // zeros: +/-j * sqrt(ws) / wi (on the imaginary axis) + // poles: sqrt(ws) * (-sig0*Vi +/- j*wi*w) / (1 + sig0^2*wi^2) + // real pole (odd N only): -sig0 * sqrt(ws) + const T ws = static_cast(1) / k; + const T sqrt_ws = gcem::sqrt(ws); + const T w = gcem::sqrt((static_cast(1) + k * sig0 * sig0) * + (static_cast(1) + sig0 * sig0 / k)); Complex poly_a[N + 1u]{}; Complex poly_b[N + 1u]{}; @@ -400,12 +409,36 @@ class Elliptic poly_b[0] = Complex{static_cast(1), static_cast(0)}; consteig::Size deg_a = 0u; - for (consteig::Size j = 0u; j < pole_cnt; ++j) - poly_mul_root(poly_a, ++deg_a, poles_proto[j]); - consteig::Size deg_b = 0u; - for (consteig::Size j = 0u; j < zero_cnt; ++j) - poly_mul_root(poly_b, ++deg_b, zeros_proto[j]); + + for (consteig::Size ii = 1u; ii <= M; ++ii) + { + const T wi = compute_wi(ii, q); + const T Vi = gcem::sqrt((static_cast(1) - k * wi * wi) * + (static_cast(1) - wi * wi / k)); + + const T omega_z = sqrt_ws / wi; + ++deg_b; + poly_mul_root(poly_b, deg_b, Complex{static_cast(0), omega_z}); + ++deg_b; + poly_mul_root(poly_b, deg_b, Complex{static_cast(0), -omega_z}); + + const T denom = static_cast(1) + sig0 * sig0 * wi * wi; + const T p_re = sqrt_ws * (-sig0 * Vi) / denom; + const T p_im = sqrt_ws * (wi * w) / denom; + + ++deg_a; + poly_mul_root(poly_a, deg_a, Complex{p_re, p_im}); + ++deg_a; + poly_mul_root(poly_a, deg_a, Complex{p_re, -p_im}); + } + + if (N % 2u == 1u) + { + ++deg_a; + poly_mul_root(poly_a, deg_a, + Complex{-sig0 * sqrt_ws, static_cast(0)}); + } // Step 5: convert ascending -> descending; take real parts (imag ~ 0). for (consteig::Size i = 0u; i <= N; ++i) @@ -483,11 +516,15 @@ class Elliptic FactoredTF ftf{}; ftf.nz = zero_cnt; for (consteig::Size i = 0u; i < pole_cnt; ++i) + { ftf.poles[i] = Complex{wc * poles_proto[i].real, wc * poles_proto[i].imag}; + } for (consteig::Size i = 0u; i < zero_cnt; ++i) + { ftf.zeros[i] = Complex{wc * zeros_proto[i].real, wc * zeros_proto[i].imag}; + } // Gain from the polynomial TF (captures the normalization from steps // 6-7). @@ -496,7 +533,9 @@ class Elliptic elliptic_tf(wc, ripple_db, attenuation_db, b_tmp, a_tmp, LowPass{}); consteig::Size d_b = 0u; while (d_b <= N && b_tmp[d_b] == static_cast(0)) + { ++d_b; + } ftf.gain = (d_b > N) ? static_cast(0) : b_tmp[d_b] / a_tmp[0]; return ftf; @@ -542,7 +581,9 @@ class Elliptic } // For odd N: LP->HP adds a zero at s=0 (from the strictly-proper LP). if (N % 2u == 1u) + { ftf.zeros[hp_nz++] = Complex{static_cast(0), static_cast(0)}; + } ftf.nz = hp_nz; // Gain from the HP polynomial TF. @@ -551,7 +592,9 @@ class Elliptic elliptic_tf(wc, ripple_db, attenuation_db, b_tmp, a_tmp, HighPass{}); consteig::Size d_b = 0u; while (d_b <= N && b_tmp[d_b] == static_cast(0)) + { ++d_b; + } ftf.gain = (d_b > N) ? static_cast(0) : b_tmp[d_b] / a_tmp[0]; return ftf; From 186746defc56399df73c40f238191bec54dcf0f6 Mon Sep 17 00:00:00 2001 From: Mitchell Thompkins Date: Thu, 25 Jun 2026 21:03:51 +0000 Subject: [PATCH 06/11] use normal var names --- include/constfilt/analog_filter.hpp | 17 ++++----- include/constfilt/butterworth.hpp | 25 +++++++------ include/constfilt/discretize.hpp | 56 ++++++++++++++--------------- include/constfilt/elliptic.hpp | 29 ++++++++------- 4 files changed, 67 insertions(+), 60 deletions(-) diff --git a/include/constfilt/analog_filter.hpp b/include/constfilt/analog_filter.hpp index c4b7452..9ce715b 100644 --- a/include/constfilt/analog_filter.hpp +++ b/include/constfilt/analog_filter.hpp @@ -38,8 +38,8 @@ namespace constfilt // method_tag - method instance; for TustinPW supply // constfilt::prewarp(warp_hz) // -// AnalogFilter(continuous_tf, ftf, sample_rate_hz, method_tag) -// ftf - factored (pole/zero/gain) form of continuous_tf; +// AnalogFilter(continuous_tf, factored_tf, sample_rate_hz, method_tag) +// factored_tf - factored (pole/zero/gain) form of continuous_tf; // used by Butterworth and Elliptic to bypass roundtrip // eigendecompositions in ZOH and MatchedZ discretization template } constexpr AnalogFilter(TransferFunction continuous_tf, - const FactoredTF &ftf, T sample_rate_hz, - BoundMethod method_tag) + const FactoredTF &factored_tf, + T sample_rate_hz, BoundMethod method_tag) : AnalogFilter(checked_discretize_factored(continuous_tf.b, - continuous_tf.a, ftf, + continuous_tf.a, factored_tf, sample_rate_hz, method_tag)) { } @@ -137,11 +137,12 @@ class AnalogFilter : public Filter static constexpr TransferFunction checked_discretize_factored(const T (&b_c)[N + 1u], const T (&a_c)[N + 1u], - const FactoredTF &ftf, T sample_rate_hz, - BoundMethod method_tag) + const FactoredTF &factored_tf, + T sample_rate_hz, BoundMethod method_tag) { return discretize_with_factored( - b_c, a_c, ftf, static_cast(1) / sample_rate_hz, method_tag); + b_c, a_c, factored_tf, static_cast(1) / sample_rate_hz, + method_tag); } }; diff --git a/include/constfilt/butterworth.hpp b/include/constfilt/butterworth.hpp index 88304c4..d150a01 100644 --- a/include/constfilt/butterworth.hpp +++ b/include/constfilt/butterworth.hpp @@ -96,17 +96,18 @@ class Butterworth static constexpr FactoredTF compute_factored_tf(T cutoff_hz, LowPass) { const T wc = static_cast(2) * static_cast(GCEM_PI) * cutoff_hz; - FactoredTF ftf{}; - ftf.nz = 0; - ftf.gain = gcem::pow(wc, static_cast(N)); + FactoredTF factored_tf{}; + factored_tf.nz = 0; + factored_tf.gain = gcem::pow(wc, static_cast(N)); for (consteig::Size k = 1u; k <= N; ++k) { const T theta = static_cast(GCEM_PI) * static_cast(2u * k + N - 1u) / static_cast(2u * N); - ftf.poles[k - 1u] = {wc * gcem::cos(theta), wc * gcem::sin(theta)}; + factored_tf.poles[k - 1u] = {wc * gcem::cos(theta), + wc * gcem::sin(theta)}; } - return ftf; + return factored_tf; } // HP: poles at wc*exp(-j*theta_k) (magnitude wc), N zeros at s=0, gain=1. @@ -114,18 +115,20 @@ class Butterworth { using Complex = consteig::Complex; const T wc = static_cast(2) * static_cast(GCEM_PI) * cutoff_hz; - FactoredTF ftf{}; - ftf.nz = N; - ftf.gain = static_cast(1); + FactoredTF factored_tf{}; + factored_tf.nz = N; + factored_tf.gain = static_cast(1); for (consteig::Size k = 1u; k <= N; ++k) { const T theta = static_cast(GCEM_PI) * static_cast(2u * k + N - 1u) / static_cast(2u * N); - ftf.poles[k - 1u] = {wc * gcem::cos(theta), -wc * gcem::sin(theta)}; - ftf.zeros[k - 1u] = Complex{static_cast(0), static_cast(0)}; + factored_tf.poles[k - 1u] = {wc * gcem::cos(theta), + -wc * gcem::sin(theta)}; + factored_tf.zeros[k - 1u] = + Complex{static_cast(0), static_cast(0)}; } - return ftf; + return factored_tf; } // Normalized Butterworth denominator coefficients (wc=1, monic). diff --git a/include/constfilt/discretize.hpp b/include/constfilt/discretize.hpp index 2cd18ec..51661f7 100644 --- a/include/constfilt/discretize.hpp +++ b/include/constfilt/discretize.hpp @@ -190,8 +190,8 @@ constexpr consteig::Matrix matrix_exp( // V[r][i] = p_i^r (Vandermonde, built from UNSCALED poles) // exp factor = exp(Ts * p_i) (spectral scaling) template -constexpr consteig::Matrix matrix_exp_ccf(const FactoredTF &ftf, - T Ts) +constexpr consteig::Matrix matrix_exp_ccf( + const FactoredTF &factored_tf, T Ts) { using Complex = consteig::Complex; using ComplexMat_NN = consteig::Matrix; @@ -205,7 +205,7 @@ constexpr consteig::Matrix matrix_exp_ccf(const FactoredTF &ftf, for (consteig::Size r = 0; r < N; ++r) { V(r, i) = power; - power = power * ftf.poles[i]; + power = power * factored_tf.poles[i]; } } @@ -225,8 +225,8 @@ constexpr consteig::Matrix matrix_exp_ccf(const FactoredTF &ftf, ComplexMat_NN result_c{}; for (consteig::Size i = 0; i < N; ++i) { - const Complex exp_lambda = consteig::exp( - Complex{Ts * ftf.poles[i].real, Ts * ftf.poles[i].imag}); + const Complex exp_lambda = consteig::exp(Complex{ + Ts * factored_tf.poles[i].real, Ts * factored_tf.poles[i].imag}); for (consteig::Size r = 0; r < N; ++r) { for (consteig::Size c = 0; c < N; ++c) @@ -289,9 +289,9 @@ constexpr StateSpace zoh_discretize(const StateSpace &sys_c, T Ts, // by exploiting the Vandermonde structure of the CCF eigenvectors. template constexpr StateSpace zoh_discretize_with_evals( - const StateSpace &sys_c, T Ts, const FactoredTF &ftf) + const StateSpace &sys_c, T Ts, const FactoredTF &factored_tf) { - const consteig::Matrix Ad = matrix_exp_ccf(ftf, Ts); + const consteig::Matrix Ad = matrix_exp_ccf(factored_tf, Ts); const consteig::Matrix AdmI = Ad - consteig::eye(); const consteig::Matrix rhs = AdmI * sys_c.B; @@ -684,16 +684,16 @@ constexpr TransferFunction matched_z_discretize_tf( // Matched-Z discretization from a FactoredTF (precomputed poles/zeros). // Bypasses both companion-matrix eigendecompositions in -// matched_z_discretize_tf. Uses ftf.poles, ftf.zeros[0..ftf.nz-1], and ftf.gain -// directly. +// matched_z_discretize_tf. Uses factored_tf.poles, +// factored_tf.zeros[0..factored_tf.nz-1], and factored_tf.gain directly. template constexpr TransferFunction matched_z_discretize_factored( - const FactoredTF &ftf, T Ts) + const FactoredTF &factored_tf, T Ts) { using Complex = consteig::Complex; - const consteig::Size nz = ftf.nz; - const T k_c = ftf.gain; + const consteig::Size nz = factored_tf.nz; + const T k_c = factored_tf.gain; // Map poles to z-domain; build monic denominator polynomial. Complex p_d_vals[N]{}; @@ -701,8 +701,8 @@ constexpr TransferFunction matched_z_discretize_factored( pole_poly[0] = Complex{static_cast(1), static_cast(0)}; for (consteig::Size k = 0; k < N; ++k) { - p_d_vals[k] = - consteig::exp(ftf.poles[k] * Complex{Ts, static_cast(0)}); + p_d_vals[k] = consteig::exp(factored_tf.poles[k] * + Complex{Ts, static_cast(0)}); const Complex zk = p_d_vals[k]; pole_poly[k + 1u] = Complex{static_cast(0), static_cast(0)} - zk * pole_poly[k]; @@ -719,8 +719,8 @@ constexpr TransferFunction matched_z_discretize_factored( zero_poly[0] = Complex{static_cast(1), static_cast(0)}; for (consteig::Size k = 0; k < nz; ++k) { - z_d_finite[k] = - consteig::exp(ftf.zeros[k] * Complex{Ts, static_cast(0)}); + z_d_finite[k] = consteig::exp(factored_tf.zeros[k] * + Complex{Ts, static_cast(0)}); const Complex zk = z_d_finite[k]; zero_poly[k + 1u] = Complex{static_cast(0), static_cast(0)} - zk * zero_poly[k]; @@ -751,8 +751,8 @@ constexpr TransferFunction matched_z_discretize_factored( bool collision = false; for (consteig::Size i = 0; i < N && !collision; ++i) { - const T dr = w_c - ftf.poles[i].real; - const T di = ftf.poles[i].imag; + const T dr = w_c - factored_tf.poles[i].real; + const T di = factored_tf.poles[i].imag; if (dr * dr + di * di < tol * tol) { collision = true; @@ -760,8 +760,8 @@ constexpr TransferFunction matched_z_discretize_factored( } for (consteig::Size i = 0; i < nz && !collision; ++i) { - const T dr = w_c - ftf.zeros[i].real; - const T di = ftf.zeros[i].imag; + const T dr = w_c - factored_tf.zeros[i].real; + const T di = factored_tf.zeros[i].imag; if (dr * dr + di * di < tol * tol) { collision = true; @@ -782,13 +782,13 @@ constexpr TransferFunction matched_z_discretize_factored( Complex num_c_cx{static_cast(1), static_cast(0)}; for (consteig::Size i = 0; i < nz; ++i) { - num_c_cx = num_c_cx * (w_c_cx - ftf.zeros[i]); + num_c_cx = num_c_cx * (w_c_cx - factored_tf.zeros[i]); } Complex den_c_cx{static_cast(1), static_cast(0)}; for (consteig::Size i = 0; i < N; ++i) { - den_c_cx = den_c_cx * (w_c_cx - ftf.poles[i]); + den_c_cx = den_c_cx * (w_c_cx - factored_tf.poles[i]); } Complex num_d_cx{static_cast(1), static_cast(0)}; @@ -961,25 +961,25 @@ constexpr TransferFunction analog_to_digital( template constexpr TransferFunction discretize_with_factored( - const T (&b_c)[N + 1u], const T (&a_c)[N + 1u], const FactoredTF &ftf, - T Ts, ZOH) + const T (&b_c)[N + 1u], const T (&a_c)[N + 1u], + const FactoredTF &factored_tf, T Ts, ZOH) { return ss_to_tf( - zoh_discretize_with_evals(tf_to_ss(b_c, a_c), Ts, ftf)); + zoh_discretize_with_evals(tf_to_ss(b_c, a_c), Ts, factored_tf)); } template constexpr TransferFunction discretize_with_factored( const T (& /*b_c*/)[N + 1u], const T (& /*a_c*/)[N + 1u], - const FactoredTF &ftf, T Ts, MatchedZ) + const FactoredTF &factored_tf, T Ts, MatchedZ) { - return matched_z_discretize_factored(ftf, Ts); + return matched_z_discretize_factored(factored_tf, Ts); } template constexpr TransferFunction discretize_with_factored( const T (&b_c)[N + 1u], const T (&a_c)[N + 1u], - const FactoredTF & /*ftf*/, T Ts, M method_tag) + const FactoredTF & /*factored_tf*/, T Ts, M method_tag) { return analog_to_digital(b_c, a_c, Ts, method_tag); } diff --git a/include/constfilt/elliptic.hpp b/include/constfilt/elliptic.hpp index 0568365..af052af 100644 --- a/include/constfilt/elliptic.hpp +++ b/include/constfilt/elliptic.hpp @@ -513,16 +513,16 @@ class Elliptic compute_prototype_poles_zeros(q, sig0, k, poles_proto, pole_cnt, zeros_proto, zero_cnt); - FactoredTF ftf{}; - ftf.nz = zero_cnt; + FactoredTF factored_tf{}; + factored_tf.nz = zero_cnt; for (consteig::Size i = 0u; i < pole_cnt; ++i) { - ftf.poles[i] = + factored_tf.poles[i] = Complex{wc * poles_proto[i].real, wc * poles_proto[i].imag}; } for (consteig::Size i = 0u; i < zero_cnt; ++i) { - ftf.zeros[i] = + factored_tf.zeros[i] = Complex{wc * zeros_proto[i].real, wc * zeros_proto[i].imag}; } @@ -536,9 +536,10 @@ class Elliptic { ++d_b; } - ftf.gain = (d_b > N) ? static_cast(0) : b_tmp[d_b] / a_tmp[0]; + factored_tf.gain = + (d_b > N) ? static_cast(0) : b_tmp[d_b] / a_tmp[0]; - return ftf; + return factored_tf; } // HP FactoredTF: derive from LP prototype via LP-to-HP transform (s -> @@ -558,14 +559,14 @@ class Elliptic const FactoredTF lp = compute_factored_tf( norm_cutoff, ripple_db, attenuation_db, LowPass{}); - FactoredTF ftf{}; + FactoredTF factored_tf{}; // HP poles: wc / lp_pole (complex division) for (consteig::Size i = 0u; i < N; ++i) { const Complex &p = lp.poles[i]; const T denom_sq = p.real * p.real + p.imag * p.imag; - ftf.poles[i] = + factored_tf.poles[i] = Complex{wc * p.real / denom_sq, -wc * p.imag / denom_sq}; } @@ -576,15 +577,16 @@ class Elliptic { const Complex &z = lp.zeros[i]; const T denom_sq = z.real * z.real + z.imag * z.imag; - ftf.zeros[hp_nz++] = + factored_tf.zeros[hp_nz++] = Complex{wc * z.real / denom_sq, -wc * z.imag / denom_sq}; } // For odd N: LP->HP adds a zero at s=0 (from the strictly-proper LP). if (N % 2u == 1u) { - ftf.zeros[hp_nz++] = Complex{static_cast(0), static_cast(0)}; + factored_tf.zeros[hp_nz++] = + Complex{static_cast(0), static_cast(0)}; } - ftf.nz = hp_nz; + factored_tf.nz = hp_nz; // Gain from the HP polynomial TF. T b_tmp[N + 1u]{}; @@ -595,9 +597,10 @@ class Elliptic { ++d_b; } - ftf.gain = (d_b > N) ? static_cast(0) : b_tmp[d_b] / a_tmp[0]; + factored_tf.gain = + (d_b > N) ? static_cast(0) : b_tmp[d_b] / a_tmp[0]; - return ftf; + return factored_tf; } }; From dcf298be17baffc084822f5c2ef6a6dc37e3fb77 Mon Sep 17 00:00:00 2001 From: Mitchell Thompkins Date: Fri, 26 Jun 2026 02:26:51 +0000 Subject: [PATCH 07/11] address MR feedback --- include/constfilt/analog_filter.hpp | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/include/constfilt/analog_filter.hpp b/include/constfilt/analog_filter.hpp index 4b0a74a..3611079 100644 --- a/include/constfilt/analog_filter.hpp +++ b/include/constfilt/analog_filter.hpp @@ -40,13 +40,14 @@ namespace constfilt // // AnalogFilter(continuous_tf, factored_tf, sample_rate_hz, method_tag) // factored_tf - factored (pole/zero/gain) form of continuous_tf; -// ZOH and MatchedZ use the known poles to build the -// Vandermonde eigenvector matrix analytically, bypassing -// the polynomial->eigenvalue roundtrip that degrades -// accuracy at high orders. Both continuous_tf (polynomial) -// and factored_tf must be supplied; Butterworth and -// Elliptic compute both and use this path internally. -// Direct callers must also provide both forms. +// ZOH builds the Vandermonde eigenvector matrix +// analytically from the known poles, bypassing the +// polynomial->eigenvalue roundtrip that degrades accuracy +// at high orders. MatchedZ maps the supplied poles/zeros +// directly via z=exp(p*Ts), skipping companion-matrix +// root-finding. Both continuous_tf and factored_tf must +// be provided; Butterworth and Elliptic compute both and +// use this path internally. template class AnalogFilter : public Filter From d11d187784f7442b8c0c2a8f1427553cae7b1597 Mon Sep 17 00:00:00 2001 From: Mitchell Thompkins Date: Fri, 26 Jun 2026 05:38:36 +0000 Subject: [PATCH 08/11] DRY --- include/constfilt/discretize.hpp | 88 +++++++++++++------------------- 1 file changed, 36 insertions(+), 52 deletions(-) diff --git a/include/constfilt/discretize.hpp b/include/constfilt/discretize.hpp index 51661f7..4d0bbb4 100644 --- a/include/constfilt/discretize.hpp +++ b/include/constfilt/discretize.hpp @@ -119,25 +119,19 @@ template struct FactoredTF // Matrix exponential -// matrix_exp(A) via eigendecomposition: -// Ad = V * diag(exp(lam_i)) * V^{-1} -// where V = eigenvectors, lam_i = eigenvalues (complex). -// Real part extracted at the end (imaginary parts cancel for real A). +// Shared kernel: given a pre-built eigenvector matrix V and per-eigenvalue +// exponential factors, compute V * diag(exp_factors) * V^{-1} and return +// the real part. template -constexpr consteig::Matrix matrix_exp( - const consteig::Matrix &A) +constexpr consteig::Matrix spectral_matrix_exp( + const consteig::Matrix, N, N> &V, + const consteig::Complex (&exp_factors)[N]) { using Complex = consteig::Complex; using ComplexMat_NN = consteig::Matrix; using ComplexMat_N1 = consteig::Matrix; - // 1. Eigenvalues and eigenvectors - const auto evals = consteig::eigenvalues(A); // Matrix - const auto V = consteig::eigenvectors(A, evals); // Matrix - - // 2. Invert V column-by-column via LU const auto lu_V = consteig::lu(V); - ComplexMat_NN V_inv{}; for (consteig::Size col = 0; col < N; ++col) { @@ -150,23 +144,19 @@ constexpr consteig::Matrix matrix_exp( } } - // 3. Accumulate: matrix_exp = sum_i exp(lam_i) * v_i * w_i^T - // where v_i = column i of V, w_i^T = row i of V_inv ComplexMat_NN result_c{}; for (consteig::Size i = 0; i < N; ++i) { - Complex exp_lambda = consteig::exp(evals(i, 0)); for (consteig::Size r = 0; r < N; ++r) { for (consteig::Size c = 0; c < N; ++c) { result_c(r, c) = - result_c(r, c) + exp_lambda * V(r, i) * V_inv(i, c); + result_c(r, c) + exp_factors[i] * V(r, i) * V_inv(i, c); } } } - // 4. Extract real part consteig::Matrix result{}; for (consteig::Size r = 0; r < N; ++r) { @@ -178,6 +168,31 @@ constexpr consteig::Matrix matrix_exp( return result; } +// matrix_exp(A) via eigendecomposition: +// Ad = V * diag(exp(lam_i)) * V^{-1} +// where V = eigenvectors, lam_i = eigenvalues (complex). +// Real part extracted at the end (imaginary parts cancel for real A). +template +constexpr consteig::Matrix matrix_exp( + const consteig::Matrix &A) +{ + using Complex = consteig::Complex; + + // 1. Eigenvalues and eigenvectors + const auto evals = consteig::eigenvalues(A); // Matrix + const auto V = consteig::eigenvectors(A, evals); // Matrix + + // 2. Per-eigenvalue exponential factors + Complex exp_factors[N]{}; + for (consteig::Size i = 0; i < N; ++i) + { + exp_factors[i] = consteig::exp(evals(i, 0)); + } + + // 3-4. Invert V, accumulate spectral sum, extract real part + return spectral_matrix_exp(V, exp_factors); +} + // matrix_exp for controllable-canonical-form (CCF) matrices given analytic // poles. // @@ -195,9 +210,7 @@ constexpr consteig::Matrix matrix_exp_ccf( { using Complex = consteig::Complex; using ComplexMat_NN = consteig::Matrix; - using ComplexMat_N1 = consteig::Matrix; - // Build Vandermonde V: V[r][i] = poles[i]^r (unscaled analog poles) ComplexMat_NN V{}; for (consteig::Size i = 0; i < N; ++i) { @@ -209,43 +222,14 @@ constexpr consteig::Matrix matrix_exp_ccf( } } - const auto lu_V = consteig::lu(V); - ComplexMat_NN V_inv{}; - for (consteig::Size col = 0; col < N; ++col) - { - ComplexMat_N1 e_col{}; - e_col(col, 0) = Complex{static_cast(1), static_cast(0)}; - auto col_vec = consteig::lu_solve(lu_V, e_col); - for (consteig::Size row = 0; row < N; ++row) - { - V_inv(row, col) = col_vec(row, 0); - } - } - - ComplexMat_NN result_c{}; + Complex exp_factors[N]{}; for (consteig::Size i = 0; i < N; ++i) { - const Complex exp_lambda = consteig::exp(Complex{ - Ts * factored_tf.poles[i].real, Ts * factored_tf.poles[i].imag}); - for (consteig::Size r = 0; r < N; ++r) - { - for (consteig::Size c = 0; c < N; ++c) - { - result_c(r, c) = - result_c(r, c) + exp_lambda * V(r, i) * V_inv(i, c); - } - } + exp_factors[i] = consteig::exp(Complex{Ts * factored_tf.poles[i].real, + Ts * factored_tf.poles[i].imag}); } - consteig::Matrix result{}; - for (consteig::Size r = 0; r < N; ++r) - { - for (consteig::Size c = 0; c < N; ++c) - { - result(r, c) = result_c(r, c).real; - } - } - return result; + return spectral_matrix_exp(V, exp_factors); } // ZOH discretization From 9acd2a01b4a3d21f196c578fe3fce1a307ec25ae Mon Sep 17 00:00:00 2001 From: Mitchell Thompkins Date: Sat, 27 Jun 2026 00:22:35 +0000 Subject: [PATCH 09/11] refactor matched z for DRY --- include/constfilt/discretize.hpp | 317 ++++++++++--------------------- 1 file changed, 102 insertions(+), 215 deletions(-) diff --git a/include/constfilt/discretize.hpp b/include/constfilt/discretize.hpp index 4d0bbb4..79459ca 100644 --- a/include/constfilt/discretize.hpp +++ b/include/constfilt/discretize.hpp @@ -457,40 +457,23 @@ constexpr StateSpace tf_to_ss(const T (&b)[N + 1u], const T (&a)[N + 1u]) // for strictly proper systems, and matches gain at a test frequency w_c. // Reference: Octave control pkg @tf/__c2d__.m, lines 32-66. template -constexpr TransferFunction matched_z_discretize_tf( - const T (&b_c)[N + 1u], const T (&a_c)[N + 1u], T Ts, MatchedZ /*tag*/) +// Shared kernel: given analog poles, finite zeros, zero count, continuous gain, +// and sample period, maps to z-domain, pads missing zeros at z=-1, matches +// gain, and assembles the discrete TF. +template +constexpr TransferFunction matched_z_assemble( + const consteig::Complex (&poles)[N], + const consteig::Complex (&zeros)[N], consteig::Size nz, T k_c, T Ts) { using Complex = consteig::Complex; - // Step 1: count leading zeros in b_c; derive nz (finite analog zero count). - consteig::Size d_b = 0; - while (d_b <= N && b_c[d_b] == static_cast(0)) - { - ++d_b; - } - const consteig::Size nz = (d_b > N) ? 0u : N - d_b; - - // Step 2: continuous leading-coefficient gain. - const T k_c = (d_b > N) ? static_cast(0) : b_c[d_b] / a_c[0]; - - // Step 3: find analog poles (roots of a_c) via companion matrix. - // consteig exposes eigenvalues, not polynomial roots; a companion matrix - // has eigenvalues equal to the polynomial roots by construction. - consteig::Matrix A_pole{}; - for (consteig::Size i = 0; i < N - 1u; ++i) - A_pole(i + 1u, i) = static_cast(1); - for (consteig::Size i = 0; i < N; ++i) - A_pole(i, N - 1u) = -(a_c[N - i] / a_c[0]); - const auto p_c_evals = consteig::eigenvalues(A_pole); - - // Step 4: map poles to z-domain; build monic denominator polynomial. + // Step 1: map poles to z-domain; build monic denominator polynomial. Complex p_d_vals[N]{}; Complex pole_poly[N + 1u]{}; pole_poly[0] = Complex{static_cast(1), static_cast(0)}; for (consteig::Size k = 0; k < N; ++k) { - p_d_vals[k] = - consteig::exp(p_c_evals(k, 0) * Complex{Ts, static_cast(0)}); + p_d_vals[k] = consteig::exp(poles[k] * Complex{Ts, static_cast(0)}); const Complex zk = p_d_vals[k]; pole_poly[k + 1u] = Complex{static_cast(0), static_cast(0)} - zk * pole_poly[k]; @@ -500,71 +483,24 @@ constexpr TransferFunction matched_z_discretize_tf( } } - // Steps 5-7: find analog zeros (roots of b_c) via companion matrix, - // same trick as above. Embedded in NxN so consteig::eigenvalues can be - // called with a fixed size; the top-left nz x nz block holds the actual - // companion, the remaining entries are 0, producing d_b spurious - // eigenvalues at the origin that are discarded afterward. - Complex z_c_finite[N]{}; + // Step 2: map finite zeros to z-domain; build monic numerator polynomial. Complex z_d_finite[N]{}; Complex zero_poly[N + 1u]{}; zero_poly[0] = Complex{static_cast(1), static_cast(0)}; - - if (nz > 0u) + for (consteig::Size k = 0; k < nz; ++k) { - consteig::Matrix A_zero{}; - for (consteig::Size i = 0; i + 1u < nz; ++i) - A_zero(i + 1u, i) = static_cast(1); - for (consteig::Size i = 0; i < nz; ++i) - A_zero(i, nz - 1u) = -(b_c[d_b + (nz - i)] / b_c[d_b]); - const auto z_c_evals = consteig::eigenvalues(A_zero); - - // Mark d_b spurious eigenvalues (smallest magnitude = from zero block). - bool spurious[N]{}; - for (consteig::Size m = 0; m < d_b; ++m) - { - consteig::Size min_idx = 0; - T min_mag_sq = static_cast(-1); - for (consteig::Size i = 0; i < N; ++i) - { - if (!spurious[i]) - { - const T mag_sq = - z_c_evals(i, 0).real * z_c_evals(i, 0).real + - z_c_evals(i, 0).imag * z_c_evals(i, 0).imag; - if (min_mag_sq < static_cast(0) || mag_sq < min_mag_sq) - { - min_mag_sq = mag_sq; - min_idx = i; - } - } - } - spurious[min_idx] = true; - } - - // Collect finite zeros and build zero polynomial prod(z - z_d_k). - consteig::Size nz_cnt = 0u; - for (consteig::Size k = 0; k < N; ++k) + z_d_finite[k] = + consteig::exp(zeros[k] * Complex{Ts, static_cast(0)}); + const Complex zk = z_d_finite[k]; + zero_poly[k + 1u] = + Complex{static_cast(0), static_cast(0)} - zk * zero_poly[k]; + for (consteig::Size i = k; i > 0u; --i) { - if (!spurious[k]) - { - z_c_finite[nz_cnt] = z_c_evals(k, 0); - z_d_finite[nz_cnt] = consteig::exp( - z_c_evals(k, 0) * Complex{Ts, static_cast(0)}); - const Complex zk = z_d_finite[nz_cnt]; - zero_poly[nz_cnt + 1u] = - Complex{static_cast(0), static_cast(0)} - - zk * zero_poly[nz_cnt]; - for (consteig::Size i = nz_cnt; i > 0u; --i) - { - zero_poly[i] = zero_poly[i] - zk * zero_poly[i - 1u]; - } - ++nz_cnt; - } + zero_poly[i] = zero_poly[i] - zk * zero_poly[i - 1u]; } } - // Step 8: pad with zeros at z = -1 to reach numerator degree N-1. + // Step 3: pad with zeros at z = -1 to reach numerator degree N-1. const consteig::Size n_extra = (nz + 1u < N) ? (N - nz - 1u) : 0u; for (consteig::Size e = 0; e < n_extra; ++e) { @@ -577,7 +513,7 @@ constexpr TransferFunction matched_z_discretize_tf( } const consteig::Size num_deg = nz + n_extra; - // Step 9: find matching frequency w_c (avoid collision with poles/zeros). + // Step 4: find matching frequency w_c (avoid collision with poles/zeros). const T tol = gcem::sqrt(consteig::epsilon()); T w_c = static_cast(0); for (consteig::Size attempt = 0; attempt < 1000u; ++attempt) @@ -585,8 +521,8 @@ constexpr TransferFunction matched_z_discretize_tf( bool collision = false; for (consteig::Size i = 0; i < N && !collision; ++i) { - const T dr = w_c - p_c_evals(i, 0).real; - const T di = p_c_evals(i, 0).imag; + const T dr = w_c - poles[i].real; + const T di = poles[i].imag; if (dr * dr + di * di < tol * tol) { collision = true; @@ -594,8 +530,8 @@ constexpr TransferFunction matched_z_discretize_tf( } for (consteig::Size i = 0; i < nz && !collision; ++i) { - const T dr = w_c - z_c_finite[i].real; - const T di = z_c_finite[i].imag; + const T dr = w_c - zeros[i].real; + const T di = zeros[i].imag; if (dr * dr + di * di < tol * tol) { collision = true; @@ -608,7 +544,7 @@ constexpr TransferFunction matched_z_discretize_tf( w_c += static_cast(0.1) / Ts; } - // Step 10: compute discrete gain k_d matching H_d(w_d) = H_c(w_c). + // Step 5: compute discrete gain k_d matching H_d(w_d) = H_c(w_c). const Complex w_c_cx{w_c, static_cast(0)}; const Complex w_d_cx = consteig::exp(w_c_cx * Complex{Ts, static_cast(0)}); @@ -616,13 +552,13 @@ constexpr TransferFunction matched_z_discretize_tf( Complex num_c_cx{static_cast(1), static_cast(0)}; for (consteig::Size i = 0; i < nz; ++i) { - num_c_cx = num_c_cx * (w_c_cx - z_c_finite[i]); + num_c_cx = num_c_cx * (w_c_cx - zeros[i]); } Complex den_c_cx{static_cast(1), static_cast(0)}; for (consteig::Size i = 0; i < N; ++i) { - den_c_cx = den_c_cx * (w_c_cx - p_c_evals(i, 0)); + den_c_cx = den_c_cx * (w_c_cx - poles[i]); } Complex num_d_cx{static_cast(1), static_cast(0)}; @@ -651,7 +587,7 @@ constexpr TransferFunction matched_z_discretize_tf( (gain_num.real * gain_den.real + gain_num.imag * gain_den.imag) / gain_den_sq; - // Step 11: assemble output TF. + // Step 6: assemble output TF. TransferFunction tf{}; for (consteig::Size i = 0; i <= N; ++i) { @@ -666,154 +602,105 @@ constexpr TransferFunction matched_z_discretize_tf( return tf; } -// Matched-Z discretization from a FactoredTF (precomputed poles/zeros). -// Bypasses both companion-matrix eigendecompositions in -// matched_z_discretize_tf. Uses factored_tf.poles, -// factored_tf.zeros[0..factored_tf.nz-1], and factored_tf.gain directly. +// Steps 1-3: extract poles and zeros from polynomial coefficients via companion +// matrix eigenvalues, then delegate to matched_z_assemble. template -constexpr TransferFunction matched_z_discretize_factored( - const FactoredTF &factored_tf, T Ts) +constexpr TransferFunction matched_z_discretize_tf( + const T (&b_c)[N + 1u], const T (&a_c)[N + 1u], T Ts, MatchedZ /*tag*/) { using Complex = consteig::Complex; - const consteig::Size nz = factored_tf.nz; - const T k_c = factored_tf.gain; - - // Map poles to z-domain; build monic denominator polynomial. - Complex p_d_vals[N]{}; - Complex pole_poly[N + 1u]{}; - pole_poly[0] = Complex{static_cast(1), static_cast(0)}; - for (consteig::Size k = 0; k < N; ++k) + // Step 1: count leading zeros in b_c; derive nz and k_c. + consteig::Size d_b = 0; + while (d_b <= N && b_c[d_b] == static_cast(0)) { - p_d_vals[k] = consteig::exp(factored_tf.poles[k] * - Complex{Ts, static_cast(0)}); - const Complex zk = p_d_vals[k]; - pole_poly[k + 1u] = - Complex{static_cast(0), static_cast(0)} - zk * pole_poly[k]; - for (consteig::Size i = k; i > 0u; --i) - { - pole_poly[i] = pole_poly[i] - zk * pole_poly[i - 1u]; - } + ++d_b; } + const consteig::Size nz = (d_b > N) ? 0u : N - d_b; + const T k_c = (d_b > N) ? static_cast(0) : b_c[d_b] / a_c[0]; - // Map finite zeros to z-domain; build monic numerator polynomial. - // No companion matrix or spurious-eigenvalue filtering needed. - Complex z_d_finite[N]{}; - Complex zero_poly[N + 1u]{}; - zero_poly[0] = Complex{static_cast(1), static_cast(0)}; - for (consteig::Size k = 0; k < nz; ++k) + // Step 2: find analog poles (roots of a_c) via companion matrix. + // consteig exposes eigenvalues, not polynomial roots; a companion matrix + // has eigenvalues equal to the polynomial roots by construction. + consteig::Matrix A_pole{}; + for (consteig::Size i = 0; i < N - 1u; ++i) { - z_d_finite[k] = consteig::exp(factored_tf.zeros[k] * - Complex{Ts, static_cast(0)}); - const Complex zk = z_d_finite[k]; - zero_poly[k + 1u] = - Complex{static_cast(0), static_cast(0)} - zk * zero_poly[k]; - for (consteig::Size i = k; i > 0u; --i) - { - zero_poly[i] = zero_poly[i] - zk * zero_poly[i - 1u]; - } + A_pole(i + 1u, i) = static_cast(1); + } + for (consteig::Size i = 0; i < N; ++i) + { + A_pole(i, N - 1u) = -(a_c[N - i] / a_c[0]); + } + const auto p_c_evals = consteig::eigenvalues(A_pole); + Complex poles[N]{}; + for (consteig::Size i = 0; i < N; ++i) + { + poles[i] = p_c_evals(i, 0); } - // Pad with zeros at z = -1 to reach numerator degree N-1. - const consteig::Size n_extra = (nz + 1u < N) ? (N - nz - 1u) : 0u; - for (consteig::Size e = 0; e < n_extra; ++e) + // Step 3: find analog zeros (roots of b_c) via companion matrix. + // Embedded in NxN so consteig::eigenvalues can be called with a fixed + // size; the top-left nz x nz block holds the actual companion, the + // remaining entries are 0, producing d_b spurious eigenvalues at the + // origin that are discarded afterward. + Complex zeros[N]{}; + if (nz > 0u) { - const consteig::Size cur_deg = nz + e; - zero_poly[cur_deg + 1u] = zero_poly[cur_deg + 1u] + zero_poly[cur_deg]; - for (consteig::Size i = cur_deg; i > 0u; --i) + consteig::Matrix A_zero{}; + for (consteig::Size i = 0; i + 1u < nz; ++i) { - zero_poly[i] = zero_poly[i] + zero_poly[i - 1u]; + A_zero(i + 1u, i) = static_cast(1); } - } - const consteig::Size num_deg = nz + n_extra; + for (consteig::Size i = 0; i < nz; ++i) + { + A_zero(i, nz - 1u) = -(b_c[d_b + (nz - i)] / b_c[d_b]); + } + const auto z_c_evals = consteig::eigenvalues(A_zero); - // Find matching frequency w_c (avoid collision with poles/zeros). - const T tol = gcem::sqrt(consteig::epsilon()); - T w_c = static_cast(0); - for (consteig::Size attempt = 0; attempt < 1000u; ++attempt) - { - bool collision = false; - for (consteig::Size i = 0; i < N && !collision; ++i) + // Mark d_b spurious eigenvalues (smallest magnitude = from zero block). + bool spurious[N]{}; + for (consteig::Size m = 0; m < d_b; ++m) { - const T dr = w_c - factored_tf.poles[i].real; - const T di = factored_tf.poles[i].imag; - if (dr * dr + di * di < tol * tol) + consteig::Size min_idx = 0; + T min_mag_sq = static_cast(-1); + for (consteig::Size i = 0; i < N; ++i) { - collision = true; + if (!spurious[i]) + { + const T mag_sq = + z_c_evals(i, 0).real * z_c_evals(i, 0).real + + z_c_evals(i, 0).imag * z_c_evals(i, 0).imag; + if (min_mag_sq < static_cast(0) || mag_sq < min_mag_sq) + { + min_mag_sq = mag_sq; + min_idx = i; + } + } } + spurious[min_idx] = true; } - for (consteig::Size i = 0; i < nz && !collision; ++i) + + consteig::Size nz_cnt = 0u; + for (consteig::Size k = 0; k < N; ++k) { - const T dr = w_c - factored_tf.zeros[i].real; - const T di = factored_tf.zeros[i].imag; - if (dr * dr + di * di < tol * tol) + if (!spurious[k]) { - collision = true; + zeros[nz_cnt++] = z_c_evals(k, 0); } } - if (!collision) - { - break; - } - w_c += static_cast(0.1) / Ts; } - // Compute discrete gain k_d matching H_d(w_d) = H_c(w_c). - const Complex w_c_cx{w_c, static_cast(0)}; - const Complex w_d_cx = - consteig::exp(w_c_cx * Complex{Ts, static_cast(0)}); - - Complex num_c_cx{static_cast(1), static_cast(0)}; - for (consteig::Size i = 0; i < nz; ++i) - { - num_c_cx = num_c_cx * (w_c_cx - factored_tf.zeros[i]); - } - - Complex den_c_cx{static_cast(1), static_cast(0)}; - for (consteig::Size i = 0; i < N; ++i) - { - den_c_cx = den_c_cx * (w_c_cx - factored_tf.poles[i]); - } - - Complex num_d_cx{static_cast(1), static_cast(0)}; - for (consteig::Size i = 0; i < N; ++i) - { - num_d_cx = num_d_cx * (w_d_cx - p_d_vals[i]); - } - - Complex den_d_cx{static_cast(1), static_cast(0)}; - for (consteig::Size i = 0; i < nz; ++i) - { - den_d_cx = den_d_cx * (w_d_cx - z_d_finite[i]); - } - for (consteig::Size e = 0; e < n_extra; ++e) - { - den_d_cx = den_d_cx * - (w_d_cx - Complex{static_cast(-1), static_cast(0)}); - } - - const Complex gain_num = - Complex{k_c, static_cast(0)} * num_c_cx * num_d_cx; - const Complex gain_den = den_c_cx * den_d_cx; - const T gain_den_sq = - gain_den.real * gain_den.real + gain_den.imag * gain_den.imag; - const T k_d = - (gain_num.real * gain_den.real + gain_num.imag * gain_den.imag) / - gain_den_sq; - - // Assemble output TF. - TransferFunction tf{}; - for (consteig::Size i = 0; i <= N; ++i) - { - tf.a[i] = pole_poly[i].real; - } - const consteig::Size pad = N - num_deg; - for (consteig::Size i = 0; i <= num_deg; ++i) - { - tf.b[pad + i] = k_d * zero_poly[i].real; - } + return matched_z_assemble(poles, zeros, nz, k_c, Ts); +} - return tf; +// Matched-Z discretization from a FactoredTF: poles and zeros are known +// analytically, so companion-matrix eigendecompositions are not needed. +template +constexpr TransferFunction matched_z_discretize_factored( + const FactoredTF &factored_tf, T Ts) +{ + return matched_z_assemble(factored_tf.poles, factored_tf.zeros, + factored_tf.nz, factored_tf.gain, Ts); } // Tustin (bilinear) discretization From bb20ca2ec3d1fff6761cf3ec52a1179bfacc8e1b Mon Sep 17 00:00:00 2001 From: Mitchell Thompkins Date: Sat, 27 Jun 2026 06:06:52 +0000 Subject: [PATCH 10/11] code cleanup --- include/constfilt/analog_filter.hpp | 12 +++--------- include/constfilt/discretize.hpp | 9 ++++----- 2 files changed, 7 insertions(+), 14 deletions(-) diff --git a/include/constfilt/analog_filter.hpp b/include/constfilt/analog_filter.hpp index 3611079..ca4913f 100644 --- a/include/constfilt/analog_filter.hpp +++ b/include/constfilt/analog_filter.hpp @@ -39,15 +39,9 @@ namespace constfilt // constfilt::prewarp(warp_hz) // // AnalogFilter(continuous_tf, factored_tf, sample_rate_hz, method_tag) -// factored_tf - factored (pole/zero/gain) form of continuous_tf; -// ZOH builds the Vandermonde eigenvector matrix -// analytically from the known poles, bypassing the -// polynomial->eigenvalue roundtrip that degrades accuracy -// at high orders. MatchedZ maps the supplied poles/zeros -// directly via z=exp(p*Ts), skipping companion-matrix -// root-finding. Both continuous_tf and factored_tf must -// be provided; Butterworth and Elliptic compute both and -// use this path internally. +// factored_tf - the poles, zeros, and gain of continuous_tf in +// factored form. Required alongside continuous_tf. +// Used internally by Butterworth and Elliptic. template class AnalogFilter : public Filter diff --git a/include/constfilt/discretize.hpp b/include/constfilt/discretize.hpp index 79459ca..666154b 100644 --- a/include/constfilt/discretize.hpp +++ b/include/constfilt/discretize.hpp @@ -452,11 +452,6 @@ constexpr StateSpace tf_to_ss(const T (&b)[N + 1u], const T (&a)[N + 1u]) // Matched-Z discretization (TF entry point) -// Full matched-Z from analog transfer function coefficients. -// Maps each finite analog zero via z = exp(s*Ts), pads with zeros at z = -1 -// for strictly proper systems, and matches gain at a test frequency w_c. -// Reference: Octave control pkg @tf/__c2d__.m, lines 32-66. -template // Shared kernel: given analog poles, finite zeros, zero count, continuous gain, // and sample period, maps to z-domain, pads missing zeros at z=-1, matches // gain, and assembles the discrete TF. @@ -602,6 +597,10 @@ constexpr TransferFunction matched_z_assemble( return tf; } +// Full matched-Z from analog transfer function coefficients. +// Maps each finite analog zero via z = exp(s*Ts), pads with zeros at z = -1 +// for strictly proper systems, and matches gain at a test frequency w_c. +// Reference: Octave control pkg @tf/__c2d__.m, lines 32-66. // Steps 1-3: extract poles and zeros from polynomial coefficients via companion // matrix eigenvalues, then delegate to matched_z_assemble. template From e7f844b4bc912b521a825d046539253211ee1658 Mon Sep 17 00:00:00 2001 From: Mitchell Thompkins Date: Sat, 27 Jun 2026 15:46:48 +0000 Subject: [PATCH 11/11] mr feedback --- include/constfilt/analog_filter.hpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/include/constfilt/analog_filter.hpp b/include/constfilt/analog_filter.hpp index ca4913f..f6ce156 100644 --- a/include/constfilt/analog_filter.hpp +++ b/include/constfilt/analog_filter.hpp @@ -41,7 +41,9 @@ namespace constfilt // AnalogFilter(continuous_tf, factored_tf, sample_rate_hz, method_tag) // factored_tf - the poles, zeros, and gain of continuous_tf in // factored form. Required alongside continuous_tf. -// Used internally by Butterworth and Elliptic. +// For ZOH, the poles are used to build the Vandermonde +// eigenvector matrix analytically. For MatchedZ, the +// poles and zeros are mapped directly to the z-domain. template class AnalogFilter : public Filter