File content as of revision 5:184a0b023f5c:

//------------------------------------------------------------------------------
//  Design of FIR filter of LPF and HPF using window method
//
//  2014/11/09, Copyright (c) 2014 MIKAMI, Naoki
//------------------------------------------------------------------------------

#include "WindowingDesignLH.hpp"

namespace Mikami
{
    WindowingDesign::WindowingDesign(int order, float fs)
            : FS_(fs), PI_FS_(3.1415926536f/fs)
    {
        order_ = order;
        if ((order % 2) != 0)
        {
            fprintf(stderr, "order must be even.");
            return;
        }
        hm_ = new float[order/2+1];
        wn_ = new float[order/2+1];    
        
        HammWindow();
    }
        
    void WindowingDesign::Design(int order, Type pb, float fc,
                                 float hk[])
    {       
        if (pb == LPF) fC_ = fc;
        if (pb == HPF) fC_ = 0.5f*FS_ - fc;
    
        if (order != order_)
        {
            order_ = order_;
            if (hm_ != NULL) delete[] hm_;
            hm_ = new float[order/2+1];
            if (wn_ != NULL) delete[] wn_;
            wn_ = new float[order/2+1];
            HammWindow();
        }

        // Calculate coefficients for LPF
        LpfCoefficients();
        // If HPF, transform coefficients
        if (pb != LPF) ToHpf();
        
        for (int k=0; k<=order/2; k++)
            hk[k] = hm_[order/2-k];

        // For compensation of gain
        float sum = hk[order_/2];
        if (pb == LPF)
            for (int k=0; k<order_/2; k++)
                sum = sum + 2.0f*hk[k];
        if (pb == HPF)
        {
            float sign = -1.0;
            for (int k=order_/2-1; k>=0; k--)
            {
                sum = sum + sign*2.0f*hk[k];
                sign = -sign;
            }
        }
        for (int k=0; k<=order_/2; k++)
            hk[k] = hk[k]/sum;
    }

    // Calculation of coefficients for LPF
    void WindowingDesign::LpfCoefficients()
    {
        float w = 2.0f*PI_FS_*fC_;
        hm_[0] = 2.0f*fC_/FS_;
        for (int k=1; k<=order_/2; k++)
            hm_[k] = (sinf(k*w)/(PI_*k))*wn_[k];  
    }
    
    // Transform LPF to HPF
    void WindowingDesign::ToHpf()
    {
        for (int k=1; k<=order_/2; k+=2)
            hm_[k] = -hm_[k];
    }    

    // Hamming window
    void WindowingDesign::HammWindow()
    {
        float pi2OvM = 2.0f*PI_/(float)(order_ + 1);
        for (int n=0; n<=order_/2; n++)
            wn_[n] = 0.54f + 0.46f*cosf(pi2OvM*n);
    }
}