Where can I get an open source implementation of the Parks-McClellan optimal FIR filter design algorithm?

Question

Background: Often times I am doing some sort of signal processing task that requires a unique filter. Usually at this point I go to MATLAB and generate a new unique filter using $\tt firpm()$. The MATLAB firpm() function implements that Parks-McClellan algorithm. Now I have a filter and I put the filter into a hardcoded array. But here's the problem I now have a hardcoded filter that only works for one scenario.

The problem: I can now solve my signal processing problem du-jour... but only for a very SPECIFIC single sample rate or SPECIFIC scenario.

The goal: I want to be able to call $\tt firpm()$ from C code or some other language and make my signal processing code more generic. I can't find an open source implementation of firpm() !

Where can I get an open source implementation of the Parks-McClellan optimal FIR filter design algorithm (aka $\tt firpm()$ in MATLAB)?

• P.S. I am aware that I can design filters differently using windowing or something else... feel free to mention those in the comments. But the point of this question is not to ask "what are other filter design techniques?" the point is to find an open source implementation of the VERY VERY useful $\tt firpm()$... or something similar.

• P.P.S. One of the goals of this question is to learn how the Parks-McClellan algorithm works by looking at the code first and then I plan on reading some background theory.

Answer 1

There is an open-source implementation of Parks-McClellan (also known as the Remez exchange algorithm) in GNU Octave, a free-software implementation of a MATLAB-like environment. The function, called "remez", is contained in the "signal" package, which is hosted at Octave-Forge. If you download the package, you'll find "remez.cc", a C++ implementation of the algorithm.

One nice thing about Octave is that it is almost code-compatible with MATLAB, so you can easily port code over to use there if you like. It's a good way to peek under the hood at implementations of algorithms that are provided in MEX form in MATLAB.

Answer 2

Here's an LGPL version of the Remez exchange algorithm. The octave code seems to be derived from it. It was linked from the wikipedia page Parks McClellan page.
The original Janovetz code might be more easily usable in your project, since it doesn't have the octave calls, but it would be wise to dig through the octave-forge svn changelog for any info about bugfixes or speedups in the remez.cc file.

Answer 3

A convenient version can be found in Python's scipy.signal.remez. Nice if using numpy/scipy.

Answer 4

Here is another source for the Parks McClellan algorithm in C. This code is different from the SciPy code mentioned above in that it has 61 of the original 69 goto statements removed (the SciPy code still has about 37 goto's). It also fixes the code in 3 places where divide by zero can occur and it has some additional code that range checks the band edge values.

http://www.iowahills.com/A7ExampleCodePage.html

Answer 5

maybe you already know this, but if you have matlab you can use the matlab coder, and create a simple function that uses the feature you want to examine. Then run it and see the created C code. I tried with with the FFT, and with QR decomposition, and while it is a bit messy, it can be understood just fine.

Answer 6

Here is a paper which does an actual Matlab version of the "core" remez algorithm. “A MATLAB based optimum multiband FIR filters design program following the original idea of the Remez multiple exchange algorithm” -2011 IEEE International Symposium on Circuits and Systems (ISCAS) - Authors (Ahsan, Saramaki)

This paper does a good job at explaining the basic algorithm. The goal of the paper was to avoid using the original Fortran code - which does not explain the algorithm very well and often just gets translated into various other languages directly.

One thing I will comment on. One of the core ideas of the algorithm is to fit a curve and then find the extremal points. Usually a Lagrange Interpolation is used to explain this idea, but Lagrange Interpolation has poor numerical properties. In the original algorithm the use Barycentric Implementation of Lagrange Interpolation - which avoids many of the associated pitfalls of Lagrangian interpolation. So if you are trying to fully understand the code, you may want to look up Barycentric Interpolation.

Answer 7

Beware of the differences between Matlab's firpm and Scipy.signal's remez. For example, these two statements are equivalent:

% Matlab
firpm(10,[.2 .8],[1 1],'Hilbert')

# Python
from scipy.signal import remez

remez(11, [0.1, 0.4], [1], type='hilbert')