Как реализовать фильтр низких частот с использованием Java

Question

Я пытаюсь реализовать фильтр нижних частот в Java.Мое требование очень простое, я должен устранить сигналы за пределами определенной частоты (одно измерение).Похоже, фильтр Баттерворта подойдет мне.

Теперь важно, чтобы время процессора было как можно меньше.Было бы около миллиона выборок, которые фильтр должен обработать, и наши пользователи не любят ждать слишком долго.Есть ли готовая реализация фильтров Баттерворта, которая имеет оптимальные алгоритмы фильтрации.

Answer 1

У меня есть страница, описывающая очень простой фильтр низких частот с очень низкой загрузкой ЦП, который также может быть независимым от частоты кадров. Я использую его для сглаживания пользовательского ввода, а также для отображения частоты кадров.

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

A 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

Вот фильтр нижних частот, который использует преобразование Фурье в математической библиотеке 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 4

Дизайн фильтра - это искусство компромиссов, и для его успешного выполнения необходимо принять во внимание некоторые детали.

Какова максимальная частота, которую необходимо передать "без особого" ослабления, и каковамаксимальное значение «без много»?

Какова минимальная частота, которая должна быть ослаблена «много» и каково минимальное значение «много»?

Сколько пульсаций (т.е.вариация ослабления) приемлема в пределах частот, которые должен пропускать фильтр?

У вас есть широкий выбор вариантов, которые обойдутся вам в различные объемы вычислений. Такая программа, как matlab или scilab, может помочь вам сравнить компромиссы .Вы захотите ознакомиться с такими понятиями, как выражение частот в виде десятичной дроби частоты дискретизации и взаимозаменяемость между линейными и логарифмическими (дБ) измерениями ослабления.

Например, «идеальный» фильтр нижних частотявляется прямоугольным в частотной области.Выраженная во временной области как импульсная характеристика, это будет функция синуса (sin x / x) с хвостами, достигающими как положительной, так и отрицательной бесконечности.Очевидно, что вы не можете рассчитать это, поэтому возникает вопрос, приближаете ли вы функцию sinc к конечной продолжительности, которую вы можете вычислить, насколько это ухудшает ваш фильтр?

Альтернативно, если вы хотите получить конечную импульсную характеристикуФильтр, который очень дешев для расчета, вы можете использовать «коробочную машину» или прямоугольный фильтр, где все коэффициенты равны 1. (Это можно сделать еще дешевле, если вы внедрите его как фильтр CIC, использующий двоичное переполнение для создания «круговых» аккумуляторов, так как вы будете принимать производную позже в любом случае).Но прямоугольный по времени фильтр выглядит как sinc-функция по частоте - он имеет спад sin x / x в полосе пропускания (часто повышенный до некоторой мощности, поскольку у вас обычно будет многоступенчатая версия), а некоторые «возвращаются в норму»в стоп-группе.Тем не менее, в некоторых случаях это полезно, либо само по себе, либо после использования фильтра другого типа.

Answer 5

Недавно я разработал простую функцию Баттерворта (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 6

Как сказал Марк Петерс в своем комментарии: фильтр, который должен фильтровать много, должен быть написан на C или C ++. Но вы все равно можете использовать Java. Просто взгляните на Собственный интерфейс Java (JNI) . Из-за того, что C / C ++ компилируется в собственный машинный код, он будет работать намного быстрее, чем запуск вашего байт-кода на виртуальной машине Java (JVM), которая на самом деле является виртуальным процессором, который переводит байт-код на локальную машину своим собственным кодом (в зависимости от на ЦП, например, x86, x64, ARM , ....)