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

How to implement low pass filter using java

提问by chai

I am trying to implement a low pass filter in Java. My requirement is very simple,I have to eliminate signals beyond a particular frequency (Single dimension). Looks like Butterworth filter would suit my need.

我正在尝试用 Java 实现低通滤波器。我的要求很简单，我必须消除超出特定频率（单维）的信号。看起来 Butterworth 过滤器适合我的需要。

Now the important thing is that CPU time should be as low as possible. There would be close to a million sample the filter would have to process and our users don't like waiting too long. Are there any readymade implementation of Butterworth filters which has optimal algorithms for filtering.

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

回答by Martijn Courteaux

Like Mark Peters said in his comment: A filter which needs to filter a lot should be written in C or C++. But you can still make use of Java. Just take a look at Java Native Interface (JNI). Because of C/C++ compiles to native machine code, it will run a lot faster than running your bytecode in the Java Virtual Machine (JVM), which is in fact a virtual processor that translates the bytecode to the local machine its native code (depending on CPU instruction set like x86, x64, ARM, ....)

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

回答by Boom

I have designed a simple butterworth function recently (http://baumdevblog.blogspot.com/2010/11/butterworth-lowpass-filter-coefficients.html). They are easy to code in Java and should be fast enough if you ask me (you'd just have to change filter(double* samples, int count) to filter(double[] samples, int count), I guess).

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

The problem with JNI is that it costs platform independence, may confuse the hotspot compiler and the JNI method calls within your code may still slow things down. So I would recommend trying Java and see if it is fast enough.

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

In some cases it might be beneficial to use a fast fourier transform first and apply the filtering in the frequency domain but I doubt that this is faster than about 6 multiplies and a few additions per sample for a simple lowpass filter.

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

回答by Chris Stratton

Filter design is an art of tradeoffs, and to do it well you need to take some details into account.

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

What is the maximum frequency which must be passed "without much" attentuation, and what is the maximum value of "without much" ?

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

What is the minimum frequency which must be attenuated "a lot" and what is the minimum value of "a lot" ?

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

How much ripple (ie variation in attenuation) is acceptable within the frequencies the filter is supposed to pass?

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

You have a wide range of choices, which will cost you a variety of amounts of computation. A program like matlab or scilab can help you compare the tradeoffs. You'll want to become familiar with concepts like expressing frequencies as a decimal fraction of a sample rate, and interchanging between linear and log (dB) measurements of attenuation.

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

For example, a "perfect" low pass filter is rectangular in the frequency domain. Expressed in the time domain as an impulse response, that would be a sinc function (sin x/x) with the tails reaching to both positive and negative infinity. Obviously you can't calculate that, so the question becomes if you approximate the sinc function to a finite duration which you can calculate, how much does that degrade your filter?

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

Alternately, if you want a finite impulse response filter that is very cheap to calculate, you can use a "box car" or rectangular filter where all the coefficients are 1. (This can be made even cheaper if you implement it as a CIC filter exploiting binary overflow to do 'circular' accumulators, since you'll be taking the derivative later anyway). But a filter that is rectangular in time looks like a sinc function in frequency - it has a sin x/x rolloff in the passband (often raised to some power since you would typically have a multi stage version), and some "bounce back" in the stop band. Still in some cases it's useful, either by itself or when followed up by another type of filter.

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

回答by Phrogz

I have a page describing a very simple, very low-CPU low-pass filter that is also able to be framerate-independent. I use it for smoothing out user input and also for graphing frame rates often.

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

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

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

In short, in your update loop:

简而言之，在您的更新循环中：

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

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

A smoothingvalue of 1causes no smoothing to occur, while higher values increasingly smooth out the result.

甲smoothing的值1的原因没有平滑发生，而较大的值越来越平滑结果。

The page has a couple of functions written in JavaScript, but the formula is language agnostic.

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

回答by Mark

Here is a lowpass filter that uses a fourier transform in the apache math library.

这是一个使用 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;
}

回答by 2sloth

I adopted this from http://www.dspguide.com/I'm quite new to java, so it's not pretty, but it works

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


}

回答by Jo?o Manuel Rodrigues

I understand this is an old question, but I'd like to add something that does not seem to have been mentioned before.

我知道这是一个老问题，但我想补充一些以前似乎没有提到的内容。

The first thing you should realize is that:

你首先应该意识到的是：

1. There are several different types of (digital) filters.
2. There are several different ways to design filters.
3. There are several different implementations of filters.

1. 有几种不同类型的（数字）滤波器。
2. 有几种不同的方法来设计滤波器。
3. 过滤器有几种不同的实现方式。

When you need to use a filter in your application, you'll have to choose a certain type of filter, choose a specific design method for that kind of filter, apply that method to find the filter coefficients that satisfy your constraints and, finally, copy those coefficients to your filter implementation.

当您需要在应用程序中使用滤波器时，您必须选择某种类型的滤波器，为该种滤波器选择一种特定的设计方法，应用该方法来找到满足您的约束条件的滤波器系数，最后，将这些系数复制到您的滤波器实现中。

Choosing the type of filter and applying the design method is something you do once, using an appropriate software, such as Matlab or Octave or Scilab. This may involve a bit of experimentation and observation of the obtained characteristics of the filter, such as the frequency response (both amplitude and phase) and/or the impulse response, to see if they match your specifications. Once you settle on the solution, you will have a set of constant coefficients. These numbers, or some linear combination of these numbers, is all you need to copy to your program (Java or otherwise) as a table of constants.

选择滤波器类型和应用设计方法只需使用适当的软件（例如 Matlab、Octave 或 Scilab）即可完成。这可能涉及对获得的滤波器特性（例如频率响应（幅度和相位）和/或脉冲响应）的一些实验和观察，以查看它们是否符合您的规格。一旦您确定了解决方案，您将拥有一组常数系数。这些数字，或这些数字的某种线性组合，是您需要将其作为常量表复制到您的程序（Java 或其他）中的全部内容。

Then, in your program, you just need a function that applies some filter implementation that uses these coefficients to do linear combinations of the input stream samples (and possibly previous output samples) to produce the new output sample in each time instant.

然后，在您的程序中，您只需要一个应用一些滤波器实现的函数，该函数使用这些系数对输入流样本（可能还有之前的输出样本）进行线性组合，以在每个时刻产生新的输出样本。

I assume you will probably be interested in Linear Time-Invariant filters, either with finite impulse response (FIR) or infinite impulse response (IIR). The design methods differ between these two subclasses. Butterworth, Chebyshev, elliptic filters are merely the result of different techniques (or optimizations) for designing IIR filters. These can achieve very good frequency responses with a small order, meaning you need few coefficients and therefore a small number of multiplications/additions per sample in the implementation. But if you really want a linear-phase response, say, then you'll need FIR filters. These have different design techniques and generally require higher orders for similar frequency characteristics. There are efficient implementations using Fast Fourier Transforms, but I doubt you would need such a thing.

我假设您可能会对具有有限脉冲响应 (FIR) 或无限脉冲响应 (IIR) 的线性时不变滤波器感兴趣。这两个子类的设计方法不同。Butterworth、Chebyshev、椭圆滤波器仅仅是设计 IIR 滤波器的不同技术（或优化）的结果。这些可以以很小的阶数实现非常好的频率响应，这意味着您需要很少的系数，因此在实现中每个样本需要少量的乘法/加法。但是，如果您真的想要线性相位响应，那么您将需要 FIR 滤波器。它们具有不同的设计技术，并且对于类似的频率特性通常需要更高的阶数。有使用快速傅立叶变换的有效实现，但我怀疑您是否需要这样的东西。

The possible different filter implementations differ mainly in numerical stability. You probably won't notice the difference unless you're using very low precision arithmetic and/or quite exotic coefficients. I believe the Matlab/Octave filter function uses the "Direct-form II" implementation, which is quite simple. You'll find a description on DSP books or the web, I'm sure.

可能的不同滤波器实现主要在数值稳定性方面有所不同。除非您使用非常低精度的算术和/或非常奇特的系数，否则您可能不会注意到差异。相信Matlab/Octave滤波器函数使用的是“Direct-form II”实现，非常简单。我敢肯定，您会在 DSP 书籍或网络上找到描述。

Jo?o Manuel Rodrigues

乔·曼努埃尔·罗德里格斯