Currently, pass/stopbands are specified using std::vectors for frequency, amplitude, and weight. In principle, this allows arbitrary FIRs to be constructed using a very high number of subbands. However, it doesn't work in practice (I can substantiate "doesn't work" with a little effort: from memory, however, it's numerically unstable at subband boundaries.)

Specifying pass/stopbands as a lambda function accepting frequency and returning an (amplitude, weight) tuple would allow an arbitrary amplitude specifier and may actually reduce the number of special cases in the algorithm.

Question for @sfilip: is this an interesting refactor? If it were submitted as a pull request, what would your gut reaction be?

This test case has two specifications both with identical pass and stop bands ... the difference being the first specification has its pass band partitioned.

With the test case built like so:

  clang++ -I/usr/local/include -O2 -march=amdfam10
    -o firpm-tc-10ka firpm-tc-10ka.cxx
    -L/usr/local/lib -lfirpm -lmpfr -lgmp

Running it results in:

  START Parks-McClellan with uniform initialization
  WARNING: The exchange algorithm did not converge.
  TRIGGER: Not enough alternating extrema!
  POSSIBLE CAUSE: Nmax too small
  WARNING: The exchange algorithm did not converge.
  TRIGGER: numerical instability
  POSSIBLE CAUSES: poor starting reference and/or a too small value for Nmax.
  Iteration count = NC
  FINISH Parks-McClellan with uniform initialization

  START Parks-McClellan with uniform initialization
  Final delta     = 1.4049544582785901791e-05
  Iteration count = 55
  FINISH Parks-McClellan with uniform initialization

What's interesting is the two filter specifications are very similar ... firpm fails for the first specification and succeeds for the second.

Perhaps this is just the nature of things ... anyway I pass this along in case comparing / contrasting the two behaviors suggests something.
firpm-tc-10ka.cxx.txt

With the test case built like so:

  clang++ -I/usr/local/include -O2 -march=amdfam10
    -o firpm-tc-10kb firpm-tc-10kb.cxx
    -L/usr/local/lib -lfirpm -lmpfr -lgmp

Running it results in:

  START Parks-McClellan with reference scaling
  WARNING: The exchange algorithm did not converge.
  TRIGGER: Not enough alternating extrema!
  POSSIBLE CAUSE: Nmax too small
  WARNING: The exchange algorithm did not converge.
  TRIGGER: numerical instability
  POSSIBLE CAUSES: poor starting reference and/or a too small value for Nmax.
  Segmentation fault (core dumped)

which is to say it crashes in addition to not converging.

firpm-tc-10kb.cxx.txt

My understanding is that ideally firpmRS should succeed for any
specification where firpm succeeded.

Attach is a test case based on a CIC compensation filter where
firpm succeeds and firpmRS fails. The code is built like so:

clang++ -I/usr/local/include -O2 -march=amdfam10 \
  -o firpm_tc firpm_tc.cxx -L/usr/local/lib -lfirpm

and running it results in:

  START Parks-McClellan with uniform initialization
  Final delta = 0.04104066019515401947
  Iteration count = 12
  FINISH Parks-McClellan with uniform initialization
  START Parks-McClellan with reference scaling
  The exchange algorithm did not converge.
  TRIGGER: Not enough alternating extrema!
  POSSIBLE CAUSE: Nmax too small
  Failed to do reference scaling since x is empty!
[firpm_tc.cxx.txt](https://github.com/sfilip/firpm/files/3822819/firpm_tc.cxx.txt)

This is using firpm_ld.

It's a bit inconvenient to have firpm terminate the calling program upon failure. Some of my samples converge successfully using firpmRS and not firpmAFP, others converge successfully using firpmAFP and not firpmRS.

The attached patch changes the firpm library so that it internally throws exceptions instead of exiting. This allows the caller to see that a failure occurred (indicated by the value of output.q), giving the caller a chance to decide how it should handle the failure (i.e. perhaps by retrying with a different initialization approach).

The supplied change is just a quick minimal approach ... it might be useful to provide the failure code to the caller. Another useful change might be to accept an error stream parameter instead of assuming that it's okay to write to std::cerr.
firpm-sfilip-PatchJLW02.txt

In trying a particular filter specification with your firpm fix-refactor
I was able to trigger various crashes. The platform I'm using is FreeBSD 11.3
with the OS supplied compiler (clang version 8.0.1).

Issues fixed by the enclosed patch:

1. There are lambdas of the form:

   fbands[i].weight = [i, w](space_t, double x) ...
   fbands[i].amplitude = [i, a, fbands](space_t space, double x) ...

   Unfortunately capturing arrays by value causes the process to consume
   excessive memory before finally dying.

2. uniform can be called with an omega array that's smaller than the
   band array. This results in avgDistance being negative so latter
   extremas calculations wrap to a large number causing a crash.

3. If firpmRS calls exchange which fails, then output.x will be empty.
   This causes referenceScaling to crash.

Notes:

a) These same changes probably should be propagated to both firpm_ld
and firpm_mp.

-- John
PatchJLW01-fr.txt

Interesting work. Thank you for sharing it.

I was trying out the double precision variant of firpm and the following causes a seg fault on Linux with gcc.

PMOutput result = firpm(20, {0.0, 0.4, 0.6, 0.64, 0.69, 0.74, 0.79, 0.83, 0.88, 1.0}, {1.0, 1.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 1.0, 1.0}, {1.0, 1.0, 1.0, 1.0, 1.0});

And this code causes runaway memory usage. (Make sure to ctrl-c before you run out of memory)

PMOutput result = firpm(12, {0.0, 0.4, 0.6, 0.6089285714285714, 0.6178571428571429, 0.6267857142857143, 0.6357142857142857, 0.6446428571428571, 0.6535714285714286, 0.6625, 0.6714285714285714, 0.6803571428571428, 0.6892857142857143, 0.6982142857142857, 0.7071428571428571, 0.7160714285714285, 0.725, 0.7339285714285714, 0.7428571428571429, 0.7517857142857143, 0.7607142857142857, 0.7696428571428571, 0.7785714285714285, 0.7874999999999999, 0.7964285714285714, 0.8053571428571428, 0.8142857142857143, 0.8232142857142857, 0.8321428571428571, 0.8410714285714285, 0.85, 0.8589285714285714, 0.8678571428571429, 0.8767857142857143, 0.8857142857142857, 0.8946428571428571, 0.9035714285714285, 0.9124999999999999, 0.9214285714285714, 0.9303571428571428, 0.9392857142857143, 0.9482142857142857, 0.9571428571428571, 0.9660714285714285, 0.975, 1.0}, {1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, {87.35985578898647, 1000.0, 1029.7619047619048, 1059.5238095238094, 1089.2857142857144, 1119.047619047619, 1148.8095238095239, 1178.5714285714284, 1208.3333333333333, 1238.0952380952383, 1267.857142857143, 1297.6190476190475, 1327.3809523809523, 1357.142857142857, 1386.904761904762, 1416.6666666666667, 1446.4285714285716, 1476.1904761904761, 1505.952380952381, 1535.7142857142858, 1565.4761904761906, 1595.2380952380952, 1625.0});

Hi,
I have some experience with C++ and Python, but never written C/C++ extensions for Python. Do you think it could be possible to create Python package without much changes to current codebase? Or perhaps the best solution is to build the library and use ctypes for interfacing?

With the test case built like so:

clang++ -I/usr/local/include -O2 -march=amdfam10
-o firpm-tc-10k firpm-tc-10k.cxx
-L/usr/local/lib -lfirpm -lmpfr -lgmp

Running it results in:

./firpm-tc-10k
START Parks-McClellan with AFP
Segmentation fault (core dumped)

This is simplified test case using firpmAFP ... it's from
a larger version which crashes the regular firpm.
firpm-tc-10k.cxx.txt

Compiling firpm-tc-50ka like so:

  clang++ -I/usr/local/include -O2 -march=amdfam10 -DHAVE_MPFR \
    -o firpm-tc-50ka firpm-tc-50ka.cxx -L/usr/local/lib -lfirpm -lmpfr -lgmp

and then running it results in:

  Filter degree 482
  START Parks-McClellan with uniform initialization
  Final delta     = 5.88012e-05
  Iteration count = 32
  FINISH Parks-McClellan with uniform initialization
  START Parks-McClellan with reference scaling
  Final delta     = 5.90794e-05
  Iteration count = 9
  FINISH Parks-McClellan with reference scaling
  START Parks-McClellan with AFP
  Final delta     = 5.94789e-05
  Iteration count = 8
  FINISH Parks-McClellan with AFP
  Iteration reduction for final filter  RS: 0.71875
  Iteration reduction for final filter AFP: 0.75

however looking at the resulting files:

  firpm-tc-50ka-UNI-coeff.txt
  firpm-tc-50ka-RS-coeff.txt
  firpm-tc-50ka-AFP-coeff.txt

shows all the RS coefficients are NaN.

firpm-tc-50ka.cxx.txt