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

Question

我试图在Java中实现低通滤波器.我的要求非常简单,我必须消除特定频率以外的信号(单维).看起来Butterworth过滤器可以满足我的需要.

现在重要的是CPU时间应该尽可能低.过滤器必须处理近百万个样本,我们的用户不喜欢等待太久.是否有任何现成的Butterworth滤波器实现,它具有最佳的滤波算法.

Answer 1

我有一个页面描述了一个非常简单,非常低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

甲smoothing的值1的原因没有发生,而较大的值越来越顺利出结果平滑.

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

Answer 2

我从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;
    }


}

Answer 3

我最近设计了一个简单的巴特沃斯函数（http://baumdevblog.blogspot.com/2010/11/butterworth-lowpass-filter-coefficients.html）。它们很容易用 Java 编码，并且如果你问我的话应该足够快（我猜你只需将 filter(double*样本，int count) 更改为 filter(double[]样本，int count)）。

JNI 的问题在于它会影响平台独立性，可能会混淆热点编译器，并且代码中的 JNI 方法调用可能仍然会减慢速度。所以我建议尝试 Java 并看看它是否足够快。

在某些情况下，首先使用快速傅立叶变换并在频域中应用滤波可能会有所帮助，但我怀疑这比简单低通滤波器的每个样本进行约 6 次乘法和一些加法要快。

Answer 4

过滤器设计是一门权衡的艺术，要想做得好，你需要考虑一些细节。

必须通过“不多”衰减的最大频率是多少，“不多”的最大值是多少？

必须“大量”衰减的最小频率是多少，“大量”的最小值是多少？

在滤波器应该通过的频率范围内，可接受的纹波（即衰减变化）有多大？

您有多种选择，这将花费您各种计算量。 像 matlab 或 scilab 这样的程序可以帮助您比较权衡。您需要熟悉一些概念，例如将频率表示为采样率的小数部分，以及衰减的线性和对数 (dB) 测量值之间的互换。

例如，“完美”的低通滤波器在频域中是矩形的。在时域中表示为脉冲响应，这将是一个 sinc 函数 (sin x/x)，其尾部达到正无穷大和负无穷大。显然你无法计算出来，所以问题就变成了如果你将 sinc 函数近似为你可以计算的有限持续时间，这会降低你的滤波器多少？

或者，如果您想要一个计算成本非常低的有限脉冲响应滤波器，您可以使用所有系数都为 1 的“箱式车”或矩形滤波器。利用二进制溢出来执行“循环”累加器，因为无论如何您稍后都会采用导数）。但是一个时间上呈矩形的滤波器在频率上看起来像一个 sinc 函数——它在通带中有一个 sin x/x 滚降（通常会提高到某个功率，因为​​你通常会有一个多级版本），还有一些“反弹”在阻带中。在某些情况下，它还是有用的，无论是单独使用还是与其他类型的过滤器配合使用。

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

就像 Mark Peters 在他的评论中所说：需要进行大量过滤的过滤器应该用 C 或 C++ 编写。但您仍然可以使用 Java。只需看一下Java 本机接口 (JNI)即可。由于 C/C++ 编译为本机代码，因此它的运行速度比在 Java 虚拟机 (JVM) 中运行字节码快得多，Java 虚拟机 (JVM) 实际上是一个虚拟处理器，它将字节码转换为本地机器的本机代码（取决于在 CPU 指令集上，如 x86、x64、ARM ......）