1
$\begingroup$

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: Actual 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.

$\endgroup$

closed as off-topic by Jason R, A_A, MBaz, Gilles, Laurent Duval Sep 24 '16 at 12:10

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "General programming questions are off-topic here, but can be asked on Stack Overflow." – Jason R, A_A, MBaz, Gilles, Laurent Duval
If this question can be reworded to fit the rules in the help center, please edit the question.

  • $\begingroup$ Can I please ask you to rephrase the question for clarity? Also, instead of focusing it so much about the coding bit, do you think you could isolate the parts that refer to the methodology and ask about them? Have you tried the standard approach first? Also, please note that even when things seem to "work", different packages will give slightly different numbers due to differences in data types, rounding off errors, integration methods and other factors that is impossible to control.Do you want to design a filter or work on the algorithm itself? $\endgroup$ – A_A Sep 21 '16 at 7:23
  • $\begingroup$ @A_A I have found the solution to my problem after investigating the frontend design for Octave & C++ interface in the remez.cc code. I will rephrase my question and answer it myself. This is somewhat a programming question where the design goal is to use existing library for filter generation so I am confused if I should post in DSP community or stackoverflow programming community. Sorry for the confusion I cause. $\endgroup$ – yc2986 Sep 21 '16 at 16:34
  • $\begingroup$ No worries, it does look like a better fit to SO. $\endgroup$ – A_A Sep 21 '16 at 16:47
  • $\begingroup$ @A_A I have to move a link down here because I lack the reputation to place 2+ link in the post. Could you please help me edit the question? Thanks. Octave example $\endgroup$ – yc2986 Sep 21 '16 at 17:07
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Thanks for the reply @SluethEye. I checked my remez.cc. It is downloaded from the latest sourceforge depot. The link to remez.cc you provide has the same documentation bug as mine. And I have tried some really easy test case like high pass filter. But the remez function still return weird result even it converge. I will update my case and output in the original post. $\endgroup$ – yc2986 Sep 21 '16 at 2:05
1
$\begingroup$

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! frequency response

$\endgroup$

Not the answer you're looking for? Browse other questions tagged or ask your own question.