title: Simple bandpass filter in java
body: <p>I am trying to write a simple band pass filter following the instructions in <a href="http://www.dspguide.com/" rel="nofollow">this book</a>. My code creates a blackman window, and combines two low pass filter kernels to create a band pass filter kernel using spectral inversion, as described in the second example <a href="http://www.dspguide.com/ch16/4.htm" rel="nofollow">here</a> (table 16-2).</p>
  
  <p>I am testing my code by comparing it with the results I get in matlab. When I test the methods that create a blackman window and a low pass filter kernel separately, I get results that are close to what I see in matlab (up to some digits after the decimal point - I attribute the error to java double variables rounding issues), but my band pass filter kernel is incorrect.</p>
  
  <p>Tests I ran:</p>
  
  <ul>
  <li>Created a blackman window and compared it with what I get in matlab - all good.</li>
  <li>Created a low pass filter using this window using my code and <code>fir1(N, Fc1/(Fs/2), win, flag);</code> in matlab (see full code below). I think the results are correct, although I get bigger error the bigger Fc1 is (<strong>why?</strong>)</li>
  <li>Created a pand pass filter using my code and <code>fir1(N, [Fc1 Fc2]/(Fs/2), 'bandpass', win, flag);</code> in matlab - results are completely off.</li>
  <li>Filtered my data using my code and the kernel generated by matlab - all good.</li>
  </ul>
  
  <p>So - why is my band pass filter kernel off? What did I do wrong?
  I think I either have a bug or fir1 uses a different algorithm, but I can't check because the <a href="http://www.mathworks.com/help/signal/ref/fir1.html#bqih_e6-1" rel="nofollow">article referenced in its documentation</a> is not publicly available.</p>
  
  <p>This is my matlab code:</p>
  
  <pre><code>Fs = 200;       % Sampling Frequency
  N    = 10;          % Order
  Fc1  = 1.5;         % First Cutoff Frequency
  Fc2  = 7.5;         % Second Cutoff Frequency
  flag = 'scale';     % Sampling Flag
  
  % Create the window vector for the design algorithm.
  win = blackman(N+1);
  
  % Calculate the coefficients using the FIR1 function.
  b  = fir1(N, [Fc1 Fc2]/(Fs/2), 'bandpass', win, flag);
  Hd = dfilt.dffir(b);
  res = filter(Hd, data);
  </code></pre>
  
  <p>This is my java code (I believe the bug is in bandPassKernel):</p>
  
  <pre><code>/**
       * See - http://www.mathworks.com/help/signal/ref/blackman.html
       * @param length
       * @return
       */
      private static double[] blackmanWindow(int length) {
  
          double[] window = new double[length];
          double factor = Math.PI / (length - 1);
  
          for (int i = 0; i &lt; window.length; ++i) {
              window[i] = 0.42d - (0.5d * Math.cos(2 * factor * i)) + (0.08d * Math.cos(4 * factor * i));
          }
  
          return window;
      }
  
  private static double[] lowPassKernel(int length, double cutoffFreq, double[] window) {
  
      double[] ker = new double[length + 1];
      double factor = Math.PI * cutoffFreq * 2; 
      double sum = 0;
  
      for (int i = 0; i &lt; ker.length; i++) {
          double d = i - length/2; 
          if (d == 0) ker[i] = factor;
          else ker[i] =  Math.sin(factor * d) / d;
          ker[i] *= window[i];
          sum += ker[i];
      }
  
      // Normalize the kernel
      for (int i = 0; i &lt; ker.length; ++i) {
          ker[i] /= sum;
      }
  
      return ker;
  }
  
  private static double[] bandPassKernel(int length, double lowFreq, double highFreq) {
  
      double[] ker = new double[length + 1];
      double[] window = blackmanWindow(length + 1);
  
      // Create a band reject filter kernel using a high pass and a low pass filter kernel 
      double[] lowPass = lowPassKernel(length, lowFreq, window);
  
      // Create a high pass kernel for the high frequency
      // by inverting a low pass kernel
      double[] highPass = lowPassKernel(length, highFreq, window);
      for (int i = 0; i &lt; highPass.length; ++i) highPass[i] = -highPass[i];
      highPass[length / 2] += 1;
  
      // Combine the filters and invert to create a bandpass filter kernel
      for (int i = 0; i &lt; ker.length; ++i) ker[i] = -(lowPass[i] + highPass[i]);
      ker[length / 2] += 1;
  
      return ker;
  }
  
  private static double[] filter(double[] signal, double[] kernel) {
  
      double[] res = new double[signal.length];
  
      for (int r = 0; r &lt; res.length; ++r) {
  
          int M = Math.min(kernel.length, r + 1);
          for (int k = 0; k &lt; M; ++k) {
              res[r] += kernel[k] * signal[r - k];
          }
      }
  
      return res;
  }
  </code></pre>
  
  <p>And this is how I use my code:</p>
  
  <pre><code>double[] kernel = bandPassKernel(10, 1.5d / (200/2), 7.5d / (200/2));
  double[] res = filter(data, kernel);
  </code></pre>