Voici une implémentation rudimentaire d'un filtrage passe-bas sur un signal 1D. Ce code fait suite à une discussion du forum "Algorithmes".
Résultat du filtrage d'un signal carré bruité par la méthode #2 avec le filtre de Butterworth.
- Calcul de la DFT (Discret Fourier Transform) et de son inverse, par la méthode brutale:
Code java : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22 /** * Compute the DFT of real data contained in input[] * The result is stored in rDFT[] and iDFT[] * * @param input array of size N containing the data * @param rDFT empty array of size N/2 (real part of the DFT) * @param iDFT empty array of size N/2 (imaginary part of the DFT) */ public static void DFT(double[] input, double[] rDFT, double[] iDFT) { int N = input.length; for (int f = 0; f < N/2; f++) { rDFT[f] = 0; iDFT[f] = 0; for (int i = 0; i < N; i++) { double w = 2 * Math.PI * (double) i / N; rDFT[f] += input[i] * Math.cos(f * w); iDFT[f] -= input[i] * Math.sin(f * w); } rDFT[f] /= N/2; iDFT[f] /= N/2; } }
Code java : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 /** * Compute the inverse DFT of complex data contained in rDFT[] and iDFT[] * The result is stored in output[] * * @param rDFT array of size N/2 containing the real part of the data * @param iDFT array of size N/2 containing the real part of the data * @param output empty array of size N */ public static void invDFT(double[] rDFT, double[] iDFT, double[] output) { int N = output.length; for (int i = 0; i < N; i++) { output[i] = 0; for (int f = 0; f < N/2; f++) { double w = 2 * Math.PI * (double) i / N; output[i] += rDFT[f] * Math.cos(f * w) - iDFT[f] * Math.sin(f * w); } } }
- Calcul des valeurs complexes d'un filtre de circuit RC:
Code java : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15 double[] ReFil = new double[N/2]; double[] ImFil = new double[N/2]; double rc = 1.0/40; for(int f=0;f<(N/2);f++) { double w = 2*Math.PI*f; double wrc = w*rc; double gain = 1.0/(1+Math.pow(wrc,2)); double real = 1; double imag = -wrc; ReFil[f]=real*gain; ImFil[f]=imag*gain; }
- Calcul des valeurs complexes d'un filtre de Butterworth (ordre 2):
Code java : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15 double[] ReFil = new double[N/2]; double[] ImFil = new double[N/2]; double c = 100.0; // cutoff frequency for (int f = 0; f < N / 2; f++) { double w = 2 * Math.PI * f; double wc = w/c; double gain = 1.0/(1+Math.pow(wc,4)); double real = 1-Math.pow(wc,2); double imag = -Math.sqrt(2)*wc; ReFil[f] = real * gain; ImFil[f] = imag * gain; }
- Méthode de Filtrage 1 : DFT(signal) --> output = produit complexe filtre*signal --> invDFT(output)
Code java : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24 int N = size_of_your_signal; double[] signal = new double[N]; // ... fill signal[] with data here ... double[] ReFil = new double[N/2]; double[] ImFil = new double[N/2]; // ... compute filter here (see above) ... // DFT of the signal double[] ReSig = new double[N/2]; double[] ImSig = new double[N/2]; DFT(signal, ReSig, ImSig); // output = filter * signal (warning: multiplication of two complex) double[] ReOutput = new double[N/2]; double[] ImOutput = new double[N/2]; for(int f=0;f<N/2;f++) { ReOutput[f]=ReSig[f]*ReFil[f]-ImSig[f]*ImFil[f]; ImOutput[f]=ReSig[f]*ImFil[f]+ImSig[f]*ReFil[f]; } // inverse DFT of the output double[] output = new double[N]; invDFT(ReOutput, ImOutput, output);
- Méthode de Filtrage 2 : invDFT(filtre) -> réduction du filtre -> output = convolution filtre*signal
Code java : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28 int N = size_of_your_signal; double[] signal = new double[N]; // ... fill signal[] with data here ... double[] ReFil = new double[N/2]; double[] ImFil = new double[N/2]; // ... compute filter here ... // inverse DFT of filter double[] fullfilter = new double[N]; invDFT(ReFil, ImFil, fullfilter); // reduce filter size (50 elements) and normalize double[] filter = new double[50]; double filternorm = 0; for (int f = 0; f < filter.length; f++) filternorm += fullfilter[f]; for (int f = 0; f < filter.length; f++) filter[f] = fullfilter[f] / filternorm; // convolution filter*signal double[] output = new double[N]; for (int i = 0; i < N; i++) { output[i] = 0; for (int f = 0; f < filter.length; f++) output[i] += signal[(N + i - f) % N] * filter[f]; }
Partager