Constructing FIR Digital Filters with valarray

Figure 2: Partial listing of file FIR.h, the FIR filter implementation

const double pi = 3.14159265359;
const complex<double> j(0,1);            // Imaginary unit

#include "filter_types.h"
#include "conversions.h"
#include "Matrix.h"

template <class T>
class basic_FIR
{
public:
    struct Constraint_point { double w, gain; };    
        // For the frequency-sampling filter design method

        // ***************  Constructors  ********************
    basic_FIR (valarray<T> coefficients)
        : h (coefficients)
    { reverse (&h[0], &h[0] + h.size()); }

    basic_FIR (const T * coefficients, size_t N);
        // Not shown here -- similar to the previous one

    basic_FIR (struct Low_pass, double scale_factor = 1);
    basic_FIR (struct High_pass, double scale_factor = 1);
    basic_FIR (struct Band_pass, double scale_factor = 1);

    basic_FIR (double (*frequency_response) (double), 
               int N, double scale_factor = 1);
    basic_FIR (const vector <struct Constraint_point> & H, 
               double scale_factor = 1);

        // *************  Filtering operations  ***************
    template <class Iterator, class Type>
    void output (Iterator sample, Type & result) const
    {
            // Use convert in case rounding is required
        convert (inner_product(sample-h.size()+1, sample+1, 
                               &(const_cast<basic_FIR *>(this)->h[0]), 
                               Type()), 
                 result);
    }

    template <class Iterator, class IteratorOut>
    void operator() (Iterator first, Iterator last, 
                     IteratorOut OutputSignal) const
    {
        for (Iterator current = first; current != last; ++current)
            output (current, *OutputSignal++);
    }

        // ************  Misc. member functions  *************

    complex<double> H (double w) const;   // Frequency response

    int length () const  { return h.size(); }

private:
    valarray<T> h;      // filter coefficients

    void compute_coefficients(const vector<Constraint_point>&);
};

typedef basic_FIR<double> FIR;


// *************************************
//      Member functions definitions
// *************************************

template <class T>
basic_FIR<T>::basic_FIR (struct Low_pass filter_descriptor, 
                         double scale_factor)
{
    if (filter_descriptor.N < 2)
        filter_descriptor.N = 2;    // N can not be < 2

    vector<Constraint_point> H (filter_descriptor.N);

    const double Wc = filter_descriptor.Fc * 2 * pi / 
                      filter_descriptor.Fs;

    for (int i = 0; i < H.size(); i++)
    {
        const double delta_w = pi / (H.size() - 1);     
                // delta_w is the difference in w between two 
                // consecutive frequency samples

        H[i].w = i * delta_w;

        if (H[i].w <= Wc - 2 * delta_w)     // In the pass-band
        {
            H[i].gain = 1 * scale_factor;
        }
        else if (H[i].w >= Wc + 2 * delta_w)    // stop-band
        {
            H[i].gain = 0;
        }
        else        // Transition band (near Wc)
        {           // Use a raised cosine to calculate values 
            H[i].gain = (1 + cos(pi * (H[i].w - (Wc - 2 * delta_w)) / 
                        (4 * delta_w))) * scale_factor / 2;
        }
    }

    compute_coefficients (H);
}

template <class T>
basic_FIR<T>::basic_FIR (double (*frequency_response) (double), int N, double scale_factor)
{
    if (N < 2)  N = 2;      // N can't be < 2

    vector <struct Constraint_point> H(N);

    for (int i = 0; i < N; i++)
    {
        H[i].w = pi * i / (N - 1);
        H[i].gain = scale_factor * frequency_response (H[i].w);
    }

    compute_coefficients (H);  // compute_coefficients resizes h
}

// *** Other constructors not shown here -- see note below 

template <class T>
complex<double> basic_FIR<T>::H (double w) const
{
    complex<double> result = complex<double> (0,0);
    for (int n = 0; n < h.size(); n++)
        result +=  static_cast<double>(h[n]) * exp (-j * (w * n));

    return result;
}

template <class T>
void basic_FIR<T>::compute_coefficients 
                   (const vector<struct Constraint_point> & H)
{
        // Use Gauss method to solve the linear set of eqs.

    const int N = H.size();

    Matrix A (N, N);
    valarray<long double> x(N), b(N);

    for (int n = 0; n < N; n++)
    {
        A[n][0] = 1;
        for (int k = 1; k < N; k++)
            A[n][k] = 2 * cos(k * H[n].w);

        b[n] = H[n].gain;
    }

    Gauss (A,x,b);

    h.resize (2*N - 1);

    convert (x[0], h[N-1]);
    for (int i = 1; i < N; i++)
    {
        convert (x[i], h[N-1 + i]);
        h[N-1 - i] = h[N-1 + i];    // Symmety condition
    }
}

// Some of the functions are omited here (as well as in the 
// next listings) to save magazine space.  The complete version 
// of the code can be downloaded from the Journal's web site, 
// or from my my web site, at http://www.mochima.com/downloads