How to implement low pass filter using java

Question

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.

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.

Benefice

Audacity is open source and contains many audio filters. They will be written in C/C++, but that's pretty simple to translate into equivalent Java code. – Wallford

Maybe you could show some code so that we know what you are trying to filter ? – Chevaldefrise

I have a tutorial here that includes second order Butterworth filters. It should be easy to implement this in Java: blog.bjornroche.com/2012/08/basic-audio-eqs.html – Barina

Answer 1

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.

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 smoothing value of 1 causes no smoothing to occur, while higher values increasingly smooth out the result.

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

Curlpaper

Love the algorithm, but is there a way to convert the smoothing value to a cutoff frequency? Thanks – Arduous

Answer 2

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

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

Colbycolbye

Answer 3

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.

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?

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.

Other

Answer 4

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).

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.

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.

Buskirk

I would NOT try to filter a million data points (as the OP suggested) using Fourier methods: blog.bjornroche.com/2012/08/why-eq-is-done-in-time-domain.html – Barina

Answer 5

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, ....)

Woofer

Before you rewrite - benchmark, you would be surprised that the difference is not as great as you think. In many instances Java is actually faster than C/C++. – Chevaldefrise

Answer 6

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

/*
* 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;
    }


}

Amp

What's the filter order? @2sloth – Billfish

Filter order is one of the parameters for the method, so whatever you want :) – Amp

Answer 7

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.

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.

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.

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.

João Manuel Rodrigues

Sherrell