Come implementare il filtro passa-basso usando java

Question

Sto cercando di implementare un filtro passa-basso in Java. Il mio requisito è molto semplice, devo eliminare i segnali oltre una particolare frequenza (dimensione singola). Sembra che il filtro Butterworth soddisfi le mie esigenze.

Ora la cosa importante è che il tempo di CPU dovrebbe essere il più basso possibile. Ci sarebbe quasi un milione di campioni che il filtro dovrebbe elaborare e ai nostri utenti non piace aspettare troppo a lungo. Ci sono implementazioni readymade di filtri Butterworth con algoritmi ottimali per il filtraggio.

saluti,

Chaitannya

Answer 1

come Mark Peters ha detto nel suo commento: Un filtro che deve filtrare un sacco dovrebbe essere scritto in C o C++. Ma puoi ancora fare uso di Java. Basta dare un'occhiata a Java Native Interface (JNI). Grazie alla compilazione di C/C++ al codice macchina nativo, verrà eseguito molto più rapidamente rispetto all'esecuzione del bytecode nella JVM (Java Virtual Machine), che è in realtà un processore virtuale che traduce il bytecode sulla macchina locale il suo codice nativo (in base sull'istruzione della CPU impostata come x86, x64, ARM, ....)

Answer 2

Ho progettato una funzione Butterworth semplice (http://baumdevblog.blogspot.com/2010/11/butterworth-lowpass-filter-coefficients.html). Sono facili da codificare in Java e dovrebbero essere abbastanza veloci se me lo chiedi (dovresti semplicemente cambiare filtro (doppio * sample, int count) per filtrare (double [] samples, int count), suppongo).

Il problema con JNI è che costa l'indipendenza dalla piattaforma, potrebbe confondere il compilatore di hotspot e le chiamate al metodo JNI all'interno del codice potrebbero comunque rallentare. Quindi consiglierei di provare Java e vedere se è abbastanza veloce.

In alcuni casi potrebbe essere utile utilizzare prima una trasformata di Fourier veloce e applicare il filtro nel dominio della frequenza, ma dubito che questo sia più veloce di circa 6 multipli e alcune aggiunte per campione per un semplice filtro passa-basso.

Answer 3

La progettazione del filtro è un'arte di compromessi, e per farlo bene è necessario prendere in considerazione alcuni dettagli.

Qual è la frequenza massima che deve essere passata "senza molto" attenzione, e qual è il valore massimo di "senza molto"?

Qual è la frequenza minima che deve essere attenuata "molto" e qual è il valore minimo di "molto"?

Quanta ondulazione (cioè variazione nell'attenuazione) è accettabile all'interno delle frequenze che il filtro dovrebbe passare?

Hai una vasta gamma di scelte, che ti costeranno una varietà di quantità di calcolo. Un programma come matlab o scilab può aiutarti a confrontare i compromessi. Ti consigliamo di acquisire familiarità con concetti come l'espressione delle frequenze come frazione decimale di una frequenza di campionamento e lo scambio tra le misure di attenuazione lineare e log (dB).

Ad esempio, un filtro passa basso "perfetto" è rettangolare nel dominio della frequenza. Espresso nel dominio del tempo come risposta all'impulso, sarebbe una funzione sinc (sin x/x) con la coda che raggiunge l'infinito sia positivo che negativo. Ovviamente non è possibile calcolarlo, quindi la domanda diventa se si approssima la funzione sinc a una durata finita che è possibile calcolare, quanto può degradare il filtro?

In alternativa, se si desidera un filtro di risposta impulsivo finito che è molto economico da calcolare, è possibile utilizzare un "box car" o un filtro rettangolare in cui tutti i coefficienti sono 1. (Questo può essere reso ancora più economico se lo si implementa come un filtro CIC che sfrutta l'overflow binario per fare accumulatori "circolari", dato che in seguito si utilizzerà comunque la derivata). Ma un filtro che è rettangolare nel tempo assomiglia ad una funzione sinc in frequenza - ha un sin x/x rolloff nella banda passante (spesso elevato a qualche potenza dato che di solito si dispone di una versione multistadio), e alcuni "rimbalzano" nella banda di stop. Ancora in alcuni casi è utile, da solo o quando seguito da un altro tipo di filtro.

Answer 4

Ho una pagina che descrive un filtro passa-basso molto semplice, a bassissima CPU che è anche in grado di essere indipendente dal framerate. Lo uso per appianare l'input dell'utente e anche per rappresentare graficamente i frame rate.

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

In breve, nel loop aggiornamento:

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

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

Un valore smoothing di 1 cause non lisciatura a verificarsi, mentre valori più alti sempre appianare il risultato.

La pagina ha un paio di funzioni scritte in JavaScript, ma la formula è indipendente dal linguaggio.

Answer 5

Ecco un filtro passa-basso che utilizza una trasformata di Fourier nella libreria matematica 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

ho adottato questo da http://www.dspguide.com/ Sono abbastanza nuovo a Java, quindi non è abbastanza, ma funziona

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


}