DewDSPMasterNET
remez Routine
Summary
Design an optimal equiripple FIR filter with Parks-McClellan algorithm.

Unit
OptimalFir

Declaration
Function remez(h: TVec; const bands, gains, weights: array of TSample; FilterType: TRemezType; out err: TSample; FS: TSample = 2; ConstantRipple: boolean = false): integer;
 Parameter  Description 
h An array of length Length(h) on entry and containes the filter impulse response on exit. 
Bands Defines the frequency bands. Gains define, if the band is a stopband or a passband. 
Weights array contain the ratios between required ripples for different bands. 
FS specifies the sampling frequency. 
Err contains the maximum ripple error upon return. 
ConstantRipple If 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

 

Description
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:

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.

A few need to know things about FIR filters:

The Fortran source code can be found in [2] p. 198
Categories
FIR filter design routines
 See Also 
[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". 
KaiserImpulse 
SavGolayImpulse 
RemezImpulse 
RemezLength 

Example 1

A sample lowpass filter design with sampling frequency 8000Hz, transition band between at 1500Hz and 2000Hz with 60 dB attenuation (20 * Log10(0.001)) and 0.001 ripple in the passband. To try out other setups comment out the lowpass and comment in the desired filter. 

  uses MtxExpr, Math387, MtxVec, SignalUtils, MtxVecTee, MtxVecEdit, OptimalFIR;

  procedure TForm1.Button1Click(Sender: TObject);
  var z,Response,Weights: Vector;
      n: integer;
      err: TSample;
  begin
        n := RemezLength([0, 1500, 2000, 4000], [1, 0], [0.001, 0.001],Weights, 8000);
        z.Size(n);
        Remez(z,[0, 1500, 2000, 4000], [1, 0], TSampleArray(Weights), rmtBandpass,err,8000);

  //Design a highpass
  //      n := RemezLength([0, 1500, 2000, 4000], [0, 1], [0.001, 0.001],Weights, 8000);
  //      z.Size(n);
  //      Remez(z,[0, 1500, 2000, 4000], [0, 1], TSampleArray(Weights), rmtBandpass,err,8000);

  //Design a bandpass (Sampling frequency = 2)
  //      n := RemezLength([0, 0.15, 0.25, 0.45, 0.55, 1], [0, 1, 0], [0.001, 0.001, 0.001],Weights);
  //      z.Size(n);
  //      Remez(z,[0, 0.15, 0.25, 0.45, 0.55, 1], [0, 1, 0], TSampleArray(Weights), rmtBandpass, err);

  //Design a bandstop (Sampling frequency = 2)
  //      n := RemezLength([0, 0.15, 0.25, 0.45, 0.55, 1], [1, 0, 1], [0.001, 0.001, 0.001],Weights);
  //      z.Size(n);
  //      Remez(z,[0, 0.15, 0.25, 0.45, 0.55, 1], [1, 0, 1], TSampleArray(Weights), rmtBandpass, err);

  //Design a multiband (Sampling frequency = 2)
  //      n := RemezLength([0.00, 0.10, 0.20, 0.30, 0.40,
  //                        0.45, 0.55, 0.60, 0.70, 1.00], [1, 0, 1, 0, 1], [0.001, 0.001, 0.001, 0.001, 0.001],Weights);
  //      z.Size(n);
  //      Remez(z,[0.00, 0.10, 0.20, 0.30, 0.40,
  //                      0.45, 0.55, 0.60, 0.70, 1.00], [1, 0, 1, 0, 1], TSampleArray(Weights), rmtBandpass, err);

  //Design a hilbert (type III) transformer  (Sampling frequency = 2)
  //        n := RemezLength([0, 0, 0.1, 0.9, 1, 1], [0,1,0], [0.01,0.01,0.01],Weights);
  //        if not Odd(n) then Inc(n); //odd length type III filter
  //        z.Size(n);
  //        Remez(z,[0.1, 0.9], [1], [1], rmtHilbert,err);

  //Design a hilbert (type IV) transformer  (Sampling frequency = 2)
  //        n := RemezLength([0, 0, 0.1, 0.9, 1, 1], [0,1,0], [0.01,0.01,0.01],Weights);
  //        if Odd(n) then Inc(n); //even length type IV filter
  //        z.Size(n);
  //        Remez(z,[0.1, 1], [1], [1], rmtHilbert,err);

  //Design a differentiator (type III) (Sampling frequency = 2)
  //        n := RemezLength([0, 0, 0.1, 0.9, 1, 1], [0,1,0], [0.01,0.01,0.01],Weights);
  //        if not Odd(n) then Inc(n); //odd length type III filter
  //        z.Size(n);
  //        Remez(z,[0.1, 0.9], [1], [1], rmtDifferentiator,err);

  //Design a differentiator (type IV) (Sampling frequency = 2)
  //        n := RemezLength([0, 0, 0.1, 0.9, 1, 1], [0,1,0], [0.01,0.01,0.01],Weights);
  //        if Odd(n) then Inc(n); //even length type IV filter
  //        z.Size(n);
  //        Remez(z,[0.1, 1], [1], [1], rmtDifferentiator,err);

  //Design a double differentiator (Sampling frequency = 2)
  //        n := RemezLength([0, 0, 0.05, 0.95, 1, 1], [0,1,0], [0.01,0.01,0.01],Weights);
  //        if Odd(n) then Inc(n); //even length type IV filter
  //        z.Size(n);
  //        Remez(z,[0.05, 1], [1], [1], rmtDoubleDifferentiator,err);

  //Design an integrator (Sampling frequency = 2)
  //        n := RemezLength([0, 0, 0.05, 0.95, 1, 1], [0,1,0], [0.01,0.01,0.01],Weights);
  //        if Odd(n) then Inc(n); //even length type IV filter
  //        z.Size(n);
  //        Remez(z,[0.05, 1], [1], [1], rmtIntegrator,err);

  //Design a double integrator  (Sampling frequency = 2)
  //        n := RemezLength([0, 0, 0.05, 0.95, 1, 1], [0,1,0], [0.01,0.01,0.01],Weights);
  //        if Odd(n) then Inc(n); //even length type IV filter
  //        z.Size(n);
  //        Remez(z,[0.05, 1], [1], [1], rmtDoubleIntegrator,err);

        FrequencyResponse(z,nil,Response,8);
        DrawIt(Response);
  end;


  #include "MtxVecCPP.h" //MtxVecCPP.cpp must be included in the project
  #include "MtxVecEdit.hpp"
  #include "MtxVecTee.hpp"
  #include "SignalUtils.hpp"
  #include "OptimalFIR.hpp"

  void __fastcall TForm1::BitBtn1Click(TObject *Sender)
  {
      Vector z, Response, Weights;
      int n;
      double err;

      n = RemezLength(OPENARRAY(double,(0, 1500, 2000, 4000)), OPENARRAY(double,(1, 0)), OPENARRAY(double,(0.001, 0.001)),Weights, 8000);
      z->Size(n);
      remez(z,OPENARRAY(double,(0, 1500, 2000, 4000)), OPENARRAY(double,(1, 0)), Weights->PValues(0), Weights->Length-1, rmtBandPass,err,8000);

    //Design a highpass
  //		n = RemezLength(OPENARRAY(double,(0, 1500, 2000, 4000)), OPENARRAY(double,(0, 1)), OPENARRAY(double,(0.001, 0.001)),Weights, 8000);
  //        z->Size(n);
  //		remez(z,OPENARRAY(double,(0, 1500, 2000, 4000)), OPENARRAY(double,(0, 1)), Weights->PValues(0), Weights->Length-1, rmtBandPass,err,8000);

    //Design a bandpass (Sampling frequency = 2)
  //		n = RemezLength(OPENARRAY(double,(0, 0.15, 0.25, 0.45, 0.55, 1)), OPENARRAY(double,(0, 1, 0)), OPENARRAY(double,(0.001, 0.001, 0.001)),Weights);
  //        z->Size(n);
  //		remez(z,OPENARRAY(double,(0, 0.15, 0.25, 0.45, 0.55, 1)), OPENARRAY(double,(0, 1, 0)), Weights->PValues(0), Weights->Length-1, rmtBandPass, err);

    //Design a bandstop (Sampling frequency = 2)
  //		n = RemezLength(OPENARRAY(double,(0, 0.15, 0.25, 0.45, 0.55, 1)), OPENARRAY(double,(1, 0, 1)), OPENARRAY(double,(0.001, 0.001, 0.001)),Weights);
  //        z->Size(n);
  //		remez(z,OPENARRAY(double,(0, 0.15, 0.25, 0.45, 0.55, 1)), OPENARRAY(double,(1, 0, 1)), Weights->PValues(0), Weights->Length-1, rmtBandPass, err);

    //Design a multiband (Sampling frequency = 2)
  //		n = RemezLength(OPENARRAY(double,(0.00, 0.10, 0.20, 0.30, 0.40, 0.45, 0.55, 0.60, 0.70, 1.00)),
  //					          OPENARRAY(double,(1, 0, 1, 0, 1)), OPENARRAY(double,(0.001, 0.001, 0.001, 0.001, 0.001)),Weights);
  //		z->Size(n);
  //		remez(z,OPENARRAY(double,(0.00, 0.10, 0.20, 0.30, 0.40, 0.45, 0.55, 0.60, 0.70, 1.00)),
  //				    OPENARRAY(double,(1, 0, 1, 0, 1)), Weights->PValues(0), Weights->Length-1, rmtBandPass, err);

    //Design a hilbert (type III) transformer  (Sampling frequency = 2)
  //		  n = RemezLength(OPENARRAY(double,(0, 0, 0.1, 0.9, 1, 1)), OPENARRAY(double,(0,1,0)), OPENARRAY(double,(0.01,0.01,0.01)),Weights);
  //          if ((n%2) == 0) n++; //odd length type III filter
  //		  z->Size(n);
  //		  remez(z,OPENARRAY(double,(0.1, 0.9)), OPENARRAY(double,(1)), OPENARRAY(double,(1)), rmtHilbert,err);

    //Design a hilbert (type IV) transformer  (Sampling frequency = 2)
  //		  n = RemezLength(OPENARRAY(double,(0, 0, 0.1, 0.9, 1, 1)), OPENARRAY(double,(0,1,0)), OPENARRAY(double,(0.01,0.01,0.01)),Weights);
  //          if (!((n%2) == 0)) n++; //odd length type IV filter
  //          z->Size(n);
  //		  remez(z,OPENARRAY(double,(0.1, 1)), OPENARRAY(double,(1)), OPENARRAY(double,(1)), rmtHilbert,err);

    //Design a differentiator (type III) (Sampling frequency = 2)
  //		  n = RemezLength(OPENARRAY(double,(0, 0, 0.1, 0.9, 1, 1)), OPENARRAY(double,(0,1,0)), OPENARRAY(double,(0.01,0.01,0.01)),Weights);
  //          if ((n%2) == 0) n++; //odd length type III filter
  //          z->Size(n);
  //		  remez(z,OPENARRAY(double,(0.1, 0.9)), OPENARRAY(double,(1)), OPENARRAY(double,(1)), rmtDifferentiator,err);

    //Design a differentiator (type IV) (Sampling frequency = 2)
  //		  n = RemezLength(OPENARRAY(double,(0, 0, 0.1, 0.9, 1, 1)), OPENARRAY(double,(0,1,0)), OPENARRAY(double,(0.01,0.01,0.01)),Weights);
  //          if (!((n%2) == 0)) n++; //even length type IV filter
  //		  z->Size(n);
  //		  remez(z,OPENARRAY(double,(0.1, 1)), OPENARRAY(double,(1)), OPENARRAY(double,(1)), rmtDifferentiator,err);

    //Design a double differentiator (Sampling frequency = 2)
  //		  n = RemezLength(OPENARRAY(double,(0, 0, 0.05, 0.95, 1, 1)), OPENARRAY(double,(0,1,0)), OPENARRAY(double,(0.01,0.01,0.01)),Weights);
  //		  if (!((n%2) == 0)) n++; //even length type IV filter
  //		  z->Size(n);
  //		  remez(z,OPENARRAY(double,(0.05, 1)), OPENARRAY(double,(1)), OPENARRAY(double,(1)), rmtDoubleDifferentiator,err);

  //  Design an integrator (Sampling frequency = 2)
  //		  n = RemezLength(OPENARRAY(double,(0, 0, 0.05, 0.95, 1, 1)), OPENARRAY(double,(0,1,0)), OPENARRAY(double,(0.01,0.01,0.01)),Weights);
  //		  if (!((n%2) == 0)) n++; //even length type IV filter
  //		  z->Size(n);
  //		  remez(z,OPENARRAY(double,(0.05, 1)), OPENARRAY(double,(1)), OPENARRAY(double,(1)), rmtIntegrator,err);

    //Design a double integrator  (Sampling frequency = 2)
  //		  n = RemezLength(OPENARRAY(double,(0, 0, 0.05, 0.95, 1, 1)), OPENARRAY(double,(0,1,0)), OPENARRAY(double,(0.01,0.01,0.01)),Weights);
  //          if (!((n%2) == 0)) n++; //even length type IV filter
  //          z->Size(n);
  //		  remez(z,OPENARRAY(double,(0.05, 1)), OPENARRAY(double,(1)), OPENARRAY(double,(1)), rmtDoubleIntegrator,err);

      FrequencyResponse(z,NULL,Response,8);
      DrawIt(Response);

  }



    using Dew.Math;
    using Dew.Math.Editors;
    using Dew.Math.Units;
    using Dew.Signal;
    using Dew.Signal.Units;
    using Dew.Math.Tee;
    using Dew.Signal.Tee;

    private void button1_Click(object sender, EventArgs e)
    {
      Vector z = new Vector(0);
      Vector Response = new Vector(0);
      Vector Weights = new Vector(0);

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

//Design a lowpass
      n = OptimalFir.RemezLength(new double[4] {0, 1500, 2000, 4000},
                                 new double[2] {1, 0}, new double[2] {0.001, 0.001} ,Weights, 8000);

      z.Size(n);
      OptimalFir.remez( z, new double[4] { 0, 1500, 2000, 4000 },
                           new double[2] { 1, 0 }, (double []) Weights, TRemezType.rmtBandPass, out err, 8000, false);

      SignalUtils.FrequencyResponse(z, null, Response, 8, false, TSignalWindowType.wtRectangular, 0);
      TeeChart.DrawIt(Response, "Lowpass", false);

//Design a highpass
      n = OptimalFir.RemezLength(new double[4] { 0, 1500, 2000, 4000 }, 
                                 new double[2] { 0, 1 }, new double[2] { 0.001, 0.001 }, Weights, 8000);
      z.Size(n);
      OptimalFir.remez(z, new double[4] { 0, 1500, 2000, 4000 }, 
            new double[2] { 0, 1 }, (double[])Weights, TRemezType.rmtBandPass, out err, 8000, false);

      SignalUtils.FrequencyResponse(z, null, Response, 8, false, TSignalWindowType.wtRectangular, 0);
      TeeChart.DrawIt(Response, "Highpass", false);

//Design a bandpass (Sampling frequency = 2)
      n = OptimalFir.RemezLength(new double[6] {0, 0.15, 0.25, 0.45, 0.55, 1}, 
                                 new double[3] {0, 1, 0},new double[3] {0.001, 0.001, 0.001},Weights,FS);
      z.Size(n);
      OptimalFir.remez(z,new double[6] {0, 0.15, 0.25, 0.45, 0.55, 1},
                         new double[3] {0, 1, 0}, (double []) Weights, TRemezType.rmtBandPass, out err, FS, false);

      SignalUtils.FrequencyResponse(z, null, Response, 8, false, TSignalWindowType.wtRectangular, 0);
      TeeChart.DrawIt(Response, "Bandpass", false);

//Design a bandstop (Sampling frequency = 2)
      n = OptimalFir.RemezLength(new double[6] { 0, 0.15, 0.25, 0.45, 0.55, 1 }, 
                                 new double[3] { 1, 0, 1 }, new double[3] { 0.001, 0.001, 0.001 }, Weights,FS) ;
      z.Size(n);
      OptimalFir.remez(z, new double[6] { 0, 0.15, 0.25, 0.45, 0.55, 1 }, 
                          new double[3] { 1, 0, 1 }, (double[]) Weights, TRemezType.rmtBandPass, out err, FS, false);

      SignalUtils.FrequencyResponse(z, null, Response, 8, false, TSignalWindowType.wtRectangular, 0);
      TeeChart.DrawIt(Response, "Bandstop", false);

//Design a multiband (Sampling frequency = 2)
      n = OptimalFir.RemezLength(new double[10] {0.00, 0.10, 0.20, 0.30, 0.40, 0.45, 0.55, 0.60, 0.70, 1.00}, 
                                 new double[5] { 1, 0, 1, 0, 1 }, new double[5] { 0.001, 0.001, 0.001, 0.001, 0.001 }, Weights, FS);
      z.Size(n);
      OptimalFir.remez(z, new double[10] {0.00, 0.10, 0.20, 0.30, 0.40,0.45, 0.55, 0.60, 0.70, 1.00}, 
                          new double[5] { 1, 0, 1, 0, 1 }, (double[]) Weights, TRemezType.rmtBandPass, out err, FS, false);

      SignalUtils.FrequencyResponse(z, null, Response, 8, false, TSignalWindowType.wtRectangular, 0);
      TeeChart.DrawIt(Response, "Multibandpass", false);

//  Design a hilbert (type III) transformer  (Sampling frequency = 2)
      n = OptimalFir.RemezLength(new double[6] {0, 0, 0.1, 0.9, 1, 1},
                                 new double[3] {0,1,0},new double[3] {0.01,0.01,0.01},Weights,FS);
      if (n % 2 == 0) n++; //odd length type III filter
      z.Size(n);
      OptimalFir.remez(z,new double[2] {0.1, 0.9},
                         new double[1] {1}, new double[1] {1}, TRemezType.rmtHilbert,out err, FS, false);

      SignalUtils.FrequencyResponse(z, null, Response, 8, false, TSignalWindowType.wtRectangular, 0);
      TeeChart.DrawIt(Response, "Hilbert III", false);

//Design a hilbert (type IV) transformer  (Sampling frequency = 2)
      n = OptimalFir.RemezLength(new double[6] {0, 0, 0.1, 0.9, 1, 1},
                                 new double[3] {0,1,0},new double[3] {0.01,0.01,0.01},Weights,FS);
      if (n % 2 != 0) n++; //even length type IV filter
      z.Size(n);
      OptimalFir.remez(z, new double[2] {0.1, 1},
                          new double[1] {1},new double[1] {1}, TRemezType.rmtHilbert,out err, 2, false);

      SignalUtils.FrequencyResponse(z, null, Response, 8, false, TSignalWindowType.wtRectangular, 0);
      TeeChart.DrawIt(Response, "Hilbert IV", false);

//Design a differentiator (type III) (Sampling frequency = 2)
      n = OptimalFir.RemezLength(new double[6] {0, 0, 0.1, 0.9, 1, 1},
                                 new double[3] {0,1,0}, new double[3] {0.01,0.01,0.01},Weights,FS);
      if (n % 2 == 0) n++; //odd length type III filter
      z.Size(n);
      OptimalFir.remez(z,new double[2] {0.1, 0.9},
                         new double[1] {1},new double[1] {1}, TRemezType.rmtDifferentiator,out err, 2, false);

      SignalUtils.FrequencyResponse(z, null, Response, 8, false, TSignalWindowType.wtRectangular, 0);
      TeeChart.DrawIt(Response, "Differentiator III", false);  

//Design a differentiator (type IV) (Sampling frequency = 2)
      n = OptimalFir.RemezLength(new double[6] {0, 0, 0.1, 0.9, 1, 1},
                                 new double[3] {0,1,0},new double[3] {0.01,0.01,0.01},Weights, FS);
      if (n % 2 != 0) n++; //even length type IV filter
      z.Size(n);
      OptimalFir.remez(z,new double[2] {0.1, 1},
                         new double[1] {1},new double[1] {1}, TRemezType.rmtDifferentiator,out err, FS, false);

      SignalUtils.FrequencyResponse(z, null, Response, 8, false, TSignalWindowType.wtRectangular, 0);
      TeeChart.DrawIt(Response, "Differentiator IV", false);  

//Design a double differentiator (Sampling frequency = 2)
      n = OptimalFir.RemezLength(new double[6] {0, 0, 0.05, 0.95, 1, 1},
                                 new double[3] {0,1,0},new double[3] {0.01,0.01,0.01},Weights, FS);
      if (n %2 != 0) n++; //even length type IV filter
      z.Size(n);
      OptimalFir.remez(z,new double[2] {0.05, 1},
                         new double[1] {1},new double[1] {1}, TRemezType.rmtDoubleDifferentiator,out err, FS, false);

      SignalUtils.FrequencyResponse(z, null, Response, 8, false, TSignalWindowType.wtRectangular, 0);
      TeeChart.DrawIt(Response, "Double differentiator IV", false);  

//Design an integrator (Sampling frequency = 2)
      n = OptimalFir.RemezLength(new double[6] {0, 0, 0.05, 0.95, 1, 1},
                                 new double[3] {0,1,0},new double[3] {0.01,0.01,0.01}, Weights, FS);
      if (n %2 != 0) n++; //even length type IV filter
      z.Size(n);
      OptimalFir.remez(z,new double[2] {0.05, 1},
                         new double[1] {1},new double[1] {1}, TRemezType.rmtIntegrator,out err, FS, false);

      SignalUtils.FrequencyResponse(z, null, Response, 8, false, TSignalWindowType.wtRectangular, 0);
      TeeChart.DrawIt(Response, "Integrator IV", false);  

//Design a double integrator  (Sampling frequency = 2)
      n = OptimalFir.RemezLength(new double[6] {0, 0, 0.05, 0.95, 1, 1},
                                 new double[3] {0,1,0}, new double[3] {0.01,0.01,0.01},Weights,FS);
      if (n %2 != 0) n++; //even length type IV filter
      z.Size(n);
      OptimalFir.remez(z,new double[2] {0.05, 1},
                         new double[1] {1}, new double[1] {1}, TRemezType.rmtDoubleIntegrator,out err, FS, false);

      SignalUtils.FrequencyResponse(z,null,Response,8,false,TSignalWindowType.wtRectangular,0);
      TeeChart.DrawIt(Response,"Double integrator IV",false);
    }

Copyright 2008 Dew Research
http://www.dewresearch.com