Signal


1   package JSci.maths.wavelet;
2   
3   import JSci.GlobalSettings;
4   import JSci.maths.*;
5   import JSci.maths.wavelet.*;
6   import JSci.maths.wavelet.splines.*;
7   import java.io.*;
8   import java.util.Arrays  ;
9   
10  /****************************************
11  * This class use the linear spline as a general model for a signal.
12  * While this is a reasonnable design choice, this can certainly be overwritten
13  * if necessary.
14  * Basic operations on signal are supported.
15  * @author Daniel Lemire
16  *****************************************/
17  public class Signal extends LinearSpline implements NumericalConstants, Cloneable    {
18      private Filter filterdual;
19      final static double normalisation=1.0/SQRT2;
20      private double[] param;
21    /********************************************
22    * Return a copy of this object
23    *********************************************/
24      public Object   clone() {
25      Signal s=(Signal) super.clone();
26      if(filterdual!=null)
27        s.filterdual=filterdual;
28        //the previous line is unsafe?
29      if(param!=null) {
30        s.param=new double[param.length];
31        System.arraycopy(s.param,0,param,0,param.length);
32      }
33      return(s);
34      }
35    /************************************
36    * This method generates a copy of the
37    * current object with a different data
38    * content.
39    **************************************/
40    private Signal copy(double[] v) {
41          if(filterdual!=null) {
42              if(param!=null) {
43                  double[] p=new double[param.length];
44                                  System.arraycopy(param,0,p,0,p.length);
45                  return(new Signal(filterdual,v,p));
46              } else {
47                  return(new Signal(filterdual,v));
48              }
49          } else {
50              return(new Signal(v));
51          }
52    }
53      /******************************
54      *******************************/
55      public Signal () {
56    }
57  
58      /***********************************
59      ************************************/
60      public Signal (double[] v) {
61          super(v);
62      }
63  
64      /******************************
65      *******************************/
66      public Signal (Filter f,double[] v, double[] p) {
67          super(v);
68          filterdual=f;
69                  System.arraycopy(p,0,param,0,param.length);
70      }
71  
72      /********************************
73      *********************************/
74      public Signal(Filter f) {
75          filterdual=f;
76      }
77  
78      /*********************************
79      **********************************/
80      public Signal (Filter f,double[] v) {
81          super(v);
82          filterdual=f;
83      }
84  
85  
86      /********************************************
87    * Get the sampled values of the sample as
88    * an array.
89      *********************************************/
90      public double[] getValues() {
91          return(this.interpolate(0));
92      }
93  
94      /*********************************
95    * set the signal associated Filter
96      **********************************/
97      public void setFilter(Filter f) {
98          filterdual=f;
99      }
100 
101     /***********************************
102   * Set the parameter of the Filter (if
103   * it applies).
104     ************************************/
105     public void setParameters(double[] p) {
106         for(int k=0;k<p.length;k++) {
107             param[k]=p[k];
108         }
109     }
110 
111     /***********************************
112   * Set the parameters of the Filter (if
113   * it applies).
114     ************************************/
115     public void setParameters(Double  [] p) {
116         for(int k=0;k<p.length;k++) {
117             param[k]=p[k].doubleValue();
118         }
119     }
120 
121     /***********************************
122   * Throws away the parameters of the Filter
123     ************************************/
124     public void removeParameters() {
125         param=null;
126     }
127 
128     /******************************************************
129     * Set the Signal to the specified length scraping or
130   * padding the beginning if necessary
131     *******************************************************/
132     public void setLengthFromEnd(int longueur) {
133     double[] newvec=ArrayMath.setLengthFromEnd(this.evaluate(0),longueur);
134     this.setValues(newvec);
135     }
136 
137   /*****************************************************
138   * Resample the signal using linear interpolation
139   ******************************************************/
140     public void resample(int newl) {
141     double[] newvec=EngineerMath.resample(this.evaluate(0),newl);
142     this.setValues(newvec);
143     }
144 
145     /******************************************************
146     * Set the Signal to the specified length scraping or
147   * padding the end if necessary
148   *******************************************************/
149     public void setLengthFromBeginning(int longueur) {
150     double[] newvec=ArrayMath.setLengthFromBeginning(this.evaluate(0),longueur);
151     this.setValues(newvec);
152     }
153 
154     /***********************************
155   * Set the data for the signal
156     ************************************/
157     public void setData(double[] v) {
158         setValues(v);
159     }
160 
161     /********************************
162     * Fast Wavelet Transform
163     *********************************/
164     public FWTCoef fwt(int J) {
165         if (J>20) {
166             throw new IllegalArgumentException  ("Too many iterations.");
167         }
168         if (J<0) {
169             throw new IllegalArgumentException  ("Cannot have a negative number of iterations.");
170         }
171         double[] data=this.interpolate(0);
172         double[][] fwt=new double[J+1][];
173         for(int j=1;j<=J;j++) {
174             fwt[j]=highpassProject(data);
175             data=lowpassProject(data);
176         }
177         fwt[0]=data;
178         FWTCoef t=new FWTCoef(fwt);
179         return(t);
180     }
181 
182     /********************************
183     * The Fast Wavelet Transform
184   * with Wavelet packets
185   * @param J number of iterations
186   * @param cout cost function
187     *********************************/
188     public FWTPacketCoef fwtPacket(int J, MappingND cout) {
189         if (J>20) {
190             throw new IllegalArgumentException  ("Too many iterations.");
191         }
192         if (J<0) {
193             throw new IllegalArgumentException  ("Cannot have a negative number of iterations.");
194         }
195         double[] data=this.interpolate(0);
196         double[][] fwt=new double[J+1][];
197         double[] choix1,choix2;
198         boolean[] choixStandard=new boolean[J];
199         for(int j=0;j<J;j++) {
200             choix1=highpassProject(data);
201             choix2=lowpassProject(data);
202             if(cout.map(choix1)[0]< cout.map(choix2)[0]) {
203                 fwt[j]=choix1;
204                 data=choix2;
205                 choixStandard[j]=true;
206             } else {
207                 data=choix1;
208                 fwt[j]=choix2;
209                 choixStandard[j]=false;
210             }
211         }
212         fwt[J]=data;
213         FWTPacketCoef t=new FWTPacketCoef(fwt,choixStandard);
214         return(t);
215     }
216 
217 
218     /************************************************
219         * Project the array according to the lowpass Filter
220     * @param v data array
221     * @author Daniel Lemire
222     ************************************************/
223     private double[] lowpassProject (double[] v) {
224         int l=filterdual.previousDimension(v.length);
225         double[] Eche;
226         double[] ans=new double[l];
227       if (param!=null) {
228       for(int i=0;i<l;i++) {
229                 double[] p=new double[param.length];
230                                 System.arraycopy(param,0,p,0,p.length);
231                 Eche=filterdual.lowpass(delta(i,l,normalisation),p);
232               ans[i]=scalarProduct(v,Eche);
233           }
234     } else {
235       for(int i=0;i<l;i++) {
236                 Eche=filterdual.lowpass(delta(i,l,normalisation));
237               ans[i]=scalarProduct(v,Eche);
238           }
239     }
240 
241 /*
242                 for(int i=0;i<l;i++) {
243             if (param!=null) {
244                 double[] p=new double[param.length];
245                                 System.arraycopy(param,0,p,0,p.length);
246                 Eche=filterdual.lowpass(delta(i,l,normalisation),p);
247             } else {
248                 Eche=filterdual.lowpass(delta(i,l,normalisation));
249             }
250             ans[i]=ArrayMath.scalarProduct(v,Eche);
251         }
252 */
253         return(ans);
254     }
255 
256     private static double scalarProduct (double[] w0,double[] w1) {
257         double sortie=0.0;
258         for(int k=0;k<w0.length;k++){
259             sortie+=w0[k]*w1[k];
260         }
261         return(sortie);
262     }
263 
264 
265 
266     /****************************************************
267   * Project the data according to the lowpass Filter
268     * @author Daniel Lemire
269     ****************************************************/
270     public double[] lowpassProject() {
271         double[] data=this.interpolate(0);
272         return(lowpassProject(data));
273     }
274 
275     /******************************************
276   * Project the signal according the the
277   * highpass Filter
278     * @param v data
279     * @author Daniel Lemire
280     ********************************************/
281     private double[] highpassProject (double[] v) {
282         int l=filterdual.previousDimension(v.length);
283         int lOnd=v.length-l;
284         double[] Onde;
285         double[] ans=new double[lOnd];
286         if (param!=null) {
287       for(int i=0;i<lOnd;i++) {
288           double[] p=new double[param.length];
289           System.arraycopy(param,0,p,0,p.length);
290           Onde=filterdual.highpass(delta(i,lOnd,1),p);
291           ans[i]=scalarProduct(v,Onde);
292       }
293     } else {
294 
295       for(int i=0;i<lOnd;i++) {
296 
297 
298                   Onde=filterdual.highpass(delta(i,lOnd,1));
299 
300           ans[i]=scalarProduct(v,Onde);
301       }
302 
303         }
304 
305 /*
306         for(int i=0;i<lOnd;i++) {
307             if (param!=null) {
308                 double[] p=new double[param.length];
309                                 System.arraycopy(param,0,p,0,p.length);
310                 Onde=filterdual.highpass(delta(i,lOnd,1),p);
311             } else {
312                 Onde=filterdual.highpass(delta(i,lOnd,1));
313             }
314             ans[i]=ArrayMath.scalarProduct(v,Onde);
315         }
316 */
317         return(ans);
318     }
319 
320 
321     /********************************************
322   * Project the signal according the the
323   * highpass Filter
324     * @author Daniel Lemire
325     *********************************************/
326     public double[] highpassProject() {
327         double[] data=this.interpolate(0);
328         return(highpassProject(data));
329     }
330 
331     /*************************************************
332   * return a kronecker
333     * @param l length
334     * @param i position
335     * @author Daniel Lemire
336     *************************************************/
337     private double[] delta(int i, int l, double a) {
338         if((i<0)||(i>l)||(l<0)) {
339             throw new IllegalArgumentException  ("This Kronecker doesn't exist.");
340         }
341         double[] v=new double[l];
342         v[i]=a;
343         return(v);
344     }
345 
346     /***************************
347   * Compute the L2 norm of the
348   * signal
349     ****************************/
350     public double norm() {
351         double[] data=this.interpolate(0);
352         return(ArrayMath.norm(data));
353     }
354 
355     /***********************************
356     * @author Don Cross
357         * @author Daniel Lemire
358     ************************************/
359     public Complex[] fft() {
360         double[] data=this.interpolate(0);
361         return(fft(data));
362     }
363 
364         /**
365         * Performs the Fourier transform.
366         * Convenience method for {@link JSci.maths.FourierMath#transform(double[]) FourierMath.transform}.
367         */
368     public static Complex[] fft(double[] data) {
369                 return FourierMath.transform(data);
370         }
371         /**
372         * Performs the Fourier transform.
373         * Convenience method for {@link JSci.maths.FourierMath#transform(Complex[]) FourierMath.transform}.
374         */
375     public static Complex[] fft(Complex[] data) {
376                 return FourierMath.transform(data);
377         }
378   /*********************************
379   * Return the absolute value of
380   * the FFT
381   **********************************/
382     public double[] absFFT() {
383         Complex[] fft=fft();
384         double[] answer=new double[fft.length];
385         for(int i=0;i<fft.length;i++) {
386             answer[i]=fft[i].mod();
387         }
388         return (answer);
389     }
390     public static double[] absFFT(double[] data) {
391         Complex[] fft=fft(data);
392         double[] answer=new double[fft.length];
393         for(int i=0;i<fft.length;i++) {
394             answer[i]=fft[i].mod();
395         }
396         return (answer);
397     }
398         /**************************************
399         * Also noted iFFT in other packages.
400         * This is the inverse to the FFT.
401         ***************************************/
402         public static Complex[] fftInverse(Complex data[]) {
403                 return FourierMath.inverseTransform(data);
404         }
405 
406         /*****************************************
407         * Check if another object is equal to this
408         * Signal object
409         ******************************************/
410     public boolean equals(Signal b) {
411                 return(Arrays.equals(this.getValues(),b.getValues()));
412         }
413     /*******************************************************
414         * Will make the signal a given dimension
415         ********************************************************/
416     public void setDimensionFromEnd(int dimension) {
417         double[] data=this.interpolate(0);
418         double[] ans=new double[dimension];
419         int debut;
420         if(dimension-data.length<0) debut=data.length-dimension;
421         else debut=0;
422         for(int k=debut;k<data.length;k++) {
423             ans[k+dimension-data.length]=data[k];
424         }
425         super.setValues(ans);
426     }
427 
428   /********************************************************
429   * Will make the signal a given dimension
430   *********************************************************/
431     public void setDimensionFromBeginning(int dimension) {
432         double[] data=this.interpolate(0);
433         double[] ans=new double[dimension];
434         int debut;
435         if(dimension-data.length<0) debut=data.length-dimension;
436         else debut=0;
437         for(int k=0;k<data.length-debut;k++) {
438             ans[k]=data[k];
439         }
440         super.setValues(ans);
441     }
442 
443     /********************************************
444   * Simplistic FFT denoising.
445     * @param k frequency to denoised
446     *********************************************/
447     public void denoiseByFFT(int k){
448         if(k<1) {
449             throw new IllegalArgumentException   ("This parameter must be 1 or more : "+k);
450         }
451         double[] data=interpolate(0);
452         if(k>data.length-2) {
453             if (data.length<4) {
454                 throw new IllegalArgumentException   ("Your signal is too short to be denoised : "+data.length+" < 4");
455             }
456             throw new IllegalArgumentException   ("Since you signal has dimension "+data.length+", the parameter must be at most : "+(data.length-2));
457         }
458         Complex[] ff=Signal.fft(data);
459         ff[k+1]=Complex.ZERO;
460         ff[data.length-1-k]=Complex.ZERO;
461         Complex[] tf=Signal.fftInverse(ff);
462         for(int l=0;l<data.length;l++) {
463                 data[l]=tf[l].real();
464                 if(Math.abs(tf[l].imag())>GlobalSettings.ZERO_TOL) {
465                     throw new IllegalArgumentException  ("Complex values detected during synthesis. Please get in touch with Daniel Lemire at Daniel.Lemire@Tintin.net to report this error.");
466                 }
467         }
468     super.setValues(data);
469     }
470 
471   /**********************************************
472   * Return the entropy of the signal
473   ***********************************************/
474   public double entropy () {
475     return(EngineerMath.entropy(this.evaluate(0)));
476   }
477 
478     /*****************************************
479   * Apply the given array as a convolution
480   * Filter and return a new Signal.
481   * As one often want to compare the result
482   * to the original signal, this method is
483   * "safe", that is, it won't change the current
484   * object.
485   * @param f an array containing the coefficients
486   *   of the convolution Filter
487     ******************************************/
488     public Signal filter (double[] f) {
489         double[] data=interpolate(0);
490         if(data.length-(f.length-1)<=0) {
491             throw new IllegalArgumentException  ("Your signal is too short for this Filter : "+data.length+", "+f.length);
492         }
493         double[] ans=new double[data.length-(f.length-1)];
494         for(int k=0;k<data.length-(f.length-1);k++) {
495             for(int l=0;l<f.length;l++) {
496                 ans[k]+=f[l]*data[k+l];
497             }
498         }
499     return(copy (ans));
500     }
501   /******************************************
502   * Apply the median Filter of a window of
503   * size 2*n+1.
504   * exception IllegalArgumentException if the
505   *   parameter n is negative
506   ********************************************/
507   public Signal medianFilter (int n) {
508     if(n<0) throw new IllegalArgumentException  ("The parameter must be positive: "+n+" < 0");
509         double[] data=super.interpolate(0);
510         if(data.length-2*n<=0) {
511             throw new IllegalArgumentException  ("Your signal is too short for this Filter : "+data.length+" - "+(2*n)+" = "+(data.length-2*n));
512         }
513         double[] ans=new double[data.length-2*n];
514     double[] vtemp=new double[2*n+1];
515         for(int k=0;k<data.length-2*n;k++) {
516       for(int l=0;l<2*n+1;l++) {
517         vtemp[l]=data[k+l];
518       }
519       ans[k]=ArrayMath.median(vtemp);
520         }
521     return(copy (ans));
522 
523   }
524   /****************************************************
525   * This denoising method will identify
526   * "short peaks" in the signal and take them away.
527   * Short peaks are defined from a comparison
528   * with the median filtered signal.
529   * Only "significative" peaks are detected (see parameter
530   * p).
531   * This method won't denoise near the boundaries.
532   * "Short" refers here to the time-domain and
533   * not the amplitude.
534   * param p percentage of the range (max-min) considered
535   *   as a significative step
536   * param n length of the peak in the time domain
537   * exception IllegalArgumentException if p is not between
538   *   0 and 1
539   * exception IllegalArgumentException if the
540   *   parameter n is negative
541   *******************************************************/
542   public Signal denoiseShortPeaks (double p, int n) {
543     if((p<0)||(p>1)) {
544       throw new IllegalArgumentException  ("The parameter p must be between 0 and 1: "+p);
545     }
546     double[] values=this.interpolate(0);
547     double range=ArrayMath.max(values)-ArrayMath.min(values);
548     double threshold=range*p;
549     double[] med=(this.medianFilter(n)).interpolate(0);
550     for(int k=n;k<values.length-n;k++) {
551       if(Math.abs(values[k]-med[k-n])>threshold) {
552         values[k]=med[k-n];
553       }
554     }
555     return(copy(values));
556   }
557 }
558 
559
A to Z: JavaDoc & Examples Daily Java News & Articles Open Source Projects Open Source Codes Free Computer Books Remove Frame
Popular Tags