OptimalFir::remez Function

OptimalFir Namespace

int remez(TVec *h, const DewArray<double> &bands, const DewArray<double> &gains, const DewArray<double> &Weights, TRemezType FilterType, double &err, double FS = 2, bool ConstantRipple = false);

Design an optimal equiripple FIR filter with Parks-McClellan algorithm.

#NameTypeDescription
1hTVec *An array of length h.Length on entry and contains the filter impulse response on exit.
2bandsconst DewArray<double> &Defines the frequency bands.
3gainsconst DewArray<double> &Defines, if the band is a stopband or a passband.
4Weightsconst DewArray<double> &Array contains the ratios between required ripples for different bands.
5FilterTypeTRemezTypeChoose between bandpass, hilbert, differentator and integrator.
6errdouble &Contains the maximum ripple error upon return.
7FS = 2doubleSpecifies the sampling frequency.
8ConstantRipple = falseboolIf True, the ripple will be constant and not weighted to give constant percentage error in case of the following filter types: rmtDifferentiator, rmtIntegrator, rmtDoubleDifferentiator, rmtDoubleIntegrator
Remarks:

Designs an equiripple (optimal) FIR filter. The routine will not always converge. Parameters have to be specified, for which the filter exists. Most common causes for trouble are:

  • too wide transition bands
  • transition bands not of equal width.
  • too strong attenuation (of the stopband) or to small ripple (of the passband) specified.
  • wrong filter length (odd/even).

The length of the filter can be estimated with the RemezLength routine. RemezLength routine will also properly adjust the error weights. An extensive explanation of the algorithm can be found in [1] Ch. 7.6, p. 462.

  • The length of the filter is defined as: n = Order-1. (Order is the order of the polynomial and n is the number of coefficients.)
  • highpass or a bandstop filter or any filter with the passband at FS/2 has to have odd length (even order). RemezLength routine automatically adjusts filter length.
  • a filter with negative symmetry also shifts the phase by 90 degrees. The routine automatically assumes negative symmetry, if the FilterType is different from rmtBandpass.
  • Required stopband attenuation is usually specified in dB. To obtain the required ripple, the following formula can be used: Ripple = Exp10(Att/-20);
  • Passband ripple can also be specified in dB. To obtain the required linear ripple, the following formula can be used: Ripple = (1-Exp10(PassRippldB/-20))/2
  • The passband ripple "rp" specified for the passband is a +/-rp specification. Total ripple is 2*rp.
  • the amount of ripple is a function of filter length. (longer filters give less ripple).
  • if different bands have different required ripple, this can be addressed by adjusting the error weights array.

The Fortran source code can be found in [2] p. 198.

References:

[[1] Discrete-time signal processing, Oppenheim and Schafer, Prentice-Hall, 1989.]([1] Discrete-time signal processing, Oppenheim and Schafer, Prentice-Hall, 1989.)

[[2] Theory and application of digital signal processing, Lawrence R. Rabiner and Bernard Gold. Prentice-Hall, 1975.]([2] Theory and application of digital signal processing, Lawrence R. Rabiner and Bernard Gold. Prentice-Hall, 1975.)

Examples
#include "MtxExpr.h"
#include "Units.OptimalFIR.h"
#include "Units.SignalUtils.h"

using namespace Dew;
using namespace Dew::Math;
using namespace Dew::Signal;
using namespace Dew::Signal::Units::OptimalFir;
using namespace Dew::Signal::Units::SignalUtils;

void Example()
{
    Vector z(0);
    Vector Response(0);
    Vector Weights(0);

    int n;
    double err;
    double FS = 2; //Sampling frequency

    //Design a lowpass
    n = RemezLength(DewArray<double>{0, 1500, 2000, 4000},
        DewArray<double>{1, 0}, DewArray<double>{0.001, 0.001}, Weights, 8000);

    z.Size(n);
    remez(z, DewArray<double>{0, 1500, 2000, 4000},
        DewArray<double>{1, 0}, (DewArray<double>) Weights, TRemezType::rmtBandPass, err, 8000, false);

    FrequencyResponse(z, nullptr, Response, 8, false, TSignalWindowType::wtRectangular, 0);

    //Design a highpass
    n = RemezLength(DewArray<double>{0, 1500, 2000, 4000},
        DewArray<double>{0, 1}, DewArray<double>{0.001, 0.001}, Weights, 8000);
    z.Size(n);
    remez(z, DewArray<double>{0, 1500, 2000, 4000},
        DewArray<double>{0, 1}, (DewArray<double>) Weights, TRemezType::rmtBandPass, err, 8000, false);

    FrequencyResponse(z, nullptr, Response, 8, false, TSignalWindowType::wtRectangular, 0);

    //Design a bandpass (Sampling frequency = 2)
    n = RemezLength(DewArray<double>{0, 0.15, 0.25, 0.45, 0.55, 1},
        DewArray<double>{0, 1, 0}, DewArray<double>{0.001, 0.001, 0.001}, Weights, FS);
    z.Size(n);
    remez(z, DewArray<double>{0, 0.15, 0.25, 0.45, 0.55, 1},
        DewArray<double>{0, 1, 0}, (DewArray<double>) Weights, TRemezType::rmtBandPass, err, FS, false);

    FrequencyResponse(z, nullptr, Response, 8, false, TSignalWindowType::wtRectangular, 0);

    //Design a bandstop (Sampling frequency = 2)
    n = RemezLength(DewArray<double>{0, 0.15, 0.25, 0.45, 0.55, 1},
        DewArray<double>{1, 0, 1}, DewArray<double>{0.001, 0.001, 0.001}, Weights, FS);
    z.Size(n);
    remez(z, DewArray<double>{0, 0.15, 0.25, 0.45, 0.55, 1},
        DewArray<double>{1, 0, 1}, (DewArray<double>) Weights, TRemezType::rmtBandPass, err, FS, false);

    FrequencyResponse(z, nullptr, Response, 8, false, TSignalWindowType::wtRectangular, 0);

    //Design a multiband (Sampling frequency = 2)
    n = RemezLength(DewArray<double>{0.00, 0.10, 0.20, 0.30, 0.40, 0.45, 0.55, 0.60, 0.70, 1.00},
        DewArray<double>{1, 0, 1, 0, 1}, DewArray<double>{0.001, 0.001, 0.001, 0.001, 0.001}, Weights, FS);
    z.Size(n);
    remez(z, DewArray<double>{0.00, 0.10, 0.20, 0.30, 0.40, 0.45, 0.55, 0.60, 0.70, 1.00},
        DewArray<double>{1, 0, 1, 0, 1}, (DewArray<double>) Weights, TRemezType::rmtBandPass, err, FS, false);

    FrequencyResponse(z, nullptr, Response, 8, false, TSignalWindowType::wtRectangular, 0);

    //  Design a hilbert (type III) transformer  (Sampling frequency = 2)
    n = RemezLength(DewArray<double>{0, 0, 0.1, 0.9, 1, 1},
        DewArray<double>{0, 1, 0}, DewArray<double>{0.01, 0.01, 0.01}, Weights, FS);
    if (n % 2 == 0) n++; //odd length type III filter
    z.Size(n);
    remez(z, DewArray<double>{0.1, 0.9},
        DewArray<double>{1}, DewArray<double>{1}, TRemezType::rmtHilbert, err, FS, false);

    FrequencyResponse(z, nullptr, Response, 8, false, TSignalWindowType::wtRectangular, 0);

    //Design a hilbert (type IV) transformer  (Sampling frequency = 2)
    n = RemezLength(DewArray<double>{0, 0, 0.1, 0.9, 1, 1},
        DewArray<double>{0, 1, 0}, DewArray<double>{0.01, 0.01, 0.01}, Weights, FS);
    if (n % 2 != 0) n++; //even length type IV filter
    z.Size(n);
    remez(z, DewArray<double>{0.1, 1},
        DewArray<double>{1}, DewArray<double>{1}, TRemezType::rmtHilbert, err, 2, false);

    FrequencyResponse(z, nullptr, Response, 8, false, TSignalWindowType::wtRectangular, 0);

    //Design a differentiator (type III) (Sampling frequency = 2)
    n = RemezLength(DewArray<double>{0, 0, 0.1, 0.9, 1, 1},
        DewArray<double>{0, 1, 0}, DewArray<double>{0.01, 0.01, 0.01}, Weights, FS);
    if (n % 2 == 0) n++; //odd length type III filter
    z.Size(n);
    remez(z, DewArray<double>{0.1, 0.9},
        DewArray<double>{1}, DewArray<double>{1}, TRemezType::rmtDifferentiator, err, 2, false);

    FrequencyResponse(z, nullptr, Response, 8, false, TSignalWindowType::wtRectangular, 0);

    //Design a differentiator (type IV) (Sampling frequency = 2)
    n = RemezLength(DewArray<double>{0, 0, 0.1, 0.9, 1, 1},
        DewArray<double>{0, 1, 0}, DewArray<double>{0.01, 0.01, 0.01}, Weights, FS);
    if (n % 2 != 0) n++; //even length type IV filter
    z.Size(n);
    remez(z, DewArray<double>{0.1, 1},
        DewArray<double>{1}, DewArray<double>{1}, TRemezType::rmtDifferentiator, err, FS, false);

    FrequencyResponse(z, nullptr, Response, 8, false, TSignalWindowType::wtRectangular, 0);

    //Design a double differentiator (Sampling frequency = 2)
    n = RemezLength(DewArray<double>{0, 0, 0.05, 0.95, 1, 1},
        DewArray<double>{0, 1, 0}, DewArray<double>{0.01, 0.01, 0.01}, Weights, FS);
    if (n %2 != 0) n++; //even length type IV filter
    z.Size(n);
    remez(z, DewArray<double>{0.05, 1},
        DewArray<double>{1}, DewArray<double>{1}, TRemezType::rmtDoubleDifferentiator, err, FS, false);

    FrequencyResponse(z, nullptr, Response, 8, false, TSignalWindowType::wtRectangular, 0);

    //Design an integrator (Sampling frequency = 2)
    n = RemezLength(DewArray<double>{0, 0, 0.05, 0.95, 1, 1},
        DewArray<double>{0, 1, 0}, DewArray<double>{0.01, 0.01, 0.01}, Weights, FS);
    if (n %2 != 0) n++; //even length type IV filter
    z.Size(n);
    remez(z, DewArray<double>{0.05, 1},
        DewArray<double>{1}, DewArray<double>{1}, TRemezType::rmtIntegrator, err, FS, false);

    FrequencyResponse(z, nullptr, Response, 8, false, TSignalWindowType::wtRectangular, 0);

    //Design a double integrator  (Sampling frequency = 2)
    n = RemezLength(DewArray<double>{0, 0, 0.05, 0.95, 1, 1},
        DewArray<double>{0, 1, 0}, DewArray<double>{0.01, 0.01, 0.01}, Weights, FS);
    if (n %2 != 0) n++; //even length type IV filter
    z.Size(n);
    remez(z, DewArray<double>{0.05, 1},
        DewArray<double>{1}, DewArray<double>{1}, TRemezType::rmtDoubleIntegrator, err, FS, false);

    FrequencyResponse(z, nullptr, Response, 8, false, TSignalWindowType::wtRectangular, 0);
}
See Also: KaiserImpulse, SavGolayImpulse, OptimalFir::RemezImpulse, OptimalFir::RemezLength
Declared in Dew::Signal::Units::OptimalFir · Dew.Signal/Units.OptimalFIR.h · Cross-compiler