Jak zaimplementować filtr dolnoprzepustowy za pomocą java

Question

Próbuję zaimplementować filtr dolnoprzepustowy w Javie. Moje wymaganie jest bardzo proste, muszę wyeliminować sygnały poza określoną częstotliwość (pojedynczy wymiar). Wygląda na to, że filtr Butterworth pasowałby do moich potrzeb.

Ważne jest teraz, aby czas procesora był jak najniższy. Będzie blisko miliona próbek, które filtr będzie musiał przetworzyć, a nasi użytkownicy nie lubią zbyt długo czekać. Czy istnieje gotowa implementacja filtrów Butterwortha, która ma optymalne algorytmy filtrowania.

Pozdrawiam,

Chaitannya

Answer 1

Podobnie jak Mark Peters powiedział w swoim komentarzu: filtr, który musi filtrować dużo powinny być napisane w języku C lub C++. Ale nadal możesz korzystać z Javy. Wystarczy spojrzeć na Java Native Interface (JNI). Z powodu kompilacji C/C++ do natywnego kodu maszynowego, będzie działał o wiele szybciej niż działanie kodu bajtowego w wirtualnej maszynie Java (JVM), która jest w rzeczywistości wirtualnym procesorem, który tłumaczy kod bajtowy na lokalny komputer jego natywnym kodem (w zależności od na zestaw instrukcji procesora, takich jak x86, x64, ARM, ....)

Answer 2

Ostatnio zaprojektowałem prostą funkcję butterwortha (http://baumdevblog.blogspot.com/2010/11/butterworth-lowpass-filter-coefficients.html). Są łatwe do kodowania w Javie i powinny być wystarczająco szybkie, jeśli mnie pytasz (musisz tylko zmienić filtr (double * sample, int count) do filtrowania (double [] sample, int count), jak sądzę).

Problem z JNI polega na tym, że kosztuje niezależność od platformy, może zmylić kompilator hotspotów, a wywołania metody JNI w ramach kodu mogą jeszcze spowolnić działanie. Polecam wypróbować Javę i sprawdzić, czy jest wystarczająco szybki.

W niektórych przypadkach korzystna może być szybka transformata Fouriera i zastosowanie filtrowania w dziedzinie częstotliwości, ale wątpię, aby było to szybciej niż około 6 mnożeń i kilka dodatków na próbkę dla prostego filtra dolnoprzepustowego.

Answer 3

Projektowanie filtrów jest sztuką kompromisów, a aby to zrobić dobrze, należy wziąć pod uwagę kilka szczegółów.

Jaka jest maksymalna częstotliwość, którą należy podać "bez większego" podkreślenia, a jaka jest maksymalna wartość "bez dużo"?

Jaka jest minimalna częstotliwość, która musi być tłumiona "dużo" i jaka jest minimalna wartość "dużo"?

Ile tętnień (np. Zmiana tłumienia) jest akceptowalna w zakresie częstotliwości, które filtr ma przepuścić?

Masz szeroki wybór opcji, które będą Cię kosztowały różne ilości obliczeń. Program podobny do Matlab lub Scilab może pomóc w porównaniu kompromisów. Będziesz chciał zapoznać się z pojęciami takimi jak wyrażanie częstotliwości jako ułamka dziesiętnego częstotliwości próbkowania i zamiana między pomiarem tłumienia liniowej i logarytmicznej (dB).

Na przykład "doskonały" filtr dolnoprzepustowy jest prostokątny w dziedzinie częstotliwości. Wyrażony w dziedzinie czasu jako odpowiedź impulsowa, byłaby to funkcja sinc (sin x/x), której ogony osiągałyby zarówno dodatnią, jak i ujemną nieskończoność. Oczywiście nie można tego obliczyć, więc pojawia się pytanie, czy przybliżasz funkcję sinc do skończonego czasu, który możesz obliczyć, ile to obniża twój filtr?

Alternatywnie, jeśli potrzebujesz skończonego filtra odpowiedzi impulsowej, który jest bardzo tani do obliczenia, możesz użyć "wózka skrzyniowego" lub filtra prostokątnego, w którym wszystkie współczynniki wynoszą 1. (To może być nawet tańsze, jeśli je wdrożysz jako filtr CIC wykorzystujący przepełnienie binarne do akumulatorów "okrągłych", ponieważ i tak będziesz używał później pochodnej). Ale filtr, który jest prostokątny w czasie, wygląda jak funkcja sinc z częstotliwością - ma przesunięcie sin x/x w pasmie przenoszenia (często podniesione do pewnej mocy, ponieważ zazwyczaj masz wersję wieloetapową), a niektóre "odbijają się" w pasie zatrzymania. W niektórych przypadkach jest to użyteczne, samo lub w przypadku zastosowania innego rodzaju filtra.

Answer 4

Mam stronę opisującą bardzo prosty, bardzo niski procesor filtru dolnoprzepustowego, który również może być niezależny od liczby klatek na sekundę. Używam go do wygładzania danych wprowadzanych przez użytkownika, a także do częstego obliczania częstości klatek.

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

W skrócie, w pętli aktualizacji:

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

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

smoothing wartość 1 przyczyn nie wygładzania występuje, natomiast wyższe wartości coraz wygładzić wynik.

Strona zawiera kilka funkcji napisanych w języku JavaScript, ale formuła jest agnostyczna.

Answer 5

Oto filtr dolnoprzepustowy, który korzysta z transformaty Fouriera w bibliotece apache matematyki.

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

przyjąłem to z http://www.dspguide.com/ Jestem zupełnie nowy w Javie, więc to nie jest ładna, ale działa

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


}