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 new file mode 100644 index 0000000..55c616b --- /dev/null +++ b/docs/profiling/results/accuracy_gcc_15.2.0.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:cfa95043df54bef12e6843007eb4561f5f086a51eb90afc46c72febefc5d4868 +size 168736 diff --git a/docs/profiling/results/accuracy_gcc_15.2.0_butterworth.png b/docs/profiling/results/accuracy_gcc_15.2.0_butterworth.png index 5062f56..318c825 100644 --- a/docs/profiling/results/accuracy_gcc_15.2.0_butterworth.png +++ b/docs/profiling/results/accuracy_gcc_15.2.0_butterworth.png @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:cc8583c6a71f1637c933bf30a98f89d1e771ec8c70dff1152e3858488738b916 -size 87468 +oid sha256:6c13882e257e20acced0dacd4c846502b97c51d6711504f567046c3b2614f758 +size 88239 diff --git a/docs/profiling/results/accuracy_gcc_15.2.0_constfilt.png b/docs/profiling/results/accuracy_gcc_15.2.0_constfilt.png index 3ac4bc4..a5aca04 100644 --- a/docs/profiling/results/accuracy_gcc_15.2.0_constfilt.png +++ b/docs/profiling/results/accuracy_gcc_15.2.0_constfilt.png @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:1ed9672c26c959f27f31dc3d5db97d30510e4b96da51b34c51e4f9636d9284a9 -size 135871 +oid sha256:3d1cc27b92e2d4220a24f1cfe4fe60953760751380bc8a1d5ae49233a1b33cd7 +size 132095 diff --git a/docs/profiling/results/accuracy_gcc_15.2.0_elliptic.png b/docs/profiling/results/accuracy_gcc_15.2.0_elliptic.png index 739e30a..a2e3e29 100644 --- a/docs/profiling/results/accuracy_gcc_15.2.0_elliptic.png +++ b/docs/profiling/results/accuracy_gcc_15.2.0_elliptic.png @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:27c3e74b25db2fac6a5c77dd10b2ec434f665653629d5ba0dac403ad9134dbb1 -size 93862 +oid sha256:d40ad7b9ee5a121698312d60047748877b837538274c732e20b9753ee41fc8b0 +size 93791 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..fa5286a 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,38664,0 +butterworth,10,tustin,0.45,84572,0 +butterworth,10,zoh,0.54,96292,0 +butterworth,11,matchedz,0.07,38672,0 +butterworth,11,tustin,0.61,104092,0 +butterworth,11,zoh,0.72,119492,0 +butterworth,12,matchedz,0.07,38888,0 +butterworth,12,tustin,0.80,127832,0 +butterworth,12,zoh,0.95,146800,0 +butterworth,1,matchedz,0.06,37560,0 +butterworth,1,tustin,0.06,38080,0 +butterworth,1,zoh,0.07,38176,0 +butterworth,2,matchedz,0.06,37732,0 +butterworth,2,tustin,0.08,38680,0 +butterworth,2,zoh,0.07,39148,0 +butterworth,4,matchedz,0.06,38044,0 +butterworth,4,tustin,0.09,41548,0 +butterworth,4,zoh,0.10,42804,0 +butterworth,6,matchedz,0.06,38096,0 +butterworth,6,tustin,0.14,47264,0 +butterworth,6,zoh,0.16,50184,0 +butterworth,8,matchedz,0.07,38272,0 +butterworth,8,tustin,0.24,60016,0 +butterworth,8,zoh,0.30,66640,0 +butterworth,9,matchedz,0.07,38508,0 +butterworth,9,tustin,0.34,71000,0 +butterworth,9,zoh,0.40,80012,0 +elliptic,10,matchedz,0.28,56052,0 +elliptic,10,tustin,0.66,101920,0 +elliptic,10,zoh,0.76,114416,0 +elliptic,11,matchedz,0.27,54756,0 +elliptic,11,tustin,0.82,120820,0 +elliptic,11,zoh,0.94,135896,0 +elliptic,12,matchedz,0.29,57064,0 +elliptic,12,tustin,1.05,147304,0 +elliptic,12,zoh,1.20,166836,0 +elliptic,2,matchedz,0.16,46144,0 +elliptic,2,tustin,0.17,47260,0 +elliptic,2,zoh,0.17,47724,0 +elliptic,4,matchedz,0.19,49120,0 +elliptic,4,tustin,0.22,52524,0 +elliptic,4,zoh,0.23,53976,0 +elliptic,6,matchedz,0.21,50788,0 +elliptic,6,tustin,0.29,59564,0 +elliptic,6,zoh,0.31,62484,0 +elliptic,8,matchedz,0.25,53068,0 +elliptic,8,tustin,0.43,74600,0 +elliptic,8,zoh,0.48,81180,0 +elliptic,9,matchedz,0.26,54128,0 +elliptic,9,tustin,0.52,86224,0 +elliptic,9,zoh,0.59,95288,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..5b8b3c2 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:1c4cef175d44ca0f176b4da87356b1a7bc0113726a66806f9785886ce5f96c20 +size 132964 diff --git a/docs/profiling/results/runtime_gcc_15.2.0.csv b/docs/profiling/results/runtime_gcc_15.2.0.csv index 5b7cc99..8805a0d 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.406,156.1,1.00000000 +constfilt,Butterworth,2,tustin,10.824,92.4,1.00000000 +constfilt,Butterworth,4,tustin,18.803,53.2,1.00000000 +constfilt,Butterworth,6,tustin,26.040,38.4,1.00000000 +constfilt,Butterworth,8,tustin,34.177,29.3,1.00000000 +constfilt,Butterworth,1,zoh,6.405,156.1,1.00000000 +constfilt,Butterworth,2,zoh,10.823,92.4,1.00000000 +constfilt,Butterworth,4,zoh,18.407,54.3,1.00000000 +constfilt,Butterworth,6,zoh,25.916,38.6,1.00000000 +constfilt,Butterworth,8,zoh,34.118,29.3,1.00000000 +constfilt,Butterworth,1,matchedz,6.404,156.1,1.00000000 +constfilt,Butterworth,2,matchedz,10.821,92.4,1.00000000 +constfilt,Butterworth,4,matchedz,18.703,53.5,1.00000000 +constfilt,Butterworth,6,matchedz,26.001,38.5,1.00000000 +constfilt,Butterworth,8,matchedz,34.263,29.2,1.00000000 +constfilt,Elliptic,2,tustin,10.847,92.2,0.94406088 +constfilt,Elliptic,4,tustin,18.317,54.6,0.94406088 +constfilt,Elliptic,6,tustin,25.925,38.6,0.94406088 +constfilt,Elliptic,8,tustin,34.192,29.2,0.94406088 +constfilt,Elliptic,2,zoh,10.809,92.5,0.94406088 +constfilt,Elliptic,4,zoh,18.345,54.5,0.94406088 +constfilt,Elliptic,6,zoh,25.902,38.6,0.94406088 +constfilt,Elliptic,8,zoh,34.761,28.8,0.94406088 +constfilt,Elliptic,2,matchedz,11.046,90.5,0.94406088 +constfilt,Elliptic,4,matchedz,18.812,53.2,0.94406088 +constfilt,Elliptic,6,matchedz,25.954,38.5,0.94406088 +constfilt,Elliptic,8,matchedz,34.177,29.3,0.94406088 +iir1,Butterworth,2,runtime,10.599,94.4,1.00000000 +iir1,Butterworth,4,runtime,18.020,55.5,1.00000000 +iir1,Butterworth,6,runtime,28.847,34.7,1.00000000 +iir1,Butterworth,8,runtime,39.983,25.0,1.00000000 +kfr,Butterworth,2,runtime+simd,602.752,1.7,1.00000000 +kfr,Butterworth,4,runtime+simd,386.308,2.6,1.00000000 +kfr,Butterworth,6,runtime+simd,307.883,3.2,1.00000000 +kfr,Butterworth,8,runtime+simd,307.968,3.2,1.00000000 +kfr,Elliptic,2,runtime+simd,602.510,1.7,0.94406088 +kfr,Elliptic,4,runtime+simd,386.593,2.6,0.94406088 +kfr,Elliptic,6,runtime+simd,308.161,3.2,0.94406088 +kfr,Elliptic,8,runtime+simd,308.193,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 new file mode 100644 index 0000000..4a90f62 --- /dev/null +++ b/docs/profiling/results/runtime_gcc_15.2.0.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:897a23f8a1483607f16882fd33041d0585b73c88ff2993661cffff2a07ec50ee +size 111057 diff --git a/docs/profiling/results/runtime_gcc_15.2.0_butterworth.png b/docs/profiling/results/runtime_gcc_15.2.0_butterworth.png index 9b0aea2..5033a33 100644 --- a/docs/profiling/results/runtime_gcc_15.2.0_butterworth.png +++ b/docs/profiling/results/runtime_gcc_15.2.0_butterworth.png @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:f2addefa51e3cfde14fa556a03b9e180e9b7c2db46e5ea97c4d5b12a4b9d2fad -size 86772 +oid sha256:cdca7b8ddc070712bd7ba47d5143cdbc2d6cba0c90ce935452121823402c2677 +size 85786 diff --git a/docs/profiling/results/runtime_gcc_15.2.0_constfilt.png b/docs/profiling/results/runtime_gcc_15.2.0_constfilt.png index d38a16c..82085d0 100644 --- a/docs/profiling/results/runtime_gcc_15.2.0_constfilt.png +++ b/docs/profiling/results/runtime_gcc_15.2.0_constfilt.png @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:145d2ba5c303edad77bd65998f5377b1903c090c2147e761717ff8b35f3b0b83 -size 99376 +oid sha256:fa4541286340609555fa4bd86e1d8321007ece63c589a606d3244a386579aa7a +size 99325 diff --git a/docs/profiling/results/runtime_gcc_15.2.0_elliptic.png b/docs/profiling/results/runtime_gcc_15.2.0_elliptic.png index fcfbd5c..77e6436 100644 --- a/docs/profiling/results/runtime_gcc_15.2.0_elliptic.png +++ b/docs/profiling/results/runtime_gcc_15.2.0_elliptic.png @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:43cad1a36f040cc2acb913be8bf19104d20d0aeac9aac7842e50c5223df38610 -size 74664 +oid sha256:2510f139d8e3a6935c7a2ee6b17a5f31712a19cac8dc6671be78df05c846967a +size 74021 diff --git a/include/constfilt/analog_filter.hpp b/include/constfilt/analog_filter.hpp index 11b19e9..f6ce156 100644 --- a/include/constfilt/analog_filter.hpp +++ b/include/constfilt/analog_filter.hpp @@ -37,6 +37,13 @@ namespace constfilt // AnalogFilter(continuous_tf, sample_rate_hz, method_tag) // method_tag - method instance; for TustinPW supply // constfilt::prewarp(warp_hz) +// +// 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. +// 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 @@ -79,6 +86,15 @@ class AnalogFilter : public Filter { } + constexpr AnalogFilter(TransferFunction continuous_tf, + const FactoredTF &factored_tf, + T sample_rate_hz, BoundMethod method_tag) + : AnalogFilter(checked_discretize_factored(continuous_tf.b, + continuous_tf.a, factored_tf, + sample_rate_hz, method_tag)) + { + } + private: constexpr explicit AnalogFilter( TransferFunction digital_tf) @@ -120,6 +136,16 @@ 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 &factored_tf, + T sample_rate_hz, BoundMethod method_tag) + { + return discretize_with_factored( + b_c, a_c, factored_tf, static_cast(1) / sample_rate_hz, + method_tag); + } }; } // namespace constfilt diff --git a/include/constfilt/butterworth.hpp b/include/constfilt/butterworth.hpp index b50d4e5..d150a01 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,45 @@ 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 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); + factored_tf.poles[k - 1u] = {wc * gcem::cos(theta), + wc * gcem::sin(theta)}; + } + return factored_tf; + } + + // 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 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); + 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 factored_tf; + } + // 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..666154b 100644 --- a/include/constfilt/discretize.hpp +++ b/include/constfilt/discretize.hpp @@ -103,27 +103,35 @@ 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: -// 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) { @@ -136,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) { @@ -164,6 +168,70 @@ 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. +// +// 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_ccf( + const FactoredTF &factored_tf, T Ts) +{ + using Complex = consteig::Complex; + using ComplexMat_NN = consteig::Matrix; + + 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 * factored_tf.poles[i]; + } + } + + Complex exp_factors[N]{}; + for (consteig::Size i = 0; i < N; ++i) + { + exp_factors[i] = consteig::exp(Complex{Ts * factored_tf.poles[i].real, + Ts * factored_tf.poles[i].imag}); + } + + return spectral_matrix_exp(V, exp_factors); +} + // ZOH discretization // ZOH: Ad = matrix_exp(Ac*Ts), Bd = Ac^{-1} * (Ad - I) * Bc @@ -200,6 +268,28 @@ constexpr StateSpace zoh_discretize(const StateSpace &sys_c, T Ts, return sys_d; } +// ZOH discretization using analytically known poles from a FactoredTF. +// 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 &factored_tf) +{ + 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; + 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: @@ -362,45 +452,23 @@ 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. +// 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_discretize_tf( - const T (&b_c)[N + 1u], const T (&a_c)[N + 1u], T Ts, MatchedZ /*tag*/) +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]; @@ -410,71 +478,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) { @@ -487,7 +508,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) @@ -495,8 +516,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; @@ -504,8 +525,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; @@ -518,7 +539,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)}); @@ -526,13 +547,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)}; @@ -561,7 +582,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) { @@ -576,6 +597,111 @@ constexpr TransferFunction matched_z_discretize_tf( 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 +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; + + // 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)) + { + ++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]; + + // 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) + { + 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); + } + + // 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) + { + 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; + } + + consteig::Size nz_cnt = 0u; + for (consteig::Size k = 0; k < N; ++k) + { + if (!spurious[k]) + { + zeros[nz_cnt++] = z_c_evals(k, 0); + } + } + } + + return matched_z_assemble(poles, zeros, nz, k_c, Ts); +} + +// 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 // // Parameterized by alpha = 2/Ts (standard) or wc/tan(wc*Ts/2) (prewarped). @@ -699,6 +825,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 &factored_tf, T Ts, ZOH) +{ + return ss_to_tf( + 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 &factored_tf, T Ts, MatchedZ) +{ + 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 & /*factored_tf*/, 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..af052af 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,6 +308,45 @@ 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: @@ -449,6 +490,118 @@ 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 factored_tf{}; + factored_tf.nz = zero_cnt; + for (consteig::Size i = 0u; i < pole_cnt; ++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) + { + factored_tf.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; + } + factored_tf.gain = + (d_b > N) ? static_cast(0) : b_tmp[d_b] / a_tmp[0]; + + return factored_tf; + } + + // 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 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; + factored_tf.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; + 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) + { + factored_tf.zeros[hp_nz++] = + Complex{static_cast(0), static_cast(0)}; + } + factored_tf.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; + } + factored_tf.gain = + (d_b > N) ? static_cast(0) : b_tmp[d_b] / a_tmp[0]; + + return factored_tf; + } }; } // namespace constfilt