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
Post a Comment