001: package JSci.maths.wavelet;
002:
003: import JSci.GlobalSettings;
004: import JSci.maths.*;
005: import JSci.maths.wavelet.*;
006: import JSci.maths.wavelet.splines.*;
007: import java.io.*;
008: import java.util.Arrays;
009:
010: /****************************************
011: * This class use the linear spline as a general model for a signal.
012: * While this is a reasonnable design choice, this can certainly be overwritten
013: * if necessary.
014: * Basic operations on signal are supported.
015: * @author Daniel Lemire
016: *****************************************/
017: public class Signal extends LinearSpline implements NumericalConstants,
018: Cloneable {
019: private Filter filterdual;
020: final static double normalisation = 1.0 / SQRT2;
021: private double[] param;
022:
023: /********************************************
024: * Return a copy of this object
025: *********************************************/
026: public Object clone() {
027: Signal s = (Signal) super .clone();
028: if (filterdual != null)
029: s.filterdual = filterdual;
030: //the previous line is unsafe?
031: if (param != null) {
032: s.param = new double[param.length];
033: System.arraycopy(s.param, 0, param, 0, param.length);
034: }
035: return (s);
036: }
037:
038: /************************************
039: * This method generates a copy of the
040: * current object with a different data
041: * content.
042: **************************************/
043: private Signal copy(double[] v) {
044: if (filterdual != null) {
045: if (param != null) {
046: double[] p = new double[param.length];
047: System.arraycopy(param, 0, p, 0, p.length);
048: return (new Signal(filterdual, v, p));
049: } else {
050: return (new Signal(filterdual, v));
051: }
052: } else {
053: return (new Signal(v));
054: }
055: }
056:
057: /******************************
058: *******************************/
059: public Signal() {
060: }
061:
062: /***********************************
063: ************************************/
064: public Signal(double[] v) {
065: super (v);
066: }
067:
068: /******************************
069: *******************************/
070: public Signal(Filter f, double[] v, double[] p) {
071: super (v);
072: filterdual = f;
073: System.arraycopy(p, 0, param, 0, param.length);
074: }
075:
076: /********************************
077: *********************************/
078: public Signal(Filter f) {
079: filterdual = f;
080: }
081:
082: /*********************************
083: **********************************/
084: public Signal(Filter f, double[] v) {
085: super (v);
086: filterdual = f;
087: }
088:
089: /********************************************
090: * Get the sampled values of the sample as
091: * an array.
092: *********************************************/
093: public double[] getValues() {
094: return (this .interpolate(0));
095: }
096:
097: /*********************************
098: * set the signal associated Filter
099: **********************************/
100: public void setFilter(Filter f) {
101: filterdual = f;
102: }
103:
104: /***********************************
105: * Set the parameter of the Filter (if
106: * it applies).
107: ************************************/
108: public void setParameters(double[] p) {
109: for (int k = 0; k < p.length; k++) {
110: param[k] = p[k];
111: }
112: }
113:
114: /***********************************
115: * Set the parameters of the Filter (if
116: * it applies).
117: ************************************/
118: public void setParameters(Double[] p) {
119: for (int k = 0; k < p.length; k++) {
120: param[k] = p[k].doubleValue();
121: }
122: }
123:
124: /***********************************
125: * Throws away the parameters of the Filter
126: ************************************/
127: public void removeParameters() {
128: param = null;
129: }
130:
131: /******************************************************
132: * Set the Signal to the specified length scraping or
133: * padding the beginning if necessary
134: *******************************************************/
135: public void setLengthFromEnd(int longueur) {
136: double[] newvec = ArrayMath.setLengthFromEnd(this .evaluate(0),
137: longueur);
138: this .setValues(newvec);
139: }
140:
141: /*****************************************************
142: * Resample the signal using linear interpolation
143: ******************************************************/
144: public void resample(int newl) {
145: double[] newvec = EngineerMath.resample(this .evaluate(0), newl);
146: this .setValues(newvec);
147: }
148:
149: /******************************************************
150: * Set the Signal to the specified length scraping or
151: * padding the end if necessary
152: *******************************************************/
153: public void setLengthFromBeginning(int longueur) {
154: double[] newvec = ArrayMath.setLengthFromBeginning(this
155: .evaluate(0), longueur);
156: this .setValues(newvec);
157: }
158:
159: /***********************************
160: * Set the data for the signal
161: ************************************/
162: public void setData(double[] v) {
163: setValues(v);
164: }
165:
166: /********************************
167: * Fast Wavelet Transform
168: *********************************/
169: public FWTCoef fwt(int J) {
170: if (J > 20) {
171: throw new IllegalArgumentException("Too many iterations.");
172: }
173: if (J < 0) {
174: throw new IllegalArgumentException(
175: "Cannot have a negative number of iterations.");
176: }
177: double[] data = this .interpolate(0);
178: double[][] fwt = new double[J + 1][];
179: for (int j = 1; j <= J; j++) {
180: fwt[j] = highpassProject(data);
181: data = lowpassProject(data);
182: }
183: fwt[0] = data;
184: FWTCoef t = new FWTCoef(fwt);
185: return (t);
186: }
187:
188: /********************************
189: * The Fast Wavelet Transform
190: * with Wavelet packets
191: * @param J number of iterations
192: * @param cout cost function
193: *********************************/
194: public FWTPacketCoef fwtPacket(int J, MappingND cout) {
195: if (J > 20) {
196: throw new IllegalArgumentException("Too many iterations.");
197: }
198: if (J < 0) {
199: throw new IllegalArgumentException(
200: "Cannot have a negative number of iterations.");
201: }
202: double[] data = this .interpolate(0);
203: double[][] fwt = new double[J + 1][];
204: double[] choix1, choix2;
205: boolean[] choixStandard = new boolean[J];
206: for (int j = 0; j < J; j++) {
207: choix1 = highpassProject(data);
208: choix2 = lowpassProject(data);
209: if (cout.map(choix1)[0] < cout.map(choix2)[0]) {
210: fwt[j] = choix1;
211: data = choix2;
212: choixStandard[j] = true;
213: } else {
214: data = choix1;
215: fwt[j] = choix2;
216: choixStandard[j] = false;
217: }
218: }
219: fwt[J] = data;
220: FWTPacketCoef t = new FWTPacketCoef(fwt, choixStandard);
221: return (t);
222: }
223:
224: /************************************************
225: * Project the array according to the lowpass Filter
226: * @param v data array
227: * @author Daniel Lemire
228: ************************************************/
229: private double[] lowpassProject(double[] v) {
230: int l = filterdual.previousDimension(v.length);
231: double[] Eche;
232: double[] ans = new double[l];
233: if (param != null) {
234: for (int i = 0; i < l; i++) {
235: double[] p = new double[param.length];
236: System.arraycopy(param, 0, p, 0, p.length);
237: Eche = filterdual
238: .lowpass(delta(i, l, normalisation), p);
239: ans[i] = scalarProduct(v, Eche);
240: }
241: } else {
242: for (int i = 0; i < l; i++) {
243: Eche = filterdual.lowpass(delta(i, l, normalisation));
244: ans[i] = scalarProduct(v, Eche);
245: }
246: }
247:
248: /*
249: for(int i=0;i<l;i++) {
250: if (param!=null) {
251: double[] p=new double[param.length];
252: System.arraycopy(param,0,p,0,p.length);
253: Eche=filterdual.lowpass(delta(i,l,normalisation),p);
254: } else {
255: Eche=filterdual.lowpass(delta(i,l,normalisation));
256: }
257: ans[i]=ArrayMath.scalarProduct(v,Eche);
258: }
259: */
260: return (ans);
261: }
262:
263: private static double scalarProduct(double[] w0, double[] w1) {
264: double sortie = 0.0;
265: for (int k = 0; k < w0.length; k++) {
266: sortie += w0[k] * w1[k];
267: }
268: return (sortie);
269: }
270:
271: /****************************************************
272: * Project the data according to the lowpass Filter
273: * @author Daniel Lemire
274: ****************************************************/
275: public double[] lowpassProject() {
276: double[] data = this .interpolate(0);
277: return (lowpassProject(data));
278: }
279:
280: /******************************************
281: * Project the signal according the the
282: * highpass Filter
283: * @param v data
284: * @author Daniel Lemire
285: ********************************************/
286: private double[] highpassProject(double[] v) {
287: int l = filterdual.previousDimension(v.length);
288: int lOnd = v.length - l;
289: double[] Onde;
290: double[] ans = new double[lOnd];
291: if (param != null) {
292: for (int i = 0; i < lOnd; i++) {
293: double[] p = new double[param.length];
294: System.arraycopy(param, 0, p, 0, p.length);
295: Onde = filterdual.highpass(delta(i, lOnd, 1), p);
296: ans[i] = scalarProduct(v, Onde);
297: }
298: } else {
299:
300: for (int i = 0; i < lOnd; i++) {
301:
302: Onde = filterdual.highpass(delta(i, lOnd, 1));
303:
304: ans[i] = scalarProduct(v, Onde);
305: }
306:
307: }
308:
309: /*
310: for(int i=0;i<lOnd;i++) {
311: if (param!=null) {
312: double[] p=new double[param.length];
313: System.arraycopy(param,0,p,0,p.length);
314: Onde=filterdual.highpass(delta(i,lOnd,1),p);
315: } else {
316: Onde=filterdual.highpass(delta(i,lOnd,1));
317: }
318: ans[i]=ArrayMath.scalarProduct(v,Onde);
319: }
320: */
321: return (ans);
322: }
323:
324: /********************************************
325: * Project the signal according the the
326: * highpass Filter
327: * @author Daniel Lemire
328: *********************************************/
329: public double[] highpassProject() {
330: double[] data = this .interpolate(0);
331: return (highpassProject(data));
332: }
333:
334: /*************************************************
335: * return a kronecker
336: * @param l length
337: * @param i position
338: * @author Daniel Lemire
339: *************************************************/
340: private double[] delta(int i, int l, double a) {
341: if ((i < 0) || (i > l) || (l < 0)) {
342: throw new IllegalArgumentException(
343: "This Kronecker doesn't exist.");
344: }
345: double[] v = new double[l];
346: v[i] = a;
347: return (v);
348: }
349:
350: /***************************
351: * Compute the L2 norm of the
352: * signal
353: ****************************/
354: public double norm() {
355: double[] data = this .interpolate(0);
356: return (ArrayMath.norm(data));
357: }
358:
359: /***********************************
360: * @author Don Cross
361: * @author Daniel Lemire
362: ************************************/
363: public Complex[] fft() {
364: double[] data = this .interpolate(0);
365: return (fft(data));
366: }
367:
368: /**
369: * Performs the Fourier transform.
370: * Convenience method for {@link JSci.maths.FourierMath#transform(double[]) FourierMath.transform}.
371: */
372: public static Complex[] fft(double[] data) {
373: return FourierMath.transform(data);
374: }
375:
376: /**
377: * Performs the Fourier transform.
378: * Convenience method for {@link JSci.maths.FourierMath#transform(Complex[]) FourierMath.transform}.
379: */
380: public static Complex[] fft(Complex[] data) {
381: return FourierMath.transform(data);
382: }
383:
384: /*********************************
385: * Return the absolute value of
386: * the FFT
387: **********************************/
388: public double[] absFFT() {
389: Complex[] fft = fft();
390: double[] answer = new double[fft.length];
391: for (int i = 0; i < fft.length; i++) {
392: answer[i] = fft[i].mod();
393: }
394: return (answer);
395: }
396:
397: public static double[] absFFT(double[] data) {
398: Complex[] fft = fft(data);
399: double[] answer = new double[fft.length];
400: for (int i = 0; i < fft.length; i++) {
401: answer[i] = fft[i].mod();
402: }
403: return (answer);
404: }
405:
406: /**************************************
407: * Also noted iFFT in other packages.
408: * This is the inverse to the FFT.
409: ***************************************/
410: public static Complex[] fftInverse(Complex data[]) {
411: return FourierMath.inverseTransform(data);
412: }
413:
414: /*****************************************
415: * Check if another object is equal to this
416: * Signal object
417: ******************************************/
418: public boolean equals(Signal b) {
419: return (Arrays.equals(this .getValues(), b.getValues()));
420: }
421:
422: /*******************************************************
423: * Will make the signal a given dimension
424: ********************************************************/
425: public void setDimensionFromEnd(int dimension) {
426: double[] data = this .interpolate(0);
427: double[] ans = new double[dimension];
428: int debut;
429: if (dimension - data.length < 0)
430: debut = data.length - dimension;
431: else
432: debut = 0;
433: for (int k = debut; k < data.length; k++) {
434: ans[k + dimension - data.length] = data[k];
435: }
436: super .setValues(ans);
437: }
438:
439: /********************************************************
440: * Will make the signal a given dimension
441: *********************************************************/
442: public void setDimensionFromBeginning(int dimension) {
443: double[] data = this .interpolate(0);
444: double[] ans = new double[dimension];
445: int debut;
446: if (dimension - data.length < 0)
447: debut = data.length - dimension;
448: else
449: debut = 0;
450: for (int k = 0; k < data.length - debut; k++) {
451: ans[k] = data[k];
452: }
453: super .setValues(ans);
454: }
455:
456: /********************************************
457: * Simplistic FFT denoising.
458: * @param k frequency to denoised
459: *********************************************/
460: public void denoiseByFFT(int k) {
461: if (k < 1) {
462: throw new IllegalArgumentException(
463: "This parameter must be 1 or more : " + k);
464: }
465: double[] data = interpolate(0);
466: if (k > data.length - 2) {
467: if (data.length < 4) {
468: throw new IllegalArgumentException(
469: "Your signal is too short to be denoised : "
470: + data.length + " < 4");
471: }
472: throw new IllegalArgumentException(
473: "Since you signal has dimension " + data.length
474: + ", the parameter must be at most : "
475: + (data.length - 2));
476: }
477: Complex[] ff = Signal.fft(data);
478: ff[k + 1] = Complex.ZERO;
479: ff[data.length - 1 - k] = Complex.ZERO;
480: Complex[] tf = Signal.fftInverse(ff);
481: for (int l = 0; l < data.length; l++) {
482: data[l] = tf[l].real();
483: if (Math.abs(tf[l].imag()) > GlobalSettings.ZERO_TOL) {
484: throw new IllegalArgumentException(
485: "Complex values detected during synthesis. Please get in touch with Daniel Lemire at Daniel.Lemire@Tintin.net to report this error.");
486: }
487: }
488: super .setValues(data);
489: }
490:
491: /**********************************************
492: * Return the entropy of the signal
493: ***********************************************/
494: public double entropy() {
495: return (EngineerMath.entropy(this .evaluate(0)));
496: }
497:
498: /*****************************************
499: * Apply the given array as a convolution
500: * Filter and return a new Signal.
501: * As one often want to compare the result
502: * to the original signal, this method is
503: * "safe", that is, it won't change the current
504: * object.
505: * @param f an array containing the coefficients
506: * of the convolution Filter
507: ******************************************/
508: public Signal filter(double[] f) {
509: double[] data = interpolate(0);
510: if (data.length - (f.length - 1) <= 0) {
511: throw new IllegalArgumentException(
512: "Your signal is too short for this Filter : "
513: + data.length + ", " + f.length);
514: }
515: double[] ans = new double[data.length - (f.length - 1)];
516: for (int k = 0; k < data.length - (f.length - 1); k++) {
517: for (int l = 0; l < f.length; l++) {
518: ans[k] += f[l] * data[k + l];
519: }
520: }
521: return (copy(ans));
522: }
523:
524: /******************************************
525: * Apply the median Filter of a window of
526: * size 2*n+1.
527: * exception IllegalArgumentException if the
528: * parameter n is negative
529: ********************************************/
530: public Signal medianFilter(int n) {
531: if (n < 0)
532: throw new IllegalArgumentException(
533: "The parameter must be positive: " + n + " < 0");
534: double[] data = super .interpolate(0);
535: if (data.length - 2 * n <= 0) {
536: throw new IllegalArgumentException(
537: "Your signal is too short for this Filter : "
538: + data.length + " - " + (2 * n) + " = "
539: + (data.length - 2 * n));
540: }
541: double[] ans = new double[data.length - 2 * n];
542: double[] vtemp = new double[2 * n + 1];
543: for (int k = 0; k < data.length - 2 * n; k++) {
544: for (int l = 0; l < 2 * n + 1; l++) {
545: vtemp[l] = data[k + l];
546: }
547: ans[k] = ArrayMath.median(vtemp);
548: }
549: return (copy(ans));
550:
551: }
552:
553: /****************************************************
554: * This denoising method will identify
555: * "short peaks" in the signal and take them away.
556: * Short peaks are defined from a comparison
557: * with the median filtered signal.
558: * Only "significative" peaks are detected (see parameter
559: * p).
560: * This method won't denoise near the boundaries.
561: * "Short" refers here to the time-domain and
562: * not the amplitude.
563: * param p percentage of the range (max-min) considered
564: * as a significative step
565: * param n length of the peak in the time domain
566: * exception IllegalArgumentException if p is not between
567: * 0 and 1
568: * exception IllegalArgumentException if the
569: * parameter n is negative
570: *******************************************************/
571: public Signal denoiseShortPeaks(double p, int n) {
572: if ((p < 0) || (p > 1)) {
573: throw new IllegalArgumentException(
574: "The parameter p must be between 0 and 1: " + p);
575: }
576: double[] values = this .interpolate(0);
577: double range = ArrayMath.max(values) - ArrayMath.min(values);
578: double threshold = range * p;
579: double[] med = (this .medianFilter(n)).interpolate(0);
580: for (int k = n; k < values.length - n; k++) {
581: if (Math.abs(values[k] - med[k - n]) > threshold) {
582: values[k] = med[k - n];
583: }
584: }
585: return (copy(values));
586: }
587: }
|