Comment implémenter un filtre passe-bas en utilisant java

Question

J'essaie d'implémenter un filtre passe-bas en Java. Mon exigence est très simple, je dois éliminer les signaux au-delà d'une fréquence particulière (dimension unique). On dirait que le filtre Butterworth conviendrait à mon besoin.

Maintenant, l'important est que le temps CPU soit le plus bas possible. Il y aurait près d'un million d'échantillons que le filtre devrait traiter et nos utilisateurs n'aiment pas attendre trop longtemps. Existe-t-il une implémentation prête à l'emploi des filtres Butterworth qui a des algorithmes optimaux pour le filtrage.

Cordialement,

Chaitannya

Answer 1

Comme Mark Peters a dit dans son commentaire: Un filtre qui doit filtrer un lot doit être écrit en C ou C++. Mais vous pouvez toujours utiliser Java. Jetez un oeil à Java Native Interface (JNI). En raison de compilations C/C++ au code machine natif, il fonctionnera beaucoup plus vite que l'exécution de votre bytecode dans la machine virtuelle Java (JVM), qui est en fait un processeur virtuel qui traduit le bytecode à la machine locale sur l'ensemble d'instructions CPU comme x86, x64, ARM, ....)

Answer 2

J'ai récemment conçu une fonction butterworth simple (http://baumdevblog.blogspot.com/2010/11/butterworth-lowpass-filter-coefficients.html). Ils sont faciles à coder en Java et devraient être assez rapides si vous me le demandez (il suffirait de changer le filtre (double * samples, int count) pour filtrer (double [] samples, int count), je suppose). Le problème avec JNI est qu'il coûte l'indépendance de la plate-forme, peut confondre le compilateur de hotspot et les appels de méthode JNI dans votre code peuvent encore ralentir les choses. Donc je recommanderais d'essayer Java et de voir si c'est assez rapide. Dans certains cas, il peut être avantageux d'utiliser d'abord une transformée de Fourier rapide et d'appliquer le filtrage dans le domaine fréquentiel, mais je doute qu'il soit plus rapide qu'environ 6 multiplications et quelques ajouts par échantillon pour un simple filtre passe-bas.

Answer 3

La conception de filtres est un art de compromis, et pour bien le faire, vous devez prendre en compte certains détails.

Quelle est la fréquence maximale qui doit être dépassée "sans beaucoup" d'attention, et quelle est la valeur maximale de "sans trop"?

Quelle est la fréquence minimale qui doit être atténuée "beaucoup" et quelle est la valeur minimale de "beaucoup"?

Quelle est l'amplitude de l'ondulation (c'est-à-dire la variation de l'atténuation) dans les fréquences que le filtre est supposé passer?

Vous avez un large éventail de choix, ce qui vous coûtera une variété de montants de calcul. Un programme comme matlab ou scilab peut vous aider à comparer les compromis. Vous voudrez vous familiariser avec des concepts tels que l'expression des fréquences en tant que fraction décimale d'une fréquence d'échantillonnage, et l'échange entre les mesures linéaires et log (dB) de l'atténuation. Par exemple, un filtre passe-bas "parfait" est rectangulaire dans le domaine fréquentiel. Par exemple, un filtre passe-bas "parfait" est rectangulaire. Exprimé dans le domaine temporel comme une réponse impulsionnelle, ce serait une fonction sinc (sin x/x) avec les queues atteignant à la fois l'infini positif et négatif. Évidemment, vous ne pouvez pas calculer cela, alors la question devient si vous approximer la fonction sinc à une durée finie que vous pouvez calculer, combien cela dégrade votre filtre?Alternativement, si vous voulez un filtre à réponse impulsionnelle finie qui est très bon marché à calculer, vous pouvez utiliser un "box car" ou filtre rectangulaire où tous les coefficients sont 1. (Cela peut être fait encore moins cher si vous l'implémentez comme un filtre CIC exploitant le débordement binaire pour faire des accumulateurs 'circulaires', puisque vous prendrez le dérivé plus tard de toute façon). Mais un filtre qui est rectangulaire dans le temps ressemble à une fonction sin en fréquence - il a un sin x/x rolloff dans la bande passante (souvent élevé à une certaine puissance puisque vous auriez typiquement une version multi-étape), et certains "rebondissent" dans la bande d'arrêt. Dans certains cas, il est toujours utile, soit seul, soit lorsqu'il est suivi par un autre type de filtre.

Answer 4

J'ai une page décrivant un filtre passe-bas très simple, à très faible CPU, qui peut aussi être indépendant du framerate. Je l'utilise pour lisser l'entrée de l'utilisateur et aussi pour tracer des fréquences d'images souvent.

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

En bref, dans votre boucle de mise à jour:

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

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

Une valeur smoothing de 1 ne provoque pas de lissage de se produire, alors que des valeurs plus élevées à plus lisse le résultat.

La page a quelques fonctions écrites en JavaScript, mais la formule est indépendante du langage.

Answer 5

Voici un filtre passe-bas qui utilise une transformée de Fourier dans la bibliothèque mathématique 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

J'adopté ce billet depuis http://www.dspguide.com/ Je suis tout à fait nouveau java, donc ce n'est pas assez, mais il fonctionne

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


}