// nonrec.cpp -- non-recursive filter design
 
//     NRLowPass()          build NR low-pass filter
//     NRHighPass()         build NR high-pass filter
//     NRBandPass()         build NR band-pass filter

#include "tools.h"
#include "ltisys.h"

// sinc() or sin(x)/x function
double sinc(double x)
{
       if (x==0)
              return cos(x);
       else
              return sin(x)/x;
}

// Create non-recursive low-pass filter by Fourier Transform method
LTISystem NRLowPass(
       double freq,         // corner freq (fraction of sampling rate)
       int ncoeff           // # coefficients
       )                    // returns LTI system to specification
{
       int    i;

       // convert frequency
       double omega = 2 * PI * freq;

       // build half-sized window from sinc function
       int nhalf = ncoeff/2;
       Waveform hwv(nhalf,1.0);
       for (i=1;i<=nhalf;i++)
              hwv[i] = omega * sinc(i*omega)/PI;

       // window with (half-)hamming window
       for (i=1;i<=nhalf;i++)
              hwv[i] *= 0.54 + 0.46 * cos(i*PI/nhalf);

       // create new LTI System
       LTISystem lpfilt(2*nhalf,0);

       // copy impulse response into system
       //   (indexed -nhalf .. 0 .. nhalf)
       lpfilt.a[nhalf] = omega/PI;
       for (i=1;i<=nhalf;i++) {
              lpfilt.a[nhalf-i] = hwv[i];
              lpfilt.a[nhalf+i] = hwv[i];
       }

       return lpfilt;
}

// create high-pass non-recursive filter from low-pass prototype
LTISystem NRHighPass(
       double freq,         // corner freq (fraction of sampling rate)
       int ncoeff           // # coefficients
       )                    // returns LTI system
{
       // get low-pass version
       LTISystem hpfilt = NRLowPass(0.5-freq,ncoeff);

       // now modulate with cos(n*PI) = +1,-1,+1,-1,...
       int nhalf = hpfilt.a.count()/2;
       for (int i=1;i<=nhalf;i+=2) {
              hpfilt.a[nhalf-i] = -hpfilt.a[nhalf-i];
              hpfilt.a[nhalf+i] = -hpfilt.a[nhalf+i];
       }

       return hpfilt;
}

// create band-pass non-recursive filter from low-pass prototype
LTISystem NRBandPass(
       double lofreq,       // low corner freq (fraction of sampling rate) 
       double hifreq,       // high corner freq (fraction of sampling rate) }
       int ncoeff           // # coefficients
       )                    // returns LTI system
{
       // get low-pass prototype
       LTISystem bpfilt = NRLowPass((hifreq-lofreq)/2,ncoeff);

       // now modulate with cos(n*centrefreq)
       int nhalf = bpfilt.a.count()/2;
       double cf = 2*PI*(lofreq+hifreq)/2;
       bpfilt.a[nhalf] = 2 * bpfilt.a[nhalf];
       for (int i=1;i<=nhalf;i++) {
              bpfilt.a[nhalf-i] = 2 * bpfilt.a[nhalf-i] * cos(i*cf);
              bpfilt.a[nhalf+i] = 2 * bpfilt.a[nhalf+i] * cos(i*cf);
       }

       return bpfilt;
}