boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Code Issues Projects Releases Wiki Activity

Cardinal trigonometric: Fix bug on last coefficient for an even number of samples.

Browse Source
This commit is contained in:
NAThompson
2019-07-28 17:07:13 -04:00
parent 31fa421c46
commit 3277587ebc
3 changed files with 136 additions and 36 deletions
Download Patch File Download Diff File Expand all files Collapse all files

48
doc/graphs/fourier_bump.svg Normal file
View File

@@ -0,0 +1,48 @@
<?xml version="1.0" encoding='UTF-8' ?>
<svg xmlns='http://www.w3.org/2000/svg' width='1100' height='679'>
<style>svg { background-color: black; }
</style>
<text x='550' y='20' font-family='Palatino' font-size='25' fill='white' alignment-baseline='middle' text-anchor='middle'>A bump function (blue) interpolated via Fourier series (orange) created from 10 samples</text>
<g transform='translate(25, 40)'>
<line x1='0' y1='0' x2='0' y2='619' stroke='gray' stroke-width='1' />
<line x1='0' y1='619' x2='1055' y2='619' stroke='gray' stroke-width='1' />
<line x1='0' y1='541.625' x2='1055' y2='541.625' stroke='gray' stroke-width='1' opacity='0.5' stroke-dasharray='4' />
<text x='-7' y='538.625' font-family='times' font-size='10' fill='white' transform='rotate(-90 -1 546.625)'>0.04598</text>
<line x1='0' y1='464.2' x2='1055' y2='464.2' stroke='gray' stroke-width='1' opacity='0.5' stroke-dasharray='4' />
<text x='-7' y='461.2' font-family='times' font-size='10' fill='white' transform='rotate(-90 -1 469.2)'>0.09197</text>
<line x1='0' y1='386.9' x2='1055' y2='386.9' stroke='gray' stroke-width='1' opacity='0.5' stroke-dasharray='4' />
<text x='-7' y='383.9' font-family='times' font-size='10' fill='white' transform='rotate(-90 -1 391.9)'>0.138</text>
<line x1='0' y1='309.5' x2='1055' y2='309.5' stroke='gray' stroke-width='1' opacity='0.5' stroke-dasharray='4' />
<text x='-7' y='306.5' font-family='times' font-size='10' fill='white' transform='rotate(-90 -1 314.5)'>0.1839</text>
<line x1='0' y1='232.1' x2='1055' y2='232.1' stroke='gray' stroke-width='1' opacity='0.5' stroke-dasharray='4' />
<text x='-7' y='229.1' font-family='times' font-size='10' fill='white' transform='rotate(-90 -1 237.1)'>0.2299</text>
<line x1='0' y1='154.8' x2='1055' y2='154.8' stroke='gray' stroke-width='1' opacity='0.5' stroke-dasharray='4' />
<text x='-7' y='151.8' font-family='times' font-size='10' fill='white' transform='rotate(-90 -1 159.8)'>0.2759</text>
<line x1='0' y1='77.38' x2='1055' y2='77.38' stroke='gray' stroke-width='1' opacity='0.5' stroke-dasharray='4' />
<text x='-7' y='74.38' font-family='times' font-size='10' fill='white' transform='rotate(-90 -1 82.38)'>0.3219</text>
<line x1='0' y1='0' x2='1055' y2='0' stroke='gray' stroke-width='1' opacity='0.5' stroke-dasharray='4' />
<text x='-7' y='-3' font-family='times' font-size='10' fill='white' transform='rotate(-90 -1 5)'>0.3679</text>
<line x1='105.5' y1='0' x2='105.5' y2='619' stroke='gray' stroke-width='1' opacity='0.5' stroke-dasharray='4' />
<text x='95.5' y='629' font-family='times' font-size='10' fill='white'>-0.8</text>
<line x1='211' y1='0' x2='211' y2='619' stroke='gray' stroke-width='1' opacity='0.5' stroke-dasharray='4' />
<text x='201' y='629' font-family='times' font-size='10' fill='white'>-0.6</text>
<line x1='316.5' y1='0' x2='316.5' y2='619' stroke='gray' stroke-width='1' opacity='0.5' stroke-dasharray='4' />
<text x='306.5' y='629' font-family='times' font-size='10' fill='white'>-0.4</text>
<line x1='422' y1='0' x2='422' y2='619' stroke='gray' stroke-width='1' opacity='0.5' stroke-dasharray='4' />
<text x='412' y='629' font-family='times' font-size='10' fill='white'>-0.2</text>
<line x1='527.5' y1='0' x2='527.5' y2='619' stroke='gray' stroke-width='1' opacity='0.5' stroke-dasharray='4' />
<text x='517.5' y='629' font-family='times' font-size='10' fill='white'>5.551e-17</text>
<line x1='633' y1='0' x2='633' y2='619' stroke='gray' stroke-width='1' opacity='0.5' stroke-dasharray='4' />
<text x='623' y='629' font-family='times' font-size='10' fill='white'>0.2</text>
<line x1='738.5' y1='0' x2='738.5' y2='619' stroke='gray' stroke-width='1' opacity='0.5' stroke-dasharray='4' />
<text x='728.5' y='629' font-family='times' font-size='10' fill='white'>0.4</text>
<line x1='844' y1='0' x2='844' y2='619' stroke='gray' stroke-width='1' opacity='0.5' stroke-dasharray='4' />
<text x='834' y='629' font-family='times' font-size='10' fill='white'>0.6</text>
<line x1='949.5' y1='0' x2='949.5' y2='619' stroke='gray' stroke-width='1' opacity='0.5' stroke-dasharray='4' />
<text x='939.5' y='629' font-family='times' font-size='10' fill='white'>0.8</text>
<line x1='1055' y1='0' x2='1055' y2='619' stroke='gray' stroke-width='1' opacity='0.5' stroke-dasharray='4' />
<text x='1045' y='629' font-family='times' font-size='10' fill='white'>1</text>
<path d='M0 619 L4.137 619 L8.275 619 L12.41 619 L16.55 619 L20.69 619 L24.82 619 L28.96 618.9 L33.1 618.6 L37.24 617.9 L41.37 616.8 L45.51 615.1 L49.65 612.6 L53.78 609.4 L57.92 605.4 L62.06 600.6 L66.2 595 L70.33 588.7 L74.47 581.8 L78.61 574.2 L82.75 566.1 L86.88 557.4 L91.02 548.4 L95.16 539 L99.29 529.4 L103.4 519.4 L107.6 509.3 L111.7 499 L115.8 488.6 L120 478.1 L124.1 467.6 L128.3 457.1 L132.4 446.6 L136.5 436.1 L140.7 425.6 L144.8 415.2 L148.9 405 L153.1 394.8 L157.2 384.7 L161.4 374.7 L165.5 364.9 L169.6 355.2 L173.8 345.6 L177.9 336.2 L182 326.9 L186.2 317.8 L190.3 308.8 L194.5 300 L198.6 291.3 L202.7 282.8 L206.9 274.5 L211 266.3 L215.1 258.3 L219.3 250.4 L223.4 242.7 L227.5 235.1 L231.7 227.7 L235.8 220.5 L240 213.4 L244.1 206.5 L248.2 199.7 L252.4 193 L256.5 186.5 L260.6 180.1 L264.8 173.9 L268.9 167.8 L273.1 161.9 L277.2 156.1 L281.3 150.4 L285.5 144.9 L289.6 139.5 L293.7 134.2 L297.9 129 L302 124 L306.2 119.1 L310.3 114.3 L314.4 109.6 L318.6 105.1 L322.7 100.7 L326.8 96.34 L331 92.13 L335.1 88.03 L339.3 84.05 L343.4 80.17 L347.5 76.39 L351.7 72.73 L355.8 69.16 L359.9 65.7 L364.1 62.34 L368.2 59.08 L372.4 55.92 L376.5 52.86 L380.6 49.89 L384.8 47.01 L388.9 44.23 L393 41.55 L397.2 38.95 L401.3 36.45 L405.5 34.03 L409.6 31.71 L413.7 29.47 L417.9 27.32 L422 25.25 L426.1 23.28 L430.3 21.38 L434.4 19.57 L438.5 17.85 L442.7 16.2 L446.8 14.64 L451 13.16 L455.1 11.76 L459.2 10.45 L463.4 9.209 L467.5 8.05 L471.6 6.971 L475.8 5.971 L479.9 5.049 L484.1 4.205 L488.2 3.438 L492.3 2.75 L496.5 2.138 L500.6 1.603 L504.7 1.145 L508.9 0.7639 L513 0.4589 L517.2 0.2303 L521.3 0.07798 L525.4 0.001815 L529.6 0.001815 L533.7 0.07798 L537.8 0.2303 L542 0.4589 L546.1 0.7639 L550.3 1.145 L554.4 1.603 L558.5 2.138 L562.7 2.75 L566.8 3.438 L570.9 4.205 L575.1 5.049 L579.2 5.971 L583.4 6.971 L587.5 8.05 L591.6 9.209 L595.8 10.45 L599.9 11.76 L604 13.16 L608.2 14.64 L612.3 16.2 L616.5 17.85 L620.6 19.57 L624.7 21.38 L628.9 23.28 L633 25.25 L637.1 27.32 L641.3 29.47 L645.4 31.71 L649.5 34.03 L653.7 36.45 L657.8 38.95 L662 41.55 L666.1 44.23 L670.2 47.01 L674.4 49.89 L678.5 52.86 L682.6 55.92 L686.8 59.08 L690.9 62.34 L695.1 65.7 L699.2 69.16 L703.3 72.73 L707.5 76.39 L711.6 80.17 L715.7 84.05 L719.9 88.03 L724 92.13 L728.2 96.34 L732.3 100.7 L736.4 105.1 L740.6 109.6 L744.7 114.3 L748.8 119.1 L753 124 L757.1 129 L761.3 134.2 L765.4 139.5 L769.5 144.9 L773.7 150.4 L777.8 156.1 L781.9 161.9 L786.1 167.8 L790.2 173.9 L794.4 180.1 L798.5 186.5 L802.6 193 L806.8 199.7 L810.9 206.5 L815 213.4 L819.2 220.5 L823.3 227.7 L827.5 235.1 L831.6 242.7 L835.7 250.4 L839.9 258.3 L844 266.3 L848.1 274.5 L852.3 282.8 L856.4 291.3 L860.5 300 L864.7 308.8 L868.8 317.8 L873 326.9 L877.1 336.2 L881.2 345.6 L885.4 355.2 L889.5 364.9 L893.6 374.7 L897.8 384.7 L901.9 394.8 L906.1 405 L910.2 415.2 L914.3 425.6 L918.5 436.1 L922.6 446.6 L926.7 457.1 L930.9 467.6 L935 478.1 L939.2 488.6 L943.3 499 L947.4 509.3 L951.6 519.4 L955.7 529.4 L959.8 539 L964 548.4 L968.1 557.4 L972.3 566.1 L976.4 574.2 L980.5 581.8 L984.7 588.7 L988.8 595 L992.9 600.6 L997.1 605.4 L1001 609.4 L1005 612.6 L1009 615.1 L1014 616.8 L1018 617.9 L1022 618.6 L1026 618.9 L1030 619 L1034 619 L1038 619 L1043 619 L1047 619 L1051 619 L1055 619' stroke='steelblue' stroke-width='1' fill='none'></path>
<path d='M0 619 L4.137 618.8 L8.275 618.4 L12.41 617.6 L16.55 616.5 L20.69 615.1 L24.82 613.4 L28.96 611.4 L33.1 609 L37.24 606.3 L41.37 603.3 L45.51 600 L49.65 596.3 L53.78 592.3 L57.92 588 L62.06 583.3 L66.2 578.3 L70.33 573 L74.47 567.3 L78.61 561.3 L82.75 555 L86.88 548.3 L91.02 541.3 L95.16 534 L99.29 526.4 L103.4 518.4 L107.6 510.2 L111.7 501.8 L115.8 493 L120 484.1 L124.1 474.9 L128.3 465.5 L132.4 455.9 L136.5 446.1 L140.7 436.2 L144.8 426.1 L148.9 416 L153.1 405.7 L157.2 395.4 L161.4 385.1 L165.5 374.7 L169.6 364.4 L173.8 354.1 L177.9 343.9 L182 333.7 L186.2 323.6 L190.3 313.7 L194.5 303.9 L198.6 294.2 L202.7 284.7 L206.9 275.4 L211 266.3 L215.1 257.4 L219.3 248.7 L223.4 240.3 L227.5 232.1 L231.7 224.1 L235.8 216.3 L240 208.8 L244.1 201.6 L248.2 194.6 L252.4 187.8 L256.5 181.2 L260.6 174.9 L264.8 168.8 L268.9 162.9 L273.1 157.2 L277.2 151.8 L281.3 146.5 L285.5 141.3 L289.6 136.4 L293.7 131.5 L297.9 126.9 L302 122.3 L306.2 117.9 L310.3 113.6 L314.4 109.4 L318.6 105.3 L322.7 101.3 L326.8 97.36 L331 93.5 L335.1 89.72 L339.3 85.99 L343.4 82.33 L347.5 78.73 L351.7 75.19 L355.8 71.71 L359.9 68.28 L364.1 64.91 L368.2 61.6 L372.4 58.35 L376.5 55.16 L380.6 52.04 L384.8 48.99 L388.9 46.01 L393 43.11 L397.2 40.29 L401.3 37.55 L405.5 34.91 L409.6 32.35 L413.7 29.89 L417.9 27.52 L422 25.25 L426.1 23.09 L430.3 21.03 L434.4 19.08 L438.5 17.23 L442.7 15.48 L446.8 13.85 L451 12.32 L455.1 10.89 L459.2 9.56 L463.4 8.333 L467.5 7.203 L471.6 6.169 L475.8 5.226 L479.9 4.373 L484.1 3.605 L488.2 2.919 L492.3 2.313 L496.5 1.784 L500.6 1.327 L504.7 0.9415 L508.9 0.624 L513 0.3726 L517.2 0.1856 L521.3 0.06172 L525.4 0 L529.6 0 L533.7 0.06172 L537.8 0.1856 L542 0.3726 L546.1 0.624 L550.3 0.9415 L554.4 1.327 L558.5 1.784 L562.7 2.313 L566.8 2.919 L570.9 3.605 L575.1 4.373 L579.2 5.226 L583.4 6.169 L587.5 7.203 L591.6 8.333 L595.8 9.56 L599.9 10.89 L604 12.32 L608.2 13.85 L612.3 15.48 L616.5 17.23 L620.6 19.08 L624.7 21.03 L628.9 23.09 L633 25.25 L637.1 27.52 L641.3 29.89 L645.4 32.35 L649.5 34.91 L653.7 37.55 L657.8 40.29 L662 43.11 L666.1 46.01 L670.2 48.99 L674.4 52.04 L678.5 55.16 L682.6 58.35 L686.8 61.6 L690.9 64.91 L695.1 68.28 L699.2 71.71 L703.3 75.19 L707.5 78.73 L711.6 82.33 L715.7 85.99 L719.9 89.72 L724 93.5 L728.2 97.36 L732.3 101.3 L736.4 105.3 L740.6 109.4 L744.7 113.6 L748.8 117.9 L753 122.3 L757.1 126.9 L761.3 131.5 L765.4 136.4 L769.5 141.3 L773.7 146.5 L777.8 151.8 L781.9 157.2 L786.1 162.9 L790.2 168.8 L794.4 174.9 L798.5 181.2 L802.6 187.8 L806.8 194.6 L810.9 201.6 L815 208.8 L819.2 216.3 L823.3 224.1 L827.5 232.1 L831.6 240.3 L835.7 248.7 L839.9 257.4 L844 266.3 L848.1 275.4 L852.3 284.7 L856.4 294.2 L860.5 303.9 L864.7 313.7 L868.8 323.6 L873 333.7 L877.1 343.9 L881.2 354.1 L885.4 364.4 L889.5 374.7 L893.6 385.1 L897.8 395.4 L901.9 405.7 L906.1 416 L910.2 426.1 L914.3 436.2 L918.5 446.1 L922.6 455.9 L926.7 465.5 L930.9 474.9 L935 484.1 L939.2 493 L943.3 501.8 L947.4 510.2 L951.6 518.4 L955.7 526.4 L959.8 534 L964 541.3 L968.1 548.3 L972.3 555 L976.4 561.3 L980.5 567.3 L984.7 573 L988.8 578.3 L992.9 583.3 L997.1 588 L1001 592.3 L1005 596.3 L1009 600 L1014 603.3 L1018 606.3 L1022 609 L1026 611.4 L1030 613.4 L1034 615.1 L1038 616.5 L1043 617.6 L1047 618.4 L1051 618.8 L1055 619' stroke='orange' stroke-width='1' fill='none'></path>
</g>
</svg>
Side by Side

After

Width:  |  Height:  |  Size: 11 KiB

93
include/boost/math/interpolators/detail/cardinal_trigonometric_detail.hpp
View File

@@ -21,17 +21,17 @@ class cardinal_trigonometric_detail {
public:
cardinal_trigonometric_detail(const Real* data, size_t length, Real t0, Real h)
{
m_data = data;
m_length = length;
m_t0 = t0;
m_h = h;
throw std::domain_error("Not implemented.");
m_data = data;
m_length = length;
m_t0 = t0;
m_h = h;
throw std::domain_error("Not implemented.");
}
private:
size_t m_length;
Real m_t0;
Real m_h;
Real* m_data;
size_t m_length;
Real m_t0;
Real m_h;
Real* m_data;
};
template<>
@@ -41,11 +41,11 @@ public:
{
if (length == 0)
{
throw std::logic_error("At least one sample is required.");
throw std::logic_error("At least one sample is required.");
}
if (h <= 0)
{
throw std::logic_error("The step size must be > 0");
throw std::logic_error("The step size must be > 0");
}
// The period sadly must be stored, since the complex vector has length that cannot be used to recover the period:
m_T = m_h*length;
@@ -66,8 +66,15 @@ public:
float denom = length;
for (size_t k = 0; k < m_complex_vector_size; ++k)
{
m_gamma[k][0] /= denom;
m_gamma[k][1] /= denom;
m_gamma[k][0] /= denom;
m_gamma[k][1] /= denom;
}
if (length % 2 == 0)
{
m_gamma[m_complex_vector_size -1][0] /= 2;
// numerically, m_gamma[m_complex_vector_size -1][1] should be zero . . .
// I believe, but need to check, that FFTW guarantees that it is identically zero.
}
}
@@ -79,28 +86,28 @@ public:
float operator()(float t) const
{
using std::sin;
using std::cos;
using boost::math::constants::two_pi;
using std::exp;
float s = m_gamma[0][0];
float x = two_pi<float>()*(t - m_t0)/m_T;
fftwf_complex z;
// boost::math::cos_pi with a redefinition of x? Not now . . .
z[0] = cos(x);
z[1] = sin(x);
fftwf_complex b{0, 0};
// u = b*z
fftw_complex u;
for (size_t k = m_complex_vector_size - 1; k >= 1; --k) {
u[0] = b[0]*z[0] - b[1]*z[1];
u[1] = b[0]*z[1] + b[1]*z[0];
b[0] = m_gamma[k][0] + u[0];
b[1] = m_gamma[k][1] + u[1];
}
using std::sin;
using std::cos;
using boost::math::constants::two_pi;
using std::exp;
float s = m_gamma[0][0];
float x = two_pi<float>()*(t - m_t0)/m_T;
fftwf_complex z;
// boost::math::cos_pi with a redefinition of x? Not now . . .
z[0] = cos(x);
z[1] = sin(x);
fftwf_complex b{0, 0};
// u = b*z
fftw_complex u;
for (size_t k = m_complex_vector_size - 1; k >= 1; --k) {
u[0] = b[0]*z[0] - b[1]*z[1];
u[1] = b[0]*z[1] + b[1]*z[0];
b[0] = m_gamma[k][0] + u[0];
b[1] = m_gamma[k][1] + u[1];
}
s += 2*(b[0]*z[0] - b[1]*z[1]);
return s;
s += 2*(b[0]*z[0] - b[1]*z[1]);
return s;
}
float period() const
@@ -110,7 +117,7 @@ public:
float integrate() const
{
return m_T*m_gamma[0][0];
return m_T*m_gamma[0][0];
}
float squared_l2() const
@@ -126,6 +133,7 @@ public:
return s*m_T;
}
~cardinal_trigonometric_detail()
{
if (m_gamma)
@@ -176,6 +184,11 @@ public:
m_gamma[k][0] /= denom;
m_gamma[k][1] /= denom;
}
if (length % 2 == 0)
{
m_gamma[m_complex_vector_size -1][0] /= 2;
}
}
cardinal_trigonometric_detail(const cardinal_trigonometric_detail& old) = delete;
@@ -281,6 +294,10 @@ public:
m_gamma[k][0] /= denom;
m_gamma[k][1] /= denom;
}
if (length % 2 == 0) {
m_gamma[m_complex_vector_size -1][0] /= 2;
}
}
cardinal_trigonometric_detail(const cardinal_trigonometric_detail& old) = delete;
@@ -385,6 +402,10 @@ public:
m_gamma[k][0] /= denom;
m_gamma[k][1] /= denom;
}
if (length % 2 == 0)
{
m_gamma[m_complex_vector_size -1][0] /= 2;
}
}
cardinal_trigonometric_detail(const cardinal_trigonometric_detail& old) = delete;
@@ -433,7 +454,7 @@ public:
__float128 s = 0;
for (size_t i = m_complex_vector_size - 1; i >= 1; --i)
{
s += (m_gamma[i][0]*m_gamma[i][0] + m_gamma[i][1]*m_gamma[i][1]);
s += (m_gamma[i][0]*m_gamma[i][0] + m_gamma[i][1]*m_gamma[i][1]);
}
s *= 2;
s += m_gamma[0][0]*m_gamma[0][0];

31
test/cardinal_trigonometric_test.cpp
View File

@@ -36,6 +36,34 @@ void test_constant()
}
}
template<class Real>
void test_interpolation_condition()
{
std::mt19937 gen(1234);
std::uniform_real_distribution<Real> dis(1, 10);
for(size_t n = 1; n < 20; ++n) {
Real t0 = dis(gen);
Real h = dis(gen);
std::vector<Real> v(n);
for (size_t i = 0; i < n; ++i) {
v[i] = dis(gen);
}
auto ct = cardinal_trigonometric<decltype(v)>(v, t0, h);
for (size_t i = 0; i < n; ++i) {
Real arg = t0 + i*h;
Real expected = v[i];
Real computed = ct(arg);
if(!CHECK_ULP_CLOSE(expected, computed, 5*n))
{
std::cerr << " Samples: " << n << "\n";
}
}
}
}
#ifdef BOOST_HAS_FLOAT128
void test_constant_q()
{
@@ -137,6 +165,7 @@ int main()
test_constant<float>();
test_sampled_sine<float>();
test_bump<float>();
test_interpolation_condition<float>();
#endif
@@ -144,12 +173,14 @@ int main()
test_constant<double>();
test_sampled_sine<double>();
test_bump<double>();
test_interpolation_condition<double>();
#endif
#ifdef TEST3
test_constant<long double>();
test_sampled_sine<long double>();
test_bump<long double>();
test_interpolation_condition<long double>();
#endif
#ifdef TEST4