[Java] Filtrage passe-bas d'un signal 1D

Voici une implémentation rudimentaire d'un filtrage passe-bas sur un signal 1D. Ce code fait suite à une discussion du forum "Algorithmes".

http://xphilipp.developpez.com/contribuez/lowpass.png
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:

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:

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:

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:

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:

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:

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];
}