C++ based FIR filter design using GNU Octave remez function [closed]

Question

I am designing a system where I want a FIR filter coefficient generator feature. I try the remez.cc source file provided in the Octave Signal package.

The source code for remez.cc is here. I have removed the octave interface part (line 757 - end) and the library header (line 34) and created a simple header with only remez function at line 592 then compiled the source code into a static library.

#pragma once
// remez.h
/********************
* remez
*=======
* Calculates the optimal (in the Chebyshev/minimax sense)
* FIR filter impulse response given a set of band edges,
* the desired response on those bands, and the weight given to
* the error in those bands.
*
* INPUT:
* ------
* int     numtaps     - Number of filter coefficients
* int     numband     - Number of bands in filter specification
* double  bands[]     - User-specified band edges [2 * numband]
* double  des[]       - User-specified band responses [numband]
* double  weight[]    - User-specified error weights [numband]
* int     type        - Type of filter
*
* OUTPUT:
* -------
* double h[]      - Impulse response of final filter [numtaps]
* returns         - true on success, false on failure to converge
********************/
int remez(double h[], int numtaps,
          int numband, const double bands[],
          const double des[], const double weight[],
          int type, int griddensity);

Then I followed an Octave example for filter design and implemented a C++ main function for experimenting how to call the remez function. Here is my main function.

#include <remez.h>
#include <cmath>
#include <vector>
#include <iostream>

using namespace std;

int main()
{
    int numtaps = 52;
    vector<double> h(numtaps);
    vector<double> bands = { 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 };
    int numbands = bands.size() / 2;
    vector<double> des = { 1, 1, 0, 0, 0.5, 0.5, 0, 0, 1, 1 };  // designed response for each band

    double ripple1 = 1 - pow(10, -0.3 / 20);
    double att2 = pow(10, -60 / 20);
    double ripple3 = (1 - pow(10, -0.2 / 20)) * 0.5;
    double att4 = pow(10, -70 / 20);
    double ripple5 = 1 - pow(10, -0.1 / 20);
    vector<double> weight = { 1 / ripple1, 1 / att2, 1 / ripple3, 1 / att4, 1 / ripple5 };

    int err = remez(h.data(), numtaps, numbands, bands.data(), des.data(), weight.data(), 1, 16);
    for (auto &c : h)
        cout << c << ", ";
    cout << endl;
    cout << "error code: " << err << endl;
    return 0;
}

This function call produce some weird coefficients. All values inside h vector are 10e144 something. It is clearly an overflow so I debugged the remez.cc and traced what cause this.

The remez function documentation specified that des array should be of size numband. The des array in remez function is passed to CreateDenseGrid function on line 86 for grid generation. The documentation for CreateDenseGrid on line 73 indicate that its parameter des should be size of 2*numband. Clearly since my des array is only of size numband some out of bound memory access happened on line 110.

I changed the des vector to be {1, 1, 0, 0, 0.5, 0.5, 0, 0, 1, 1} and tries again. The coefficients looks nicer, but when I plot the frequency response in Octave it is not anything close to my intended filter design.

I get this weird frequency response:

Anyone could help me with this issue? Is anything wrong with my remez function call? Any suggestion is appreciated here. Thank you very much.

By the way. I am using vs2015, windows7 64-bits, Debug X86 build for compiler settings.

Answer 1

After some debugging and remez.cc source code investigation, I found the correct way to use the backend C++ implementation for Octave's remez function.

Thanks to @SleuthEye the latest version of remez.cc is here.

I have made several mistakes when calling remez function in C++.

1. I lost some accuracy due to integer division when initialize the passband/stopband ripple.
2. I didn't notice that numtaps of filter should return (numtaps + 1) of filter coefficients.
3. Inside the removed Octave interface part from remez.cc, there are some data processing and valuable error code which indicate how to feed data to the C++ remez function.

Here is my updated code using remez.cc with correct input parameters and comment.

#include <remez.h>
#include <cmath>
#include <vector>
#include <iostream>

using namespace std;

int main()
{
    int filterorder = 52;
    vector<double> h(filterorder + 1);  // filter coefficients buffer should be numtaps + 1
    vector<double> bands = { 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 };
    int numbands = bands.size() / 2;
    vector<double> des = { 1, 1, 0, 0, 0.5, 0.5, 0, 0, 1, 1 };  // designed response for each band

    /** 
     * casting to double before calculation 
     */
    double ripple1 = 1 - pow(10, -0.3 / 20);
    double att2 = pow(10, static_cast<double>(-60) / 20);
    double ripple3 = (1 - pow(10, -0.2 / 20)) * 0.5;
    double att4 = pow(10, static_cast<double>(-70) / 20);
    double ripple5 = 1 - pow(10, -0.1 / 20);
    vector<double> weight = { 1 / ripple1, 1 / att2, 1 / ripple3, 1 / att4, 1 / ripple5 };

    /**
     * According to remez.cc line 810
     * the frequency band should be devided by 2
     * before feeding into C++ remez function
     */
    for (auto &b : bands)
        b /= 2;

    /**
     * According to remez.cc line 787
     * numtap = filterorder + 1
     */
    int err = remez(h.data(), filterorder + 1, numbands, bands.data(), des.data(), weight.data(), 1, 16);
    cout << "coefficients: " << endl;
    for (auto &c : h)
        cout << c << ", ";
    cout << endl;
    cout << "error code: " << err << endl;
    return 0;
}

Some hidden requirement in remez.cc remez function call can only be found in the Octave interface code. These requirements are not documented or incorrectly documented in the comment for remez function in remez.cc. After fixing all of these discrepancies I am able to get the correct frequency response!

Answer 2

There are some known bugs (such as bug #38134 and bug #48078 on savannah) with the current implementation of remez found in Octave's signal package (see remez.cc on sourceforge), not to mention that the source to link to seems to lag a little from that latest version (you are thus lacking some critical fixes that have been included to address various stability and convergence issues).

Since you appear to be reasonably comfortable with modifying sources, you may want to try out the latest remez.cc source with the submitted patch attached to bug #38134, or alternatively the newly submitted firpm replacement to remez.