matlab - Simple bandpass filter in java -


i trying write simple band pass filter following instructions in this book. code creates blackman window, , combines 2 low pass filter kernels create band pass filter kernel using spectral inversion, described in second example here (table 16-2).

i testing code comparing results in matlab. when test methods create blackman window , low pass filter kernel separately, results close see in matlab (up digits after decimal point - attribute error java double variables rounding issues), band pass filter kernel incorrect.

tests ran:

  • created blackman window , compared in matlab - good.
  • created low pass filter using window using code , fir1(n, fc1/(fs/2), win, flag); in matlab (see full code below). think results correct, although bigger error bigger fc1 (why?)
  • created pand pass filter using code , fir1(n, [fc1 fc2]/(fs/2), 'bandpass', win, flag); in matlab - results off.
  • filtered data using code , kernel generated matlab - good.

so - why band pass filter kernel off? did wrong? think either have bug or fir1 uses different algorithm, can't check because article referenced in documentation not publicly available.

this matlab 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 window vector design algorithm. win = blackman(n+1);  % calculate coefficients using fir1 function. b  = fir1(n, [fc1 fc2]/(fs/2), 'bandpass', win, flag); hd = dfilt.dffir(b); res = filter(hd, data); 

this java code (i believe bug in bandpasskernel):

/**      * 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);          (int = 0; < 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;      (int = 0; < ker.length; i++) {         double d = - length/2;          if (d == 0) ker[i] = factor;         else ker[i] =  math.sin(factor * d) / d;         ker[i] *= window[i];         sum += ker[i];     }      // normalize kernel     (int = 0; < 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 band reject filter kernel using high pass , low pass filter kernel      double[] lowpass = lowpasskernel(length, lowfreq, window);      // create high pass kernel high frequency     // inverting low pass kernel     double[] highpass = lowpasskernel(length, highfreq, window);     (int = 0; < highpass.length; ++i) highpass[i] = -highpass[i];     highpass[length / 2] += 1;      // combine filters , invert create bandpass filter kernel     (int = 0; < 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];      (int r = 0; r < res.length; ++r) {          int m = math.min(kernel.length, r + 1);         (int k = 0; k < m; ++k) {             res[r] += kernel[k] * signal[r - k];         }     }      return res; } 

and how use code:

double[] kernel = bandpasskernel(10, 1.5d / (200/2), 7.5d / (200/2)); double[] res = filter(data, kernel); 

i ended implementing matlab's fir1 function in java. results quite accurate.


Comments