vinniefalco / DSPFilters

A Collection of Useful C++ Classes for Digital Signal Processing
1.87k stars 380 forks source link

Assert fails in debug code for band stop filter #31

Open dojeda opened 8 years ago

dojeda commented 8 years ago

An old bug of a colleague led me to find this code in PoleFilter.cpp (BandStopTransform::BandStopTransform method).

The issue here is that the third assert fails because the difference between zc.first and z.first is very small (1e-16) but not zero. This is not an unexpected scenario; floating-point arithmetic is not exact and this happens a lot (which is why exact comparisons of floats is not recommended).

Meanwhile, the code below makes me wonder: there a conjugate trick in NDEBUG to set z, but there is no modification at all of z in DEBUG. Is this normal? It certainly does not seem correct!

Are there any details regarding this conjugate trick?


  const int numPoles = analog.getNumPoles ();
  const int pairs = numPoles / 2;
  for (int i = 0; i < pairs; ++i)
  {
    const PoleZeroPair& pair = analog[i];
    ComplexPair p  = transform (pair.poles.first);
    ComplexPair z  = transform (pair.zeros.first);

    //
    // Optimize out the calculations for conjugates for Release builds
    //
#ifdef NDEBUG
    // trick to get the conjugate
    if (z.second == z.first)
      z.second = std::conj (z.first);

#else
    // Do the full calculation to verify correctness
    ComplexPair pc = transform (analog[i].poles.second);
    ComplexPair zc = transform (analog[i].zeros.second);

    // get the conjugates into pc and zc
    if (zc.first == z.first)
      std::swap (zc.first, zc.second);

    assert (pc.first  == std::conj (p.first));
    assert (pc.second == std::conj (p.second));
    assert (zc.first  == std::conj (z.first));
    assert (zc.second == std::conj (z.second));

#endif

    digital.addPoleZeroConjugatePairs (p.first, z.first);
    digital.addPoleZeroConjugatePairs (p.second, z.second);
  }