如何使用java实现低通滤波器

Question

我想在Java中实现低通滤波器。我的要求非常简单，我必须消除特定频率以外的信号（单维）。看起来巴特沃斯过滤器会适合我的需要。

现在重要的是CPU时间应该尽可能低。过滤器需要处理的数量将接近一百万个，我们的用户不喜欢等待太久。是否有任何现成的Butterworth滤波器实现，它具有最佳的滤波算法。

问候，

Chaitannya

Answer 1

像马克·彼得斯在他的评论中说：这需要过滤大量的过滤器应该用C或C++。但是你仍然可以使用Java。只要看看Java Native Interface (JNI)。由于C/C++编译为本地机器代码，它将比运行Java虚拟机（JVM）中的字节码快很多，它实际上是一个虚拟处理器，它将字节码转换为本地机器的本地代码（取决于在CPU指令集如x86，x64，ARM，....）

Answer 2

我设计了一个简单的butterworth函数最近（http://baumdevblog.blogspot.com/2010/11/butterworth-lowpass-filter-coefficients.html）。他们很容易在Java中编码，如果你问我（你只需要改变过滤器（double * samples，int count）来过滤（double [] samples，int count），我猜应​​该足够快）。

JNI的问题是，它会降低平台的独立性，可能会混淆热点编译器和代码中的JNI方法调用，可能会让速度变慢。所以我会建议尝试Java，看看它是否足够快。

在某些情况下，首先使用快速傅立叶变换并应用频域滤波可能是有益的，但我怀疑这比每个样本对于简单低通滤波器的约6次乘法和几次加法更快。

Answer 3

滤波器设计是一种权衡的艺术，为了做到这一点，您需要考虑一些细节。

什么是必须通过的“没有太多”的最大频率，什么是“没有太多”的最大值？

什么是必须衰减的最低频率“很多”，什么是“很多”的最小值？

滤波器应该通过的频率内有多少波纹（即衰减的变化）是可以接受的？

您有广泛的选择，这将花费您不同的计算量。 像matlab或scilab这样的程序可以帮助您比较折衷。您需要熟悉概念，例如将频率表示为采样率的小数部分，以及线性和对数（dB）衰减测量之间的交换。

例如，“完美”低通滤波器在频域中是矩形的。在时域中表示为脉冲响应，这将是一个正弦函数（sin x/x），其尾部达到正负无穷大。显然你不能计算这个问题，所以问题就变成了如果你将sinc函数逼近到你可以计算的有限持续时间，那么你的滤波器会降低多少？或者，如果你想要一个非常便宜的有限脉冲响应滤波器，你可以使用一个“盒子车”或矩形滤波器，其中所有的系数都是1.（如果你实现它，它可以更便宜作为一个CIC滤波器，利用二进制溢出来做'循环'累加器，因为无论如何你都会接受这个导数）。但是时间上呈矩形的滤波器看起来像频率上的正弦函数 - 它在通带中有一个正弦x/x滚降（通常由于通常会有多级版本而升高到某些功率），并且一些“反弹”在阻带中。仍然在某些情况下，它本身很有用，或者在其他类型的过滤器后面进行跟踪。

Answer 4

我有一个页面描述了一个非常简单的，非常低CPU的低通滤波器，它也能够独立于帧率。我使用它来平滑用户输入，并经常用于绘制帧频。

http://phrogz.net/js/framerate-independent-low-pass-filter.html

总之，在您的更新循环：

// If you have a fixed frame rate 
smoothedValue += (newValue - smoothedValue)/smoothing 

// If you have a varying frame rate 
smoothedValue += timeSinceLastUpdate * (newValue - smoothedValue)/smoothing 

的1原因没有发生，而较高的值越来越顺利进行，结果上平滑smoothing值。

该页面有几个用JavaScript编写的函数，但公式与语言无关。

Answer 5

这是一个在apache数学库中使用傅里叶变换的低通滤波器。

public double[] fourierLowPassFilter(double[] data, double lowPass, double frequency){ 
    //data: input data, must be spaced equally in time. 
    //lowPass: The cutoff frequency at which 
    //frequency: The frequency of the input data. 

    //The apache Fft (Fast Fourier Transform) accepts arrays that are powers of 2. 
    int minPowerOf2 = 1; 
    while(minPowerOf2 < data.length) 
     minPowerOf2 = 2 * minPowerOf2; 

    //pad with zeros 
    double[] padded = new double[minPowerOf2]; 
    for(int i = 0; i < data.length; i++) 
     padded[i] = data[i]; 


    FastFourierTransformer transformer = new FastFourierTransformer(DftNormalization.STANDARD); 
    Complex[] fourierTransform = transformer.transform(padded, TransformType.FORWARD); 

    //build the frequency domain array 
    double[] frequencyDomain = new double[fourierTransform.length]; 
    for(int i = 0; i < frequencyDomain.length; i++) 
     frequencyDomain[i] = frequency * i/(double)fourierTransform.length; 

    //build the classifier array, 2s are kept and 0s do not pass the filter 
    double[] keepPoints = new double[frequencyDomain.length]; 
    keepPoints[0] = 1; 
    for(int i = 1; i < frequencyDomain.length; i++){ 
     if(frequencyDomain[i] < lowPass) 
      keepPoints[i] = 2; 
     else 
      keepPoints[i] = 0; 
    } 

    //filter the fft 
    for(int i = 0; i < fourierTransform.length; i++) 
     fourierTransform[i] = fourierTransform[i].multiply((double)keepPoints[i]); 

    //invert back to time domain 
    Complex[] reverseFourier = transformer.transform(fourierTransform, TransformType.INVERSE); 

    //get the real part of the reverse 
    double[] result = new double[data.length]; 
    for(int i = 0; i< result.length; i++){ 
     result[i] = reverseFourier[i].getReal(); 
    } 

    return result; 
} 

Answer 6

我采用了这种从http://www.dspguide.com/ 我很新到Java，所以它不是很漂亮，但它的工作原理

/* 
* To change this license header, choose License Headers in Project Properties. 
* To change this template file, choose Tools | Templates 
* and open the template in the editor. 
*/ 
package SoundCruncher; 

import java.util.ArrayList; 

/** 
* 
* @author 2sloth 
* filter routine from "The scientist and engineer's guide to DSP" Chapter 20 
* filterOrder can be any even number between 2 & 20 


* cutoffFreq must be smaller than half the samplerate 
* filterType: 0=lowPass 1=highPass 
* ripplePercent is amount of ripple in Chebyshev filter (0-29) (0=butterworth) 
*/ 
public class Filtering { 
    double[] filterSignal(ArrayList<Float> signal, double sampleRate ,double cutoffFreq, double filterOrder, int filterType, double ripplePercent) { 
     double[][] recursionCoefficients = new double[22][2]; 
     // Generate double array for ease of coding 
     double[] unfilteredSignal = new double[signal.size()]; 
     for (int i=0; i<signal.size(); i++) { 
      unfilteredSignal[i] = signal.get(i); 
     } 

     double cutoffFraction = cutoffFreq/sampleRate; // convert cut-off frequency to fraction of sample rate 
     System.out.println("Filtering: cutoffFraction: " + cutoffFraction); 
     //ButterworthFilter(0.4,6,ButterworthFilter.Type highPass); 
     double[] coeffA = new double[22]; //a coeffs 
     double[] coeffB = new double[22]; //b coeffs 
     double[] tA = new double[22]; 
     double[] tB = new double[22]; 

     coeffA[2] = 1; 
     coeffB[2] = 1; 

     // calling subroutine 
     for (int i=1; i<filterOrder/2; i++) { 
      double[] filterParameters = MakeFilterParameters(cutoffFraction, filterType, ripplePercent, filterOrder, i); 

      for (int j=0; j<coeffA.length; j++){ 
       tA[j] = coeffA[j]; 
       tB[j] = coeffB[j]; 
      } 
      for (int j=2; j<coeffA.length; j++){ 
       coeffA[j] = filterParameters[0]*tA[j]+filterParameters[1]*tA[j-1]+filterParameters[2]*tA[j-2]; 
       coeffB[j] = tB[j]-filterParameters[3]*tB[j-1]-filterParameters[4]*tB[j-2]; 
      } 
     } 
     coeffB[2] = 0; 
     for (int i=0; i<20; i++){ 
      coeffA[i] = coeffA[i+2]; 
      coeffB[i] = -coeffB[i+2]; 
     } 

     // adjusting coeffA and coeffB for high/low pass filter 
     double sA = 0; 
     double sB = 0; 
     for (int i=0; i<20; i++){ 
      if (filterType==0) sA = sA+coeffA[i]; 
      if (filterType==0) sB = sB+coeffB[i]; 
      if (filterType==1) sA = sA+coeffA[i]*Math.pow(-1,i); 
      if (filterType==1) sB = sB+coeffA[i]*Math.pow(-1,i); 
     } 

     // applying gain 
     double gain = sA/(1-sB); 
     for (int i=0; i<20; i++){ 
      coeffA[i] = coeffA[i]/gain; 
     } 
     for (int i=0; i<22; i++){ 
      recursionCoefficients[i][0] = coeffA[i]; 
      recursionCoefficients[i][1] = coeffB[i]; 
     } 
     double[] filteredSignal = new double[signal.size()]; 
     double filterSampleA = 0; 
     double filterSampleB = 0; 

     // loop for applying recursive filter 
     for (int i= (int) Math.round(filterOrder); i<signal.size(); i++){ 
      for(int j=0; j<filterOrder+1; j++) { 
       filterSampleA = filterSampleA+coeffA[j]*unfilteredSignal[i-j]; 
      } 
      for(int j=1; j<filterOrder+1; j++) { 
       filterSampleB = filterSampleB+coeffB[j]*filteredSignal[i-j]; 
      } 
      filteredSignal[i] = filterSampleA+filterSampleB; 
      filterSampleA = 0; 
      filterSampleB = 0; 
     } 


     return filteredSignal; 

    } 
    /* pi=3.14... 
     cutoffFreq=fraction of samplerate, default 0.4 FC 
     filterType: 0=LowPass 1=HighPass    LH 
     rippleP=ripple procent 0-29      PR 
     iterateOver=1 to poles/2      P% 
    */ 
    // subroutine called from "filterSignal" method 
    double[] MakeFilterParameters(double cutoffFraction, int filterType, double rippleP, double numberOfPoles, int iteration) { 
     double rp = -Math.cos(Math.PI/(numberOfPoles*2)+(iteration-1)*(Math.PI/numberOfPoles)); 
     double ip = Math.sin(Math.PI/(numberOfPoles*2)+(iteration-1)*Math.PI/numberOfPoles); 
     System.out.println("MakeFilterParameters: ripplP:"); 
      System.out.println("cutoffFraction filterType rippleP numberOfPoles iteration"); 
      System.out.println(cutoffFraction + " " + filterType + " " + rippleP + " " + numberOfPoles + " " + iteration); 
     if (rippleP != 0){ 
      double es = Math.sqrt(Math.pow(100/(100-rippleP),2)-1); 
//   double vx1 = 1/numberOfPoles; 
//   double vx2 = 1/Math.pow(es,2)+1; 
//   double vx3 = (1/es)+Math.sqrt(vx2); 
//   System.out.println("VX's: "); 
//   System.out.println(vx1 + " " + vx2 + " " + vx3); 
//   double vx = vx1*Math.log(vx3); 
      double vx = (1/numberOfPoles)*Math.log((1/es)+Math.sqrt((1/Math.pow(es,2))+1)); 
      double kx = (1/numberOfPoles)*Math.log((1/es)+Math.sqrt((1/Math.pow(es,2))-1)); 
      kx = (Math.exp(kx)+Math.exp(-kx))/2; 
      rp = rp*((Math.exp(vx)-Math.exp(-vx))/2)/kx; 
      ip = ip*((Math.exp(vx)+Math.exp(-vx))/2)/kx; 
      System.out.println("MakeFilterParameters (rippleP!=0):"); 
      System.out.println("es vx kx rp ip"); 
      System.out.println(es + " " + vx*100 + " " + kx + " " + rp + " " + ip); 
     } 

     double t = 2*Math.tan(0.5); 
     double w = 2*Math.PI*cutoffFraction; 
     double m = Math.pow(rp, 2)+Math.pow(ip,2); 
     double d = 4-4*rp*t+m*Math.pow(t,2); 
     double x0 = Math.pow(t,2)/d; 
     double x1 = 2*Math.pow(t,2)/d; 
     double x2 = Math.pow(t,2)/d; 
     double y1 = (8-2*m*Math.pow(t,2))/d; 
     double y2 = (-4-4*rp*t-m*Math.pow(t,2))/d; 
     double k = 0; 
     if (filterType==1) { 
      k = -Math.cos(w/2+0.5)/Math.cos(w/2-0.5); 
     } 
     if (filterType==0) { 
      k = -Math.sin(0.5-w/2)/Math.sin(w/2+0.5); 
     } 
     d = 1+y1*k-y2*Math.pow(k,2); 
     double[] filterParameters = new double[5]; 
     filterParameters[0] = (x0-x1*k+x2*Math.pow(k,2))/d;   //a0 
     filterParameters[1] = (-2*x0*k+x1+x1*Math.pow(k,2)-2*x2*k)/d; //a1 
     filterParameters[2] = (x0*Math.pow(k,2)-x1*k+x2)/d;   //a2 
     filterParameters[3] = (2*k+y1+y1*Math.pow(k,2)-2*y2*k)/d;  //b1 
     filterParameters[4] = (-(Math.pow(k,2))-y1*k+y2)/d;   //b2 
     if (filterType==1) { 
      filterParameters[1] = -filterParameters[1]; 
      filterParameters[3] = -filterParameters[3]; 
     } 
//  for (double number: filterParameters){ 
//   System.out.println("MakeFilterParameters: " + number); 
//  } 


     return filterParameters; 
    } 


}