source: src/main/java/agents/anac/y2015/agentBuyogV2/flanagan/math/FourierTransform.java

Last change on this file was 127, checked in by Wouter Pasman, 6 years ago

#41 ROLL BACK of rev.126 . So this version is equal to rev. 125

File size: 96.3 KB
Line 
1/*
2* Fourier Transform
3*
4* This class contains the method for performing a
5* Fast Fourier Transform (FFT) and associated methods
6* e.g. for estimation of a power spectrum, for windowing data,
7* obtaining a time-frequency representation.
8* Basic FFT method is adapted from the Numerical Recipes
9* methods written in the C language:
10* Numerical Recipes in C, The Art of Scientific Computing,
11* W.H. Press, S.A. Teukolsky, W.T. Vetterling & B.P. Flannery,
12* Cambridge University Press, 2nd Edition (1992) pp 496 - 558.
13* (http://www.nr.com/).
14*
15* AUTHOR: Dr Michael Thomas Flanagan
16* DATE: 20 December 2003
17* UPDATES: 26 July 2004, 31 August 2004, 15 June 2005, 27 January 2006
18* UPDATES: 18 February 2006 method correlation correction (thanks to Daniel Mader, Universtät Freiburg -- IMTEK)
19 7 July 2008
20*
21* DOCUMENTATION:
22* See Michael Thomas Flanagan's Java library on-line web page:
23* http://www.ee.ucl.ac.uk/~mflanaga/java/FourierTranasform.html
24* http://www.ee.ucl.ac.uk/~mflanaga/java/
25*
26*
27* Copyright (c) 2003 - 2008 Michael Thomas Flanagan
28*
29* PERMISSION TO COPY:
30* Permission to use, copy and modify this software and its documentation for
31* NON-COMMERCIAL purposes is granted, without fee, provided that an acknowledgement
32* to the author, Michael Thomas Flanagan at www.ee.ucl.ac.uk/~mflanaga, appears in all copies.
33*
34* Dr Michael Thomas Flanagan makes no representations about the suitability
35* or fitness of the software for any or for a particular purpose.
36* Michael Thomas Flanagan shall not be liable for any damages suffered
37* as a result of using, modifying or distributing this software or its derivatives.
38*
39***************************************************************************************/
40
41package agents.anac.y2015.agentBuyogV2.flanagan.math;
42
43import java.awt.Canvas;
44import java.awt.Color;
45import java.awt.Graphics;
46import java.io.Serializable;
47
48import javax.swing.JFrame;
49
50import agents.anac.y2015.agentBuyogV2.flanagan.complex.Complex;
51import agents.anac.y2015.agentBuyogV2.flanagan.io.FileInput;
52import agents.anac.y2015.agentBuyogV2.flanagan.io.FileOutput;
53import agents.anac.y2015.agentBuyogV2.flanagan.plot.PlotGraph;
54
55public class FourierTransform extends Canvas implements Serializable {
56
57 private static final long serialVersionUID = 1L; // serial version unique
58 // identifier
59
60 private Complex[] complexData = null; // array to hold the input data as a
61 // set of Complex numbers
62 private Complex[] complexCorr = null; // corresponding array to hold the
63 // data to be correlated with first
64 // data set
65 private boolean complexDataSet = false; // if true - the complex data input
66 // array has been filled, if false -
67 // it has not.
68 private int originalDataLength = 0; // original data length value; the
69 // working data length may be altered by
70 // deletion or padding
71 private int fftDataLength = 0; // working data length - usually the smallest
72 // power of two that is either equal to
73 // originalDataLength or larger than
74 // originalDataLength
75 private boolean dataAltered = false; // set to true if originalDataLength
76 // altered, e.g. by point deletion
77 // or padding.
78
79 private double[] fftData = null; // array to hold a data set of complex
80 // numbers arranged as alternating
81 // real and imaginary parts, e.g. real_0
82 // imag_0, real_1 imag_1, for the fast
83 // Fourier Transform method
84 private double[] fftCorr = null; // corresponding array to hold the data to
85 // be correlated with first data set
86 private double[] fftResp = null; // corresponding array to hold the response
87 // to be convolved with first data set
88 private boolean fftDataSet = false; // if true - the fftData array has been
89 // filled, if false - it has not.
90
91 private double[] fftDataWindow = null; // array holding fftData array
92 // elements multiplied by the
93 // windowing weights
94 private double[] fftCorrWindow = null; // corresponding array to hold the
95 // data to be correlated with first
96 // data set
97
98 private int windowOption = 0; // Window Option
99 // = 0; no windowing applied (default) -
100 // equivalent to option = 1
101 // = 1; Rectangular (square, box-car)
102 // = 2; Bartlett (triangular)
103 // = 3; Welch
104 // = 4; Hann (Hanning)
105 // = 5; Hamming
106 // = 6; Kaiser
107 // = 7; Gaussian
108 // all window names
109 private String[] windowNames = { "no windowing applied", "Rectangular (square, box-car)", "Bartlett (triangular)",
110 "Welch", "Hann (Hanning)", "Hamming", "Kaiser", "Gaussian" };
111 private String windowName = windowNames[0]; // current window name
112 private double kaiserAlpha = 2.0D; // Kaiser window constant, alpha
113 private double gaussianAlpha = 2.5D; // Gaussian window constant, alpha
114 private double[] weights = null; // windowing weights
115 private boolean windowSet = false; // = true when a windowing option has
116 // been chosen, otherwise = false
117 private boolean windowApplied = false; // = true when data has been
118 // multiplied by windowing weights,
119 // otherwise = false
120 private double sumOfSquaredWeights = 0.0D; // Sum of the windowing weights
121
122 private Complex[] transformedDataComplex = null; // transformed data set of
123 // Complex numbers
124 private double[] transformedDataFft = null; // transformed data set of
125 // double adjacent real and
126 // imaginary parts
127 private boolean fftDone = false; // = false - basicFft has not been called
128 // = true - basicFft has been called
129
130 private double[][] powerSpectrumEstimate = null; // first row - array to
131 // hold frequencies
132 // second row - array to
133 // hold estimated power
134 // density (psd)
135 // spectrum
136 private boolean powSpecDone = false; // = false - PowerSpectrum has not been
137 // called
138 // = true - PowerSpectrum has been
139 // called
140 private int psdNumberOfPoints = 0; // Number of points in the estimated
141 // power spectrum
142
143 private int segmentNumber = 1; // Number of segments into which the data has
144 // been split
145 private int segmentLength = 0; // Number of of data points in a segment
146 private boolean overlap = false; // Data segment overlap option
147 // = true; overlap by half segment
148 // length - smallest spectral variance
149 // per data point
150 // good where data already recorded and
151 // data reduction is after the process
152 // = false; no overlap - smallest
153 // spectral variance per conputer
154 // operation
155 // good for real time data collection
156 // where data reduction is computer
157 // limited
158 private boolean segNumSet = false; // true if segment number has been set
159 private boolean segLenSet = false; // true of segment length has been set
160
161 private double deltaT = 1.0D; // Sampling period (needed only for true
162 // graphical output)
163 private boolean deltaTset = false; // true if sampling period has been set
164
165 private double[][] correlationArray = null; // first row - array to hold
166 // time lags
167 // second row - correlation
168 // between fftDataWindow and
169 // fftCorrWindow
170 private boolean correlateDone = false; // = false - correlation has not been
171 // called
172 // = true - correlation has been
173 // called
174
175 private int numberOfWarnings = 9; // Number of warnings
176 private boolean[] warning = new boolean[numberOfWarnings]; // warnings - if
177 // warning[x] =
178 // true
179 // warningText[x]
180 // is printed
181
182 private int plotLineOption = 0; // PlotPowerSpectrum line option
183 // = 0 points linked by straight line
184 // [default option]
185 // = 1 cubic spline interpolation
186 // = 2 no line - only points
187
188 private int plotPointOption = 0; // PlotPowerSpectrum point option
189 // = 0 no point symbols [default option]
190 // = 1 filled circles
191
192 private double[][] timeFrequency = null; // matrix of time against frequency
193 // mean square powers from
194 // shoert time FT
195 // first row = blank cell
196 // followed by time vector
197 // first column = blank cell
198 // followed by frequency vector
199 // each cell is then the mean
200 // square amplitude at that
201 // frequency and time
202 private boolean shortTimeDone = false; // = true when short time Fourier
203 // Transform has been performed
204 private int numShortFreq = 0; // number of frequency points in short time
205 // Fourier transform
206 private int numShortTimes = 0; // number of time points in short time
207 // Fourier transform
208 private String shortTitle = " "; // Short Time Fourier Transform graph title
209
210 // constructors
211 // No initialisation of the data variables
212 public FourierTransform() {
213 for (int i = 0; i < numberOfWarnings; i++)
214 warning[i] = false;
215
216 }
217
218 // constuctor entering a data array of real numbers
219 public FourierTransform(double[] realData) {
220 this.originalDataLength = realData.length;
221 this.fftDataLength = FourierTransform.nextPowerOfTwo(this.originalDataLength);
222 this.complexData = Complex.oneDarray(this.fftDataLength);
223 for (int i = 0; i < this.originalDataLength; i++) {
224 this.complexData[i].setReal(realData[i]);
225 this.complexData[i].setImag(0.0D);
226 }
227 for (int i = this.originalDataLength; i < this.fftDataLength; i++)
228 this.complexData[i].reset(0.0D, 0.0D);
229 this.complexDataSet = true;
230
231 this.fftData = new double[2 * this.fftDataLength];
232 int j = 0;
233 for (int i = 0; i < this.fftDataLength; i++) {
234 this.fftData[j] = complexData[i].getReal();
235 j++;
236 this.fftData[j] = 0.0D;
237 j++;
238 }
239 this.fftDataSet = true;
240
241 this.fftDataWindow = new double[2 * this.fftDataLength];
242 this.weights = new double[this.fftDataLength];
243 this.sumOfSquaredWeights = windowData(this.fftData, this.fftDataWindow, this.weights);
244
245 this.transformedDataFft = new double[2 * this.fftDataLength];
246 this.transformedDataComplex = Complex.oneDarray(this.fftDataLength);
247 this.segmentLength = this.fftDataLength;
248
249 for (int i = 0; i < numberOfWarnings; i++)
250 warning[i] = false;
251 }
252
253 // constuctor entering a data array of complex numbers
254 public FourierTransform(Complex[] data) {
255 this.originalDataLength = data.length;
256 this.fftDataLength = FourierTransform.nextPowerOfTwo(this.originalDataLength);
257 this.complexData = Complex.oneDarray(this.fftDataLength);
258 for (int i = 0; i < this.originalDataLength; i++) {
259 this.complexData[i] = data[i].copy();
260 }
261 for (int i = this.originalDataLength; i < this.fftDataLength; i++)
262 this.complexData[i].reset(0.0D, 0.0D);
263 this.complexDataSet = true;
264
265 this.fftData = new double[2 * this.fftDataLength];
266 int j = 0;
267 for (int i = 0; i < this.fftDataLength; i++) {
268 this.fftData[j] = complexData[i].getReal();
269 j++;
270 this.fftData[j] = complexData[i].getImag();
271 j++;
272 }
273 this.fftDataSet = true;
274
275 this.fftDataWindow = new double[2 * this.fftDataLength];
276 this.weights = new double[this.fftDataLength];
277 this.sumOfSquaredWeights = windowData(this.fftData, this.fftDataWindow, this.weights);
278
279 this.transformedDataFft = new double[2 * this.fftDataLength];
280 this.transformedDataComplex = Complex.oneDarray(this.fftDataLength);
281 this.segmentLength = this.fftDataLength;
282
283 for (int i = 0; i < numberOfWarnings; i++)
284 warning[i] = false;
285 }
286
287 // Enter a data array of real numbers
288 public void setData(double[] realData) {
289 this.originalDataLength = realData.length;
290 this.fftDataLength = FourierTransform.nextPowerOfTwo(this.originalDataLength);
291 this.complexData = Complex.oneDarray(this.fftDataLength);
292 for (int i = 0; i < this.originalDataLength; i++) {
293 this.complexData[i].setReal(realData[i]);
294 this.complexData[i].setImag(0.0D);
295 }
296 for (int i = this.originalDataLength; i < this.fftDataLength; i++)
297 this.complexData[i].reset(0.0D, 0.0D);
298 this.complexDataSet = true;
299
300 this.fftData = new double[2 * this.fftDataLength];
301 int j = 0;
302 for (int i = 0; i < this.fftDataLength; i++) {
303 this.fftData[j] = complexData[i].getReal();
304 j++;
305 this.fftData[j] = 0.0D;
306 j++;
307 }
308 this.fftDataSet = true;
309
310 this.fftDataWindow = new double[2 * this.fftDataLength];
311 this.weights = new double[this.fftDataLength];
312 this.sumOfSquaredWeights = windowData(this.fftData, this.fftDataWindow, this.weights);
313
314 this.transformedDataFft = new double[2 * this.fftDataLength];
315 this.transformedDataComplex = Complex.oneDarray(this.fftDataLength);
316
317 if (this.segNumSet) {
318 this.setSegmentNumber(this.segmentNumber);
319 } else {
320 if (this.segLenSet) {
321 this.setSegmentLength(this.segmentLength);
322 } else {
323 this.segmentLength = this.fftDataLength;
324 }
325 }
326 }
327
328 // Enter a data array of complex numbers
329 public void setData(Complex[] data) {
330 this.originalDataLength = data.length;
331 this.fftDataLength = FourierTransform.nextPowerOfTwo(this.originalDataLength);
332 this.complexData = Complex.oneDarray(this.fftDataLength);
333 for (int i = 0; i < this.originalDataLength; i++) {
334 this.complexData[i] = data[i].copy();
335 }
336 for (int i = this.originalDataLength; i < this.fftDataLength; i++)
337 this.complexData[i].reset(0.0D, 0.0D);
338 this.complexDataSet = true;
339
340 this.fftData = new double[2 * this.fftDataLength];
341 int j = 0;
342 for (int i = 0; i < this.fftDataLength; i++) {
343 this.fftData[j] = complexData[i].getReal();
344 j++;
345 this.fftData[j] = complexData[i].getImag();
346 j++;
347 }
348 this.fftDataSet = true;
349
350 this.fftDataWindow = new double[2 * this.fftDataLength];
351 this.weights = new double[this.fftDataLength];
352 this.sumOfSquaredWeights = windowData(this.fftData, this.fftDataWindow, this.weights);
353
354 this.transformedDataFft = new double[2 * this.fftDataLength];
355 this.transformedDataComplex = Complex.oneDarray(this.fftDataLength);
356
357 if (this.segNumSet) {
358 this.setSegmentNumber(this.segmentNumber);
359 } else {
360 if (this.segLenSet) {
361 this.setSegmentLength(this.segmentLength);
362 } else {
363 this.segmentLength = this.fftDataLength;
364 }
365 }
366 }
367
368 // Enter a data array of adjacent alternating real and imaginary parts for
369 // fft method, fastFourierTransform
370 public void setFftData(double[] fftdata) {
371 if (fftdata.length % 2 != 0)
372 throw new IllegalArgumentException("data length must be an even number");
373
374 this.originalDataLength = fftdata.length / 2;
375 this.fftDataLength = FourierTransform.nextPowerOfTwo(this.originalDataLength);
376 this.fftData = new double[2 * this.fftDataLength];
377 for (int i = 0; i < 2 * this.originalDataLength; i++)
378 this.fftData[i] = fftdata[i];
379 for (int i = 2 * this.originalDataLength; i < 2 * this.fftDataLength; i++)
380 this.fftData[i] = 0.0D;
381 this.fftDataSet = true;
382
383 this.complexData = Complex.oneDarray(this.fftDataLength);
384 int j = -1;
385 for (int i = 0; i < this.fftDataLength; i++) {
386 this.complexData[i].setReal(this.fftData[++j]);
387 this.complexData[i].setImag(this.fftData[++j]);
388 }
389 this.complexDataSet = true;
390
391 this.fftDataWindow = new double[2 * this.fftDataLength];
392 this.weights = new double[this.fftDataLength];
393 this.sumOfSquaredWeights = windowData(this.fftData, this.fftDataWindow, this.weights);
394
395 this.transformedDataFft = new double[2 * this.fftDataLength];
396 this.transformedDataComplex = Complex.oneDarray(this.fftDataLength);
397
398 if (this.segNumSet) {
399 this.setSegmentNumber(this.segmentNumber);
400 } else {
401 if (this.segLenSet) {
402 this.setSegmentLength(this.segmentLength);
403 } else {
404 this.segmentLength = this.fftDataLength;
405 }
406 }
407 }
408
409 // Get the input data array as Complex
410 public Complex[] getComplexInputData() {
411 if (!this.complexDataSet) {
412 System.out.println("complex data set not entered or calculated - null returned");
413 }
414 return this.complexData;
415 }
416
417 // Get the input data array as adjacent real and imaginary pairs
418 public double[] getAlternateInputData() {
419 if (!this.fftDataSet) {
420 System.out.println("fft data set not entered or calculted - null returned");
421 }
422 return this.fftData;
423 }
424
425 // Get the windowed input data array as windowed adjacent real and imaginary
426 // pairs
427 public double[] getAlternateWindowedInputData() {
428 if (!this.fftDataSet) {
429 System.out.println("fft data set not entered or calculted - null returned");
430 }
431 if (!this.fftDataSet) {
432 System.out.println("fft data set not entered or calculted - null returned");
433 }
434 if (!this.windowApplied) {
435 System.out.println("fft data set has not been multiplied by windowing weights");
436 }
437 return this.fftDataWindow;
438 }
439
440 // get the original number of data points
441 public int getOriginalDataLength() {
442 return this.originalDataLength;
443 }
444
445 // get the actual number of data points
446 public int getUsedDataLength() {
447 return this.fftDataLength;
448 }
449
450 // Set a samplimg period
451 public void setDeltaT(double deltaT) {
452 this.deltaT = deltaT;
453 this.deltaTset = true;
454 }
455
456 // Get the samplimg period
457 public double getDeltaT() {
458 double ret = 0.0D;
459 if (this.deltaTset) {
460 ret = this.deltaT;
461 } else {
462 System.out.println("detaT has not been set - zero returned");
463 }
464 return ret;
465 }
466
467 // Set a Rectangular window option
468 public void setRectangular() {
469 this.windowOption = 1;
470 this.windowSet = true;
471 if (fftDataSet) {
472 this.sumOfSquaredWeights = this.windowData(this.fftData, this.fftDataWindow, this.weights);
473 this.windowApplied = true;
474 }
475 }
476
477 // Set a Bartlett window option
478 public void setBartlett() {
479 this.windowOption = 2;
480 this.windowSet = true;
481 if (fftDataSet) {
482 this.sumOfSquaredWeights = this.windowData(this.fftData, this.fftDataWindow, this.weights);
483 this.windowApplied = true;
484 }
485 }
486
487 // Set a Welch window option
488 public void setWelch() {
489 this.windowOption = 3;
490 this.windowSet = true;
491 if (fftDataSet) {
492 this.sumOfSquaredWeights = this.windowData(this.fftData, this.fftDataWindow, this.weights);
493 this.windowApplied = true;
494 }
495 }
496
497 // Set a Hann window option
498 public void setHann() {
499 this.windowOption = 4;
500 this.windowSet = true;
501 if (fftDataSet) {
502 this.sumOfSquaredWeights = this.windowData(this.fftData, this.fftDataWindow, this.weights);
503 this.windowApplied = true;
504 }
505 }
506
507 // Set a Hamming window option
508 public void setHamming() {
509 this.windowOption = 5;
510 this.windowSet = true;
511 if (fftDataSet) {
512 this.sumOfSquaredWeights = this.windowData(this.fftData, this.fftDataWindow, this.weights);
513 this.windowApplied = true;
514 }
515 }
516
517 // Set a Kaiser window option
518 public void setKaiser(double alpha) {
519 this.kaiserAlpha = alpha;
520 this.windowOption = 6;
521 this.windowSet = true;
522 if (fftDataSet) {
523 this.sumOfSquaredWeights = this.windowData(this.fftData, this.fftDataWindow, this.weights);
524 this.windowApplied = true;
525 }
526 }
527
528 // Set a Kaiser window option
529 // default option for alpha
530 public void setKaiser() {
531 this.windowOption = 6;
532 this.windowSet = true;
533 if (fftDataSet) {
534 this.sumOfSquaredWeights = this.windowData(this.fftData, this.fftDataWindow, this.weights);
535 this.windowApplied = true;
536 }
537 }
538
539 // Set a Gaussian window option
540 public void setGaussian(double alpha) {
541 if (alpha < 2.0D) {
542 alpha = 2.0D;
543 System.out.println("setGaussian; alpha must be greater than or equal to 2 - alpha has been reset to 2");
544 }
545 this.gaussianAlpha = alpha;
546 this.windowOption = 7;
547 this.windowSet = true;
548 if (fftDataSet) {
549 this.sumOfSquaredWeights = this.windowData(this.fftData, this.fftDataWindow, this.weights);
550 this.windowApplied = true;
551 }
552 }
553
554 // Set a Gaussian window option
555 // default option for alpha
556 public void setGaussian() {
557 this.windowOption = 7;
558 this.windowSet = true;
559 if (fftDataSet) {
560 this.sumOfSquaredWeights = this.windowData(this.fftData, this.fftDataWindow, this.weights);
561 this.windowApplied = true;
562 }
563 }
564
565 // Remove windowing
566 public void removeWindow() {
567 this.windowOption = 0;
568 this.windowSet = false;
569 if (fftDataSet) {
570 this.sumOfSquaredWeights = this.windowData(this.fftData, this.fftDataWindow, this.weights);
571 this.windowApplied = false;
572 }
573 }
574
575 // Applies a window to the data
576 private double windowData(double[] data, double[] window, double[] weight) {
577 int m = data.length;
578 int n = m / 2 - 1;
579 int j = 0;
580 double sum = 0.0D;
581 switch (this.windowOption) {
582 // 0. No windowing applied or remove windowing
583 case 0:
584 // 1. Rectangular
585 case 1:
586 for (int i = 0; i <= n; i++) {
587 weight[i] = 1.0D;
588 window[j] = data[j++];
589 window[j] = data[j++];
590 }
591 sum = n + 1;
592 break;
593 // 2. Bartlett
594 case 2:
595 for (int i = 0; i <= n; i++) {
596 weight[i] = 1.0D - Math.abs((i - n / 2) / n / 2);
597 sum += weight[i] * weight[i];
598 window[j] = data[j++] * weight[i];
599 window[j] = data[j++] * weight[i];
600 }
601 break;
602 // 3. Welch
603 case 3:
604 for (int i = 0; i <= n; i++) {
605 weight[i] = 1.0D - Fmath.square((i - n / 2) / n / 2);
606 sum += weight[i] * weight[i];
607 window[j] = data[j++] * weight[i];
608 window[j] = data[j++] * weight[i];
609 }
610 break;
611 // 4. Hann
612 case 4:
613 for (int i = 0; i <= n; i++) {
614 weight[i] = (1.0D - Math.cos(2.0D * i * Math.PI / n)) / 2.0D;
615 sum += weight[i] * weight[i];
616 window[j] = data[j++] * weight[i];
617 window[j] = data[j++] * weight[i];
618 }
619 break;
620 // 5. Hamming
621 case 5:
622 for (int i = 0; i <= n; i++) {
623 weight[i] = 0.54D + 0.46D * Math.cos(2.0D * i * Math.PI / n);
624 sum += weight[i] * weight[i];
625 window[j] = data[j++] * weight[i];
626 window[j] = data[j++] * weight[i];
627 }
628 break;
629 // 6. Kaiser
630 case 6:
631 double denom = FourierTransform.modBesselIo(Math.PI * this.kaiserAlpha);
632 double numer = 0.0D;
633 for (int i = 0; i <= n; i++) {
634 numer = FourierTransform
635 .modBesselIo(Math.PI * this.kaiserAlpha * Math.sqrt(1.0D - Fmath.square(2.0D * i / n - 1.0D)));
636 weight[i] = numer / denom;
637 sum += weight[i] * weight[i];
638 window[j] = data[j++] * weight[i];
639 window[j] = data[j++] * weight[i];
640 }
641 break;
642 // 6. Kaiser
643 case 7:
644 for (int i = 0; i <= n; i++) {
645 weight[i] = Math.exp(-0.5D * Fmath.square(this.gaussianAlpha * (2 * i - n) / n));
646 sum += weight[i] * weight[i];
647 window[j] = data[j++] * weight[i];
648 window[j] = data[j++] * weight[i];
649 }
650 break;
651 }
652 return sum;
653 }
654
655 // return modified Bessel Function of the zeroth order (for Kaiser window)
656 // after numerical Recipe's bessi0
657 // - Abramowitz and Stegun coeeficients
658 public static double modBesselIo(double arg) {
659 double absArg = 0.0D;
660 double poly = 0.0D;
661 double bessel = 0.0D;
662
663 if ((absArg = Math.abs(arg)) < 3.75) {
664 poly = arg / 3.75;
665 poly *= poly;
666 bessel = 1.0D + poly * (3.5156229D + poly * (3.08989424D
667 + poly * (1.2067492D + poly * (0.2659732 + poly * (0.360768e-1 + poly * 0.45813e-2)))));
668 } else {
669 bessel = (Math.exp(absArg) / Math.sqrt(absArg)) * (0.39894228D + poly * (0.1328592e-1D
670 + poly * (0.225319e-2 + poly * (-0.157565e-2 + poly * (0.916281e-2 + poly * (-0.2057706e-1
671 + poly * (0.2635537e-1 + poly * (-0.1647633e-1 + poly * 0.392377e-2))))))));
672 }
673 return bessel;
674 }
675
676 // get window option - see above for options
677 public String getWindowOption() {
678 String option = " ";
679 switch (this.windowOption) {
680 case 0:
681 option = "No windowing applied";
682 break;
683 case 1:
684 option = "Rectangular";
685 break;
686 case 2:
687 option = "Bartlett";
688 break;
689 case 3:
690 option = "Welch";
691 break;
692 case 4:
693 option = "Hann";
694 break;
695 case 5:
696 option = "Hamming";
697 break;
698 case 6:
699 option = "Kaiser";
700 break;
701 case 7:
702 option = "Gaussian";
703 break;
704 }
705 return option;
706 }
707
708 // Get the windowing weights
709 public double[] getWeights() {
710 return this.weights;
711 }
712
713 // set the number of segments
714 public void setSegmentNumber(int sNum) {
715 this.segmentNumber = sNum;
716 this.segNumSet = true;
717 if (this.segLenSet)
718 this.segLenSet = false;
719 }
720
721 // set the segment length
722 public void setSegmentLength(int sLen) {
723 this.segmentLength = sLen;
724 this.segLenSet = true;
725 if (this.segNumSet)
726 this.segNumSet = false;
727 }
728
729 // check and set up the segments
730 private void checkSegmentDetails() {
731 if (!this.fftDataSet)
732 throw new IllegalArgumentException("No fft data has been entered or calculated");
733 if (this.fftDataLength < 2)
734 throw new IllegalArgumentException("More than one point, MANY MORE, are needed");
735
736 // check if data number is even
737 if (this.fftDataLength % 2 != 0) {
738 System.out.println("Number of data points must be an even number");
739 System.out.println("last point deleted");
740 this.fftDataLength -= 1;
741 this.dataAltered = true;
742 this.warning[0] = true;
743 }
744
745 // check segmentation with no overlap
746 if (this.segNumSet && !this.overlap) {
747 if (this.fftDataLength % this.segmentNumber == 0) {
748 int segL = this.fftDataLength / this.segmentNumber;
749 if (FourierTransform.checkPowerOfTwo(segL)) {
750 this.segmentLength = segL;
751 this.segLenSet = true;
752 } else {
753 System.out.println("segment length is not an integer power of two");
754 System.out.println("segment length reset to total data length, i.e. no segmentation");
755 warning[1] = true;
756 this.segmentNumber = 1;
757 this.segmentLength = this.fftDataLength;
758 this.segLenSet = true;
759 }
760 } else {
761 System.out.println("total data length divided by the number of segments is not an integer");
762 System.out.println("segment length reset to total data length, i.e. no segmentation");
763 warning[2] = true;
764 this.segmentNumber = 1;
765 this.segmentLength = this.fftDataLength;
766 this.segLenSet = true;
767 }
768 }
769
770 if (this.segLenSet && !this.overlap) {
771 if (this.fftDataLength % this.segmentLength == 0) {
772 if (FourierTransform.checkPowerOfTwo(this.segmentLength)) {
773 this.segmentNumber = this.fftDataLength / this.segmentLength;
774 this.segNumSet = true;
775 } else {
776 System.out.println("segment length is not an integer power of two");
777 System.out.println("segment length reset to total data length, i.e. no segmentation");
778 warning[1] = true;
779 this.segmentNumber = 1;
780 this.segmentLength = this.fftDataLength;
781 this.segNumSet = true;
782 }
783 } else {
784 System.out.println("total data length divided by the segment length is not an integer");
785 System.out.println("segment length reset to total data length, i.e. no segmentation");
786 warning[3] = true;
787 this.segmentNumber = 1;
788 this.segmentLength = this.fftDataLength;
789 this.segNumSet = true;
790 }
791 }
792
793 // check segmentation with overlap
794 if (this.segNumSet && this.overlap) {
795 if (this.fftDataLength % (this.segmentNumber + 1) == 0) {
796 int segL = 2 * this.fftDataLength / (this.segmentNumber + 1);
797 if (FourierTransform.checkPowerOfTwo(segL)) {
798 this.segmentLength = segL;
799 this.segLenSet = true;
800 } else {
801 System.out.println("segment length is not an integer power of two");
802 System.out.println("segment length reset to total data length, i.e. no segmentation");
803 warning[1] = true;
804 this.segmentNumber = 1;
805 this.segmentLength = this.fftDataLength;
806 this.segLenSet = true;
807 this.overlap = false;
808 }
809 } else {
810 System.out.println("total data length divided by the number of segments plus one is not an integer");
811 System.out.println("segment length reset to total data length, i.e. no segmentation");
812 warning[4] = true;
813 this.segmentNumber = 1;
814 this.segmentLength = this.fftDataLength;
815 this.segLenSet = true;
816 this.overlap = false;
817 }
818 }
819
820 if (this.segLenSet && this.overlap) {
821 if ((2 * this.fftDataLength) % this.segmentLength == 0) {
822 if (FourierTransform.checkPowerOfTwo(this.segmentLength)) {
823 this.segmentNumber = (2 * this.fftDataLength) / this.segmentLength - 1;
824 this.segNumSet = true;
825 } else {
826 System.out.println("segment length is not an integer power of two");
827 System.out.println("segment length reset to total data length, i.e. no segmentation");
828 warning[1] = true;
829 this.segmentNumber = 1;
830 this.segmentLength = this.fftDataLength;
831 this.segNumSet = true;
832 this.overlap = false;
833 }
834 } else {
835 System.out.println("twice the total data length divided by the segment length is not an integer");
836 System.out.println("segment length reset to total data length, i.e. no segmentation");
837 warning[5] = true;
838 this.segmentNumber = 1;
839 this.segmentLength = this.fftDataLength;
840 this.segNumSet = true;
841 this.overlap = false;
842 }
843 }
844
845 if (!this.segNumSet && !this.segLenSet) {
846 this.segmentNumber = 1;
847 this.segNumSet = true;
848 this.overlap = false;
849 }
850
851 if (this.overlap && this.segmentNumber < 2) {
852 System.out.println("Overlap is not possible with less than two segments.");
853 System.out.println("Overlap option has been reset to 'no overlap' i.e. to false.");
854 this.overlap = false;
855 this.segmentNumber = 1;
856 this.segNumSet = true;
857 warning[6] = true;
858 }
859
860 // check no segmentation option
861 int segLno = 0;
862 int segNno = 0;
863 int segLov = 0;
864 int segNov = 0;
865
866 if (this.segmentNumber == 1) {
867 // check if data number is a power of two
868 if (!FourierTransform.checkPowerOfTwo(this.fftDataLength)) {
869 boolean test0 = true;
870 boolean test1 = true;
871 boolean test2 = true;
872 int newL = 0;
873 int ii = 2;
874 // not a power of two - check segmentation options
875 // no overlap option
876 while (test0) {
877 newL = this.fftDataLength / ii;
878 if (FourierTransform.checkPowerOfTwo(newL) && (this.fftDataLength % ii) == 0) {
879 test0 = false;
880 segLno = newL;
881 segNno = ii;
882 } else {
883 if (newL < 2) {
884 test1 = false;
885 test0 = false;
886 } else {
887 ii++;
888 }
889 }
890 }
891 test0 = true;
892 ii = 2;
893 // overlap option
894 while (test0) {
895 newL = 2 * (this.fftDataLength / (ii + 1));
896 if (FourierTransform.checkPowerOfTwo(newL) && (this.fftDataLength % (ii + 1)) == 0) {
897 test0 = false;
898 segLov = newL;
899 segNov = ii;
900 } else {
901 if (newL < 2) {
902 test2 = false;
903 test0 = false;
904 } else {
905 ii++;
906 }
907 }
908 }
909 // compare overlap and no overlap options
910 boolean setSegment = true;
911 int segL = 0;
912 int segN = 0;
913 boolean ovrlp = false;
914 if (test1) {
915 if (test2) {
916 if (segLov > segLno) {
917 segL = segLov;
918 segN = segNov;
919 ovrlp = true;
920 } else {
921 segL = segLno;
922 segN = segNno;
923 ovrlp = false;
924 }
925 } else {
926 segL = segLno;
927 segN = segNno;
928 ovrlp = false;
929 }
930 } else {
931 if (test2) {
932 segL = segLov;
933 segN = segNov;
934 ovrlp = true;
935 } else {
936 setSegment = false;
937 }
938 }
939
940 // compare segmentation and zero padding
941 if (setSegment && (this.originalDataLength - segL <= this.fftDataLength - this.originalDataLength)) {
942 System.out.println("Data length is not an integer power of two");
943 System.out.println("Data cannot be transformed as a single segment");
944 System.out.print("The data has been split into " + segN + " segments of length " + segL);
945 if (ovrlp) {
946 System.out.println(" with 50% overlap");
947 } else {
948 System.out.println(" with no overlap");
949 }
950 this.segmentLength = segL;
951 this.segmentNumber = segN;
952 this.overlap = ovrlp;
953 this.warning[7] = true;
954 } else {
955 System.out.println("Data length is not an integer power of two");
956 if (this.dataAltered) {
957 System.out.println(
958 "Deleted point has been restored and the data has been padded with zeros to give a power of two length");
959 this.warning[0] = false;
960 } else {
961 System.out.println("Data has been padded with zeros to give a power of two length");
962 }
963 // this.fftDataLength = this.fftDataLength;
964 this.warning[8] = true;
965 }
966 }
967 }
968 }
969
970 private void printWarnings(FileOutput fout) {
971 if (warning[0]) {
972 fout.println("WARNING!");
973 fout.println("Number of data points must be an even number");
974 fout.println("The last point was deleted");
975 fout.println();
976 }
977
978 if (warning[1]) {
979 fout.println("WARNING!");
980 fout.println("Segment length was not an integer power of two");
981 fout.println("Segment length was reset to total data length, i.e. no segmentation");
982 fout.println();
983 }
984
985 if (warning[2]) {
986 fout.println("WARNING!");
987 fout.println("Total data length divided by the number of segments was not an integer");
988 fout.println("Segment length was reset to total data length, i.e. no segmentation");
989 fout.println();
990 }
991
992 if (warning[3]) {
993 fout.println("WARNING!");
994 fout.println("Total data length divided by the segment length was not an integer");
995 fout.println("Segment length was reset to total data length, i.e. no segmentation");
996 fout.println();
997 }
998
999 if (warning[4]) {
1000 fout.println("WARNING!");
1001 fout.println("Total data length divided by the number of segments plus one was not an integer");
1002 fout.println("Segment length was reset to total data length, i.e. no segmentation");
1003 fout.println();
1004 }
1005
1006 if (warning[5]) {
1007 fout.println("WARNING!");
1008 fout.println("Twice the total data length divided by the segment length was not an integer");
1009 fout.println("Segment length was reset to total data length, i.e. no segmentation");
1010 fout.println();
1011 }
1012
1013 if (warning[6]) {
1014 fout.println("WARNING!");
1015 fout.println("Overlap is not possible with less than two segments");
1016 fout.println("Overlap option has been reset to 'no overlap' i.e. to false");
1017 fout.println();
1018 }
1019
1020 if (warning[7]) {
1021 fout.println("WARNING!");
1022 fout.println("Data length was not an integer power of two");
1023 fout.println("The data could not be transformed as a single segment");
1024 fout.print("The data has been split into " + this.segmentNumber + " segment/s of length "
1025 + this.segmentLength);
1026 if (this.overlap) {
1027 fout.println(" with 50% overlap");
1028 } else {
1029 fout.println(" with no overlap");
1030 }
1031 fout.println();
1032 }
1033
1034 if (warning[8]) {
1035 fout.println("WARNING!");
1036 fout.println("Data length was not an integer power of two");
1037 fout.println("Data has been padded with " + (this.fftDataLength - this.originalDataLength)
1038 + " zeros to give an integer power of two length");
1039 fout.println();
1040 }
1041 }
1042
1043 // get the number of segments
1044 public int getSegmentNumber() {
1045 return this.segmentNumber;
1046 }
1047
1048 // get the segment length
1049 public int getSegmentLength() {
1050 return this.segmentLength;
1051 }
1052
1053 // set overlap option - see above (head of program comment lines) for option
1054 // description
1055 public void setOverlapOption(boolean overlapOpt) {
1056 boolean old = this.overlap;
1057 this.overlap = overlapOpt;
1058 if (old != this.overlap) {
1059 if (this.fftDataSet) {
1060 this.setSegmentNumber(this.segmentNumber);
1061 }
1062 }
1063 }
1064
1065 // get overlap option - see above for options
1066 public boolean getOverlapOption() {
1067 return this.overlap;
1068 }
1069
1070 // calculate the number of data points given the:
1071 // segment length (segLen), number of segments (segNum)
1072 // and the overlap option (overlap: true - overlap, false - no overlap)
1073 public static int calcDataLength(boolean overlap, int segLen, int segNum) {
1074 if (overlap) {
1075 return (segNum + 1) * segLen / 2;
1076 } else {
1077 return segNum * segLen;
1078 }
1079 }
1080
1081 // Method for performing a Fast Fourier Transform
1082 public void transform() {
1083
1084 // set up data array
1085 int isign = 1;
1086 if (!this.fftDataSet)
1087 throw new IllegalArgumentException("No data has been entered for the Fast Fourier Transform");
1088 if (this.originalDataLength != this.fftDataLength) {
1089 System.out.println("Fast Fourier Transform data length ," + this.originalDataLength
1090 + ", is not an integer power of two");
1091 System.out.println(
1092 "WARNING!!! Data has been padded with zeros to fill to nearest integer power of two length "
1093 + this.fftDataLength);
1094 }
1095
1096 // Perform fft
1097 double[] hold = new double[this.fftDataLength * 2];
1098 for (int i = 0; i < this.fftDataLength * 2; i++)
1099 hold[i] = this.fftDataWindow[i];
1100 basicFft(hold, this.fftDataLength, isign);
1101 for (int i = 0; i < this.fftDataLength * 2; i++)
1102 this.transformedDataFft[i] = hold[i];
1103
1104 // fill transformed data arrays
1105 for (int i = 0; i < this.fftDataLength; i++) {
1106 this.transformedDataComplex[i].reset(this.transformedDataFft[2 * i], this.transformedDataFft[2 * i + 1]);
1107 }
1108 }
1109
1110 // Method for performing an inverse Fast Fourier Transform
1111 public void inverse() {
1112
1113 // set up data array
1114 int isign = -1;
1115 if (!this.fftDataSet)
1116 throw new IllegalArgumentException("No data has been entered for the inverse Fast Fourier Transform");
1117 if (this.originalDataLength != this.fftDataLength) {
1118 System.out.println("Fast Fourier Transform data length ," + this.originalDataLength
1119 + ", is not an integer power of two");
1120 System.out.println(
1121 "WARNING!!! Data has been padded with zeros to fill to nearest integer power of two length "
1122 + this.fftDataLength);
1123 }
1124
1125 // Perform inverse fft
1126 double[] hold = new double[this.fftDataLength * 2];
1127 for (int i = 0; i < this.fftDataLength * 2; i++)
1128 hold[i] = this.fftDataWindow[i];
1129 basicFft(hold, this.fftDataLength, isign);
1130
1131 for (int i = 0; i < this.fftDataLength * 2; i++)
1132 this.transformedDataFft[i] = hold[i] / this.fftDataLength;
1133
1134 // fill transformed data arrays
1135 for (int i = 0; i < this.fftDataLength; i++) {
1136 this.transformedDataComplex[i].reset(this.transformedDataFft[2 * i], this.transformedDataFft[2 * i + 1]);
1137 }
1138 }
1139
1140 // Base method for performing a Fast Fourier Transform
1141 // Based on the Numerical Recipes procedure four1
1142 // If isign is set to +1 this method replaces fftData[0 to 2*nn-1] by its
1143 // discrete Fourier Transform
1144 // If isign is set to -1 this method replaces fftData[0 to 2*nn-1] by nn
1145 // times its inverse discrete Fourier Transform
1146 // nn MUST be an integer power of 2. This is not checked for in this method,
1147 // fastFourierTransform(...), for speed.
1148 // If not checked for by the calling method, e.g. powerSpectrum(...) does,
1149 // the method checkPowerOfTwo() may be used to check this.
1150 // The real and imaginary parts of the data are stored adjacently
1151 // i.e. fftData[0] holds the real part, fftData[1] holds the corresponding
1152 // imaginary part of a data point
1153 // data array and data array length over 2 (nn) transferred as arguments
1154 // result NOT returned to this.transformedDataFft
1155 // Based on the Numerical Recipes procedure four1
1156 public void basicFft(double[] data, int nn, int isign) {
1157 double dtemp = 0.0D, wtemp = 0.0D, tempr = 0.0D, tempi = 0.0D;
1158 double theta = 0.0D, wr = 0.0D, wpr = 0.0D, wpi = 0.0D, wi = 0.0D;
1159 int istep = 0, m = 0, mmax = 0;
1160 int n = nn << 1;
1161 int j = 1;
1162 int jj = 0;
1163 for (int i = 1; i < n; i += 2) {
1164 jj = j - 1;
1165 if (j > i) {
1166 int ii = i - 1;
1167 dtemp = data[jj];
1168 data[jj] = data[ii];
1169 data[ii] = dtemp;
1170 dtemp = data[jj + 1];
1171 data[jj + 1] = data[ii + 1];
1172 data[ii + 1] = dtemp;
1173 }
1174 m = n >> 1;
1175 while (m >= 2 && j > m) {
1176 j -= m;
1177 m >>= 1;
1178 }
1179 j += m;
1180 }
1181 mmax = 2;
1182 while (n > mmax) {
1183 istep = mmax << 1;
1184 theta = isign * (6.28318530717959D / mmax);
1185 wtemp = Math.sin(0.5D * theta);
1186 wpr = -2.0D * wtemp * wtemp;
1187 wpi = Math.sin(theta);
1188 wr = 1.0D;
1189 wi = 0.0D;
1190 for (m = 1; m < mmax; m += 2L) {
1191 for (int i = m; i <= n; i += istep) {
1192 int ii = i - 1;
1193 jj = ii + mmax;
1194 tempr = wr * data[jj] - wi * data[jj + 1];
1195 tempi = wr * data[jj + 1] + wi * data[jj];
1196 data[jj] = data[ii] - tempr;
1197 data[jj + 1] = data[ii + 1] - tempi;
1198 data[ii] += tempr;
1199 data[ii + 1] += tempi;
1200 }
1201 wr = (wtemp = wr) * wpr - wi * wpi + wr;
1202 wi = wi * wpr + wtemp * wpi + wi;
1203 }
1204 mmax = istep;
1205 }
1206 }
1207
1208 // Get the transformed data as Complex
1209 public Complex[] getTransformedDataAsComplex() {
1210 return this.transformedDataComplex;
1211 }
1212
1213 // Get the transformed data array as adjacent real and imaginary pairs
1214 public double[] getTransformedDataAsAlternate() {
1215 return this.transformedDataFft;
1216 }
1217
1218 // Performs and returns results a fft power spectrum density (psd)
1219 // estimation
1220 // of unsegmented, segmented or segemented and overlapped data
1221 // data in array fftDataWindow
1222 public double[][] powerSpectrum() {
1223
1224 this.checkSegmentDetails();
1225
1226 this.psdNumberOfPoints = this.segmentLength / 2;
1227 this.powerSpectrumEstimate = new double[2][this.psdNumberOfPoints];
1228
1229 if (!overlap && this.segmentNumber < 2) {
1230 // Unsegmented and non-overlapped data
1231
1232 // set up data array
1233 int isign = 1;
1234 if (!this.fftDataSet)
1235 throw new IllegalArgumentException("No data has been entered for the Fast Fourier Transform");
1236 if (!FourierTransform.checkPowerOfTwo(this.fftDataLength))
1237 throw new IllegalArgumentException("Fast Fourier Transform data length ," + this.fftDataLength
1238 + ", is not an integer power of two");
1239
1240 // perform fft
1241 double[] hold = new double[this.fftDataLength * 2];
1242 for (int i = 0; i < this.fftDataLength * 2; i++)
1243 hold[i] = this.fftDataWindow[i];
1244 basicFft(hold, this.fftDataLength, isign);
1245 for (int i = 0; i < this.fftDataLength * 2; i++)
1246 this.transformedDataFft[i] = hold[i];
1247
1248 // fill transformed data arrays
1249 for (int i = 0; i < this.fftDataLength; i++) {
1250 this.transformedDataComplex[i].reset(this.transformedDataFft[2 * i],
1251 this.transformedDataFft[2 * i + 1]);
1252 }
1253
1254 // obtain weighted mean square amplitudes
1255 this.powerSpectrumEstimate[1][0] = Fmath.square(hold[0]) + Fmath.square(hold[1]);
1256 for (int i = 1; i < this.psdNumberOfPoints; i++) {
1257 this.powerSpectrumEstimate[1][i] = Fmath.square(hold[2 * i]) + Fmath.square(hold[2 * i + 1])
1258 + Fmath.square(hold[2 * this.segmentLength - 2 * i])
1259 + Fmath.square(hold[2 * this.segmentLength - 2 * i + 1]);
1260 }
1261
1262 // Normalise
1263 for (int i = 0; i < this.psdNumberOfPoints; i++) {
1264 this.powerSpectrumEstimate[1][i] = 2.0D * this.powerSpectrumEstimate[1][i]
1265 / (this.fftDataLength * this.sumOfSquaredWeights);
1266 }
1267
1268 // Calculate frequencies
1269 for (int i = 0; i < this.psdNumberOfPoints; i++) {
1270 this.powerSpectrumEstimate[0][i] = (double) i / ((double) this.segmentLength * this.deltaT);
1271 }
1272 } else {
1273 // Segmented or segmented and overlapped data
1274 this.powerSpectrumEstimate = powerSpectrumSeg();
1275 }
1276
1277 this.powSpecDone = true;
1278
1279 return this.powerSpectrumEstimate;
1280 }
1281
1282 // Performs and returns results a fft power spectrum density (psd)
1283 // estimation
1284 // of unsegmented, segmented or segemented and overlaped data
1285 // data read in from a text file
1286 public double[][] powerSpectrum(String fileName) {
1287
1288 if (!FourierTransform.checkPowerOfTwo(this.segmentLength))
1289 throw new IllegalArgumentException("Fast Fourier Transform segment length ," + this.segmentLength
1290 + ", is not an integer power of two");
1291
1292 FileInput fin = new FileInput(fileName);
1293
1294 this.psdNumberOfPoints = this.segmentLength / 2;
1295 this.powerSpectrumEstimate = new double[2][this.psdNumberOfPoints];
1296 this.fftDataLength = FourierTransform.calcDataLength(this.overlap, this.segmentLength, this.segmentNumber);
1297
1298 if (!overlap && this.segmentNumber < 2) {
1299 // Unsegmented and non-overlapped data
1300
1301 // read in data
1302 this.fftData = new double[2 * this.fftDataLength];
1303 int j = -1;
1304 for (int i = 0; i < this.segmentLength; i++) {
1305 this.fftData[++j] = fin.readDouble();
1306 this.fftData[++j] = fin.readDouble();
1307 }
1308
1309 this.complexData = Complex.oneDarray(this.fftDataLength);
1310 j = -1;
1311 for (int i = 0; i < this.fftDataLength; i++) {
1312 this.complexData[i].setReal(this.fftData[++j]);
1313 this.complexData[i].setImag(this.fftData[++j]);
1314 }
1315
1316 this.fftDataWindow = new double[2 * this.fftDataLength];
1317 this.sumOfSquaredWeights = this.windowData(this.fftData, this.fftDataWindow, this.weights);
1318
1319 // perform fft
1320 int isign = 1;
1321 double[] hold = new double[this.fftDataLength * 2];
1322 for (int i = 0; i < this.fftDataLength * 2; i++)
1323 hold[i] = this.fftDataWindow[i];
1324 basicFft(hold, this.fftDataLength, isign);
1325 for (int i = 0; i < this.fftDataLength * 2; i++)
1326 this.transformedDataFft[i] = hold[i];
1327
1328 // fill transformed data arrays
1329 for (int i = 0; i < this.fftDataLength; i++) {
1330 this.transformedDataComplex[i].reset(this.transformedDataFft[2 * i],
1331 this.transformedDataFft[2 * i + 1]);
1332 }
1333
1334 // obtain weighted mean square amplitudes
1335 this.powerSpectrumEstimate[1][0] = Fmath.square(hold[0]) + Fmath.square(hold[1]);
1336 for (int i = 1; i < this.psdNumberOfPoints; i++) {
1337 this.powerSpectrumEstimate[1][i] = Fmath.square(hold[2 * i]) + Fmath.square(hold[2 * i + 1])
1338 + Fmath.square(hold[2 * this.segmentLength - 2 * i])
1339 + Fmath.square(hold[2 * this.segmentLength - 2 * i + 1]);
1340 }
1341
1342 // Normalise
1343 for (int i = 0; i < this.psdNumberOfPoints; i++) {
1344 this.powerSpectrumEstimate[1][i] = 2.0D * this.powerSpectrumEstimate[1][i]
1345 / (this.fftDataLength * this.sumOfSquaredWeights);
1346 }
1347
1348 // Calculate frequencies
1349 for (int i = 0; i < this.psdNumberOfPoints; i++) {
1350 this.powerSpectrumEstimate[0][i] = (double) i / ((double) this.segmentLength * this.deltaT);
1351 }
1352
1353 } else {
1354 // Segmented or segmented and overlapped data
1355 this.powerSpectrumEstimate = powerSpectrumSeg(fin);
1356 }
1357
1358 this.powSpecDone = true;
1359
1360 return this.powerSpectrumEstimate;
1361 }
1362
1363 // Performs and returns results a fft power spectrum density (psd)
1364 // estimation of segmented or segemented and overlaped data
1365 // Data in fftDataWindow array
1366 // Private method for PowerSpectrum (see above)
1367 private double[][] powerSpectrumSeg() {
1368
1369 // set up segment details
1370 int segmentStartIndex = 0;
1371 int segmentStartIncrement = this.segmentLength;
1372 if (this.overlap)
1373 segmentStartIncrement /= 2;
1374 double[] data = new double[2 * this.segmentLength]; // holds data and
1375 // transformed data
1376 // for working
1377 // segment
1378 this.psdNumberOfPoints = this.segmentLength / 2; // number of PSD points
1379 double[] segPSD = new double[this.psdNumberOfPoints]; // holds psd for
1380 // working
1381 // segment
1382 double[][] avePSD = new double[2][this.psdNumberOfPoints]; // first row
1383 // -
1384 // frequencies
1385 // second
1386 // row -
1387 // accumaltes
1388 // psd for
1389 // averaging
1390 // and then
1391 // the
1392 // averaged
1393 // psd
1394
1395 // initialis psd array and transform option
1396 for (int j = 0; j < this.psdNumberOfPoints; j++)
1397 avePSD[1][j] = 0.0D;
1398 int isign = 1;
1399
1400 // loop through segments
1401 for (int i = 1; i <= this.segmentNumber; i++) {
1402
1403 // collect segment data
1404 for (int j = 0; j < 2 * this.segmentLength; j++)
1405 data[j] = this.fftData[segmentStartIndex + j];
1406
1407 // window data
1408 if (i == 1) {
1409 this.sumOfSquaredWeights = this.windowData(data, data, this.weights);
1410 } else {
1411 int k = 0;
1412 for (int j = 0; j < this.segmentLength; j++) {
1413 data[k] = data[k] * this.weights[j];
1414 data[++k] = data[k] * this.weights[j];
1415 ++k;
1416 }
1417 }
1418
1419 // perform fft on windowed segment
1420 basicFft(data, this.segmentLength, isign);
1421
1422 // obtain weighted mean square amplitudes
1423 segPSD[0] = Fmath.square(data[0]) + Fmath.square(data[1]);
1424 for (int j = 1; j < this.psdNumberOfPoints; j++) {
1425 segPSD[j] = Fmath.square(data[2 * j]) + Fmath.square(data[2 * j + 1])
1426 + Fmath.square(data[2 * this.segmentLength - 2 * j])
1427 + Fmath.square(data[2 * this.segmentLength - 2 * j + 1]);
1428 }
1429
1430 // Normalise
1431 for (int j = 0; j < this.psdNumberOfPoints; j++) {
1432 segPSD[j] = 2.0D * segPSD[j] / (this.segmentLength * this.sumOfSquaredWeights);
1433 }
1434
1435 // accumalate for averaging
1436 for (int j = 0; j < this.psdNumberOfPoints; j++)
1437 avePSD[1][j] += segPSD[j];
1438
1439 // increment segment start index
1440 segmentStartIndex += segmentStartIncrement;
1441 }
1442
1443 // average all segments
1444 for (int j = 0; j < this.psdNumberOfPoints; j++)
1445 avePSD[1][j] /= this.segmentNumber;
1446
1447 // Calculate frequencies
1448 for (int i = 0; i < this.psdNumberOfPoints; i++) {
1449 avePSD[0][i] = (double) i / ((double) this.segmentLength * this.deltaT);
1450 }
1451
1452 return avePSD;
1453 }
1454
1455 // Performs and returns results a fft power spectrum density (psd)
1456 // estimation of segmented or segemented and overlaped data
1457 // Data read in from a text file
1458 // Private method for PowerSpectrum(fileName) (see above)
1459 private double[][] powerSpectrumSeg(FileInput fin) {
1460
1461 // set up segment details
1462 double[] data = new double[2 * this.segmentLength]; // holds data and
1463 // transformed data
1464 // for working
1465 // segment
1466 this.weights = new double[this.segmentLength]; // windowing weights for
1467 // segment
1468 double[] hold = new double[2 * this.segmentLength]; // working array
1469 this.psdNumberOfPoints = this.segmentLength / 2; // number of PSD points
1470 double[] segPSD = new double[this.psdNumberOfPoints]; // holds psd for
1471 // working
1472 // segment
1473 double[][] avePSD = new double[2][this.psdNumberOfPoints]; // first row
1474 // -
1475 // frequencies
1476 // second
1477 // row -
1478 // accumaltes
1479 // psd for
1480 // averaging
1481 // and then
1482 // the
1483 // averaged
1484 // psd
1485
1486 // initialise psd array and fft option
1487 for (int j = 0; j < this.psdNumberOfPoints; j++)
1488 avePSD[1][j] = 0.0D;
1489 int isign = 1;
1490
1491 // calculate window weights
1492 this.sumOfSquaredWeights = this.windowData(hold, hold, this.weights);
1493
1494 if (this.overlap) {
1495 // overlapping segments
1496
1497 // read in first half segment
1498 for (int j = 0; j < this.segmentLength; j++) {
1499 data[j] = fin.readDouble();
1500 }
1501
1502 // loop through segments
1503 for (int i = 1; i <= this.segmentNumber; i++) {
1504
1505 // read in next half segment
1506 for (int j = 0; j < this.segmentLength; j++) {
1507 data[j + this.segmentLength] = fin.readDouble();
1508 }
1509
1510 // window data
1511 int k = -1;
1512 for (int j = 0; j < this.segmentLength; j++) {
1513 data[++k] = data[k] * this.weights[j];
1514 data[++k] = data[k] * this.weights[j];
1515 }
1516
1517 // perform fft on windowed segment
1518 basicFft(data, this.segmentLength, isign);
1519
1520 // obtain weighted mean square amplitudes
1521 segPSD[0] = Fmath.square(data[0]) + Fmath.square(data[1]);
1522 for (int j = 1; j < this.psdNumberOfPoints; j++) {
1523 segPSD[j] = Fmath.square(data[2 * j]) + Fmath.square(data[2 * j + 1])
1524 + Fmath.square(data[2 * this.segmentLength - 2 * j])
1525 + Fmath.square(data[2 * this.segmentLength - 2 * j + 1]);
1526 }
1527
1528 // Normalise
1529 for (int j = 0; j < this.psdNumberOfPoints; j++) {
1530 segPSD[j] = 2.0D * segPSD[j] / (this.segmentLength * this.sumOfSquaredWeights);
1531 }
1532
1533 // accumalate for averaging
1534 for (int j = 0; j < this.psdNumberOfPoints; j++)
1535 avePSD[1][j] += segPSD[j];
1536
1537 // shift half segment
1538 for (int j = 0; j < this.segmentLength; j++) {
1539 data[j] = data[j + this.segmentLength];
1540 }
1541 }
1542 } else {
1543 // No overlap
1544
1545 // loop through segments
1546 for (int i = 1; i <= this.segmentNumber; i++) {
1547
1548 // read in segment data
1549 for (int j = 0; j < 2 * this.segmentLength; j++) {
1550 data[j] = fin.readDouble();
1551 }
1552
1553 // window data
1554 int k = -1;
1555 for (int j = 0; j < this.segmentLength; j++) {
1556 data[++k] = data[k] * this.weights[j];
1557 data[++k] = data[k] * this.weights[j];
1558 }
1559
1560 // perform fft on windowed segment
1561 basicFft(data, this.segmentLength, isign);
1562
1563 // obtain weighted mean square amplitudes
1564 segPSD[0] = Fmath.square(data[0]) + Fmath.square(data[1]);
1565 for (int j = 1; j < this.psdNumberOfPoints; j++) {
1566 segPSD[j] = Fmath.square(data[2 * j]) + Fmath.square(data[2 * j + 1])
1567 + Fmath.square(data[2 * this.segmentLength - 2 * j])
1568 + Fmath.square(data[2 * this.segmentLength - 2 * j + 1]);
1569 }
1570
1571 // Normalise
1572 for (int j = 1; j < this.psdNumberOfPoints; j++) {
1573 segPSD[j] = 2.0D * segPSD[j] / (this.segmentLength * this.sumOfSquaredWeights);
1574 }
1575
1576 // accumalate for averaging
1577 for (int j = 0; j < this.psdNumberOfPoints; j++)
1578 avePSD[1][j] += segPSD[j];
1579 }
1580 }
1581
1582 // average all segments
1583 for (int j = 0; j < this.psdNumberOfPoints; j++)
1584 avePSD[1][j] /= this.segmentNumber;
1585
1586 // Calculate frequencies
1587 for (int i = 0; i < this.psdNumberOfPoints; i++) {
1588 avePSD[0][i] = (double) i / ((double) this.segmentLength * this.deltaT);
1589 }
1590
1591 return avePSD;
1592 }
1593
1594 // Get the power spectrum
1595 public double[][] getpowerSpectrumEstimate() {
1596 if (!this.powSpecDone)
1597 System.out.println("getpowerSpectrumEstimate - powerSpectrum has not been called - null returned");
1598 return this.powerSpectrumEstimate;
1599 }
1600
1601 // get the number of power spectrum frequency points
1602 public int getNumberOfPsdPoints() {
1603 return this.psdNumberOfPoints;
1604 }
1605
1606 // Print the power spectrum to a text file
1607 // default file name
1608 public void printPowerSpectrum() {
1609 String filename = "FourierTransformPSD.txt";
1610 printPowerSpectrum(filename);
1611 }
1612
1613 // Print the power spectrum to a text file
1614 public void printPowerSpectrum(String filename) {
1615 if (!this.powSpecDone)
1616 this.powerSpectrum();
1617
1618 FileOutput fout = new FileOutput(filename);
1619 fout.println("Power Spectrum Density Estimate Output File from FourierTransform");
1620 fout.dateAndTimeln(filename);
1621 String title = "Window: " + this.windowNames[this.windowOption];
1622 if (this.windowOption == 6)
1623 title += ", alpha = " + this.kaiserAlpha;
1624 if (this.windowOption == 7)
1625 title += ", alpha = " + this.gaussianAlpha;
1626 fout.println(title);
1627 fout.printtab("Number of segments = ");
1628 fout.println(this.segmentNumber);
1629 fout.printtab("Segment length = ");
1630 fout.println(this.segmentLength);
1631 if (this.segmentNumber > 1) {
1632 if (this.overlap) {
1633 fout.printtab("Segments overlap by 50%");
1634 } else {
1635 fout.printtab("Segments do not overlap");
1636 }
1637 }
1638
1639 fout.println();
1640 printWarnings(fout);
1641
1642 fout.printtab("Frequency");
1643 fout.println("Mean Square");
1644 fout.printtab("(cycles per");
1645 fout.println("Amplitude");
1646 if (this.deltaTset) {
1647 fout.printtab("unit time)");
1648 } else {
1649 fout.printtab("gridpoint)");
1650 }
1651 fout.println(" ");
1652 int n = this.powerSpectrumEstimate[0].length;
1653 for (int i = 0; i < n; i++) {
1654 fout.printtab(Fmath.truncate(this.powerSpectrumEstimate[0][i], 4));
1655 fout.println(Fmath.truncate(this.powerSpectrumEstimate[1][i], 4));
1656 }
1657 fout.close();
1658 }
1659
1660 // Display a plot of the power spectrum from the given point number
1661 // no graph title provided
1662 public void plotPowerSpectrum(int lowPoint) {
1663 String graphTitle = "Estimation of Power Spectrum Density";
1664 this.plotPowerSpectrum(lowPoint, this.powerSpectrumEstimate[0].length - 1, graphTitle);
1665 }
1666
1667 // Display a plot of the power spectrum from the given point number
1668 // graph title provided
1669 public void plotPowerSpectrum(int lowPoint, String graphTitle) {
1670 this.plotPowerSpectrum(lowPoint, this.powerSpectrumEstimate[0].length - 1, graphTitle);
1671 }
1672
1673 // Display a plot of the power spectrum within a defined points window
1674 // no graph title provided
1675 public void plotPowerSpectrum(int lowPoint, int highPoint) {
1676 String graphTitle = "Estimation of Power Spectrum Density";
1677 this.plotPowerSpectrum(lowPoint, highPoint, graphTitle);
1678 }
1679
1680 // Display a plot of the power spectrum within a defined points window
1681 // Graph title provided
1682 public void plotPowerSpectrum(int lowPoint, int highPoint, String graphTitle) {
1683 if (!this.powSpecDone) {
1684 System.out.println("plotPowerSpectrum - powerSpectrum has not been called - no plot displayed");
1685 } else {
1686 int n = this.powerSpectrumEstimate[0].length - 1;
1687 if (lowPoint < 0 || lowPoint >= n)
1688 lowPoint = 0;
1689 if (highPoint < 0 || highPoint > n)
1690 highPoint = n;
1691 this.plotPowerSpectrumLinear(lowPoint, highPoint, graphTitle);
1692 }
1693 }
1694
1695 // Display a plot of the power spectrum from a given frequency
1696 // no graph title provided
1697 public void plotPowerSpectrum(double lowFreq) {
1698 String graphTitle = "Estimation of Power Spectrum Density";
1699 this.plotPowerSpectrum(lowFreq, graphTitle);
1700 }
1701
1702 // Display a plot of the power spectrum from a given frequency
1703 // graph title provided
1704 public void plotPowerSpectrum(double lowFreq, String graphTitle) {
1705 if (!this.powSpecDone)
1706 this.powerSpectrum();
1707
1708 double highFreq = this.powerSpectrumEstimate[1][this.powerSpectrumEstimate[0].length - 1];
1709 this.plotPowerSpectrum(lowFreq, highFreq, graphTitle);
1710 }
1711
1712 // Display a plot of the power spectrum within a defined frequency window
1713 // no graph title provided
1714 public void plotPowerSpectrum(double lowFreq, double highFreq) {
1715 if (!this.powSpecDone) {
1716 System.out.println("plotPowerSpectrum - powerSpectrum has not been called - no plot displayed");
1717 } else {
1718 String graphTitle = "Estimation of Power Spectrum Density";
1719 this.plotPowerSpectrum(lowFreq, highFreq, graphTitle);
1720 }
1721 }
1722
1723 // Display a plot of the power spectrum within a defined frequency window
1724 // graph title provided
1725 public void plotPowerSpectrum(double lowFreq, double highFreq, String graphTitle) {
1726 if (!this.powSpecDone) {
1727 System.out.println("plotPowerSpectrum - powerSpectrum has not been called - no plot displayed");
1728 } else {
1729 int low = 0;
1730 int high = 0;
1731 if (!this.deltaTset) {
1732 System.out.println("plotPowerSpectrum - deltaT has not been set");
1733 System.out.println("full spectrum plotted");
1734 } else {
1735 int ii = 0;
1736 int n = this.powerSpectrumEstimate[0].length - 1;
1737 boolean test = true;
1738 if (lowFreq == -1.0D) {
1739 low = 1;
1740 } else {
1741 while (test) {
1742 if (this.powerSpectrumEstimate[0][ii] > lowFreq) {
1743 low = ii - 1;
1744 if (low < 0)
1745 low = 0;
1746 test = false;
1747 } else {
1748 ii++;
1749 if (ii >= n) {
1750 low = 0;
1751 System.out.println("plotPowerSpectrum - lowFreq out of range - reset to zero");
1752 test = false;
1753 }
1754 }
1755 }
1756 }
1757 test = true;
1758 ii = 0;
1759 while (test) {
1760 if (this.powerSpectrumEstimate[0][ii] > highFreq) {
1761 high = ii - 1;
1762 if (high < 0) {
1763 System.out.println("plotPowerSpectrum - highFreq out of range - reset to highest value");
1764 high = n;
1765 }
1766 test = false;
1767 } else {
1768 ii++;
1769 if (ii >= n) {
1770 high = n;
1771 System.out.println("plotPowerSpectrum - highFreq out of range - reset to highest value");
1772 test = false;
1773 }
1774 }
1775 }
1776 this.plotPowerSpectrumLinear(low, high, graphTitle);
1777 }
1778 }
1779 }
1780
1781 // Display a plot of the power spectrum
1782 // no graph title provided
1783 public void plotPowerSpectrum() {
1784 if (!this.powSpecDone)
1785 this.powerSpectrum();
1786
1787 String graphTitle = "Estimation of Power Spectrum Density";
1788 this.plotPowerSpectrumLinear(0, this.powerSpectrumEstimate[0].length - 1, graphTitle);
1789 }
1790
1791 // Display a plot of the power spectrum
1792 public void plotPowerSpectrum(String graphTitle) {
1793 if (!this.powSpecDone)
1794 this.powerSpectrum();
1795
1796 this.plotPowerSpectrumLinear(0, this.powerSpectrumEstimate[0].length - 1, graphTitle);
1797 }
1798
1799 // Prepare a plot of the power spectrum (linear)
1800 private void plotPowerSpectrumLinear(int low, int high, String graphTitle) {
1801
1802 int nData = this.powerSpectrumEstimate[0].length;
1803 int nNew = high - low + 1;
1804 double[][] spectrum = new double[2][nNew];
1805 for (int i = 0; i < nNew; i++) {
1806 spectrum[0][i] = this.powerSpectrumEstimate[0][i + low];
1807 spectrum[1][i] = this.powerSpectrumEstimate[1][i + low];
1808 }
1809 String yLegend = "Mean Square Amplitude";
1810
1811 plotPowerDisplay(spectrum, low, high, graphTitle, yLegend);
1812 }
1813
1814 // Display a log plot of the power spectrum from the given point number
1815 // no graph title provided
1816 public void plotPowerLog(int lowPoint) {
1817 String graphTitle = "Estimation of Power Spectrum Density";
1818 this.plotPowerLog(lowPoint, this.powerSpectrumEstimate[0].length - 1, graphTitle);
1819 }
1820
1821 // Display a log plot of the power spectrum from the given point number
1822 // graph title provided
1823 public void plotPowerLog(int lowPoint, String graphTitle) {
1824 this.plotPowerLog(lowPoint, this.powerSpectrumEstimate[0].length - 1, graphTitle);
1825 }
1826
1827 // Display a log plot of the power spectrum within a defined points window
1828 // no graph title provided
1829 public void plotPowerLog(int lowPoint, int highPoint) {
1830 String graphTitle = "Estimation of Power Spectrum Density";
1831 this.plotPowerLog(lowPoint, highPoint, graphTitle);
1832 }
1833
1834 // Display a plot of the power spectrum within a defined points window
1835 // Graph title provided
1836 public void plotPowerLog(int lowPoint, int highPoint, String graphTitle) {
1837 if (!this.powSpecDone)
1838 this.powerSpectrum();
1839
1840 int n = this.powerSpectrumEstimate[0].length - 1;
1841 if (lowPoint < 0 || lowPoint >= n)
1842 lowPoint = 0;
1843 if (highPoint < 0 || highPoint > n)
1844 highPoint = n;
1845 this.plotPowerSpectrumLog(lowPoint, highPoint, graphTitle);
1846 }
1847
1848 // Display a plot of the power spectrum from a given frequency
1849 // no graph title provided
1850 public void plotPowerLog(double lowFreq) {
1851 String graphTitle = "Estimation of Power Spectrum Density";
1852 this.plotPowerLog(lowFreq, graphTitle);
1853 }
1854
1855 // Display a log plot of the power spectrum from a given frequency
1856 // graph title provided
1857 public void plotPowerLog(double lowFreq, String graphTitle) {
1858 if (!this.powSpecDone)
1859 this.powerSpectrum();
1860
1861 double highFreq = this.powerSpectrumEstimate[1][this.powerSpectrumEstimate[0].length - 1];
1862 this.plotPowerLog(lowFreq, highFreq, graphTitle);
1863 }
1864
1865 // Display a plot of the power spectrum within a defined frequency window
1866 // no graph title provided
1867 public void plotPowerLog(double lowFreq, double highFreq) {
1868 if (!this.powSpecDone)
1869 this.powerSpectrum();
1870
1871 String graphTitle = "Estimation of Power Spectrum Density";
1872 this.plotPowerLog(lowFreq, highFreq, graphTitle);
1873 }
1874
1875 // Display a log plot of the power spectrum within a defined frequency
1876 // window
1877 // graph title provided
1878 public void plotPowerLog(double lowFreq, double highFreq, String graphTitle) {
1879 if (!this.powSpecDone)
1880 this.powerSpectrum();
1881
1882 int low = 0;
1883 int high = 0;
1884 if (!this.deltaTset) {
1885 System.out.println("plotPowerLog - deltaT has not been set");
1886 System.out.println("full spectrum plotted");
1887 } else {
1888 int ii = 0;
1889 int n = this.powerSpectrumEstimate[0].length - 1;
1890 boolean test = true;
1891 if (lowFreq == -1.0D) {
1892 low = 1;
1893 } else {
1894 while (test) {
1895 if (this.powerSpectrumEstimate[0][ii] > lowFreq) {
1896 low = ii - 1;
1897 if (low < 0)
1898 low = 0;
1899 test = false;
1900 } else {
1901 ii++;
1902 if (ii >= n) {
1903 low = 0;
1904 System.out.println("plotPowerLog - lowFreq out of range - reset to zero");
1905 test = false;
1906 }
1907 }
1908 }
1909 }
1910 test = true;
1911 ii = 0;
1912 while (test) {
1913 if (this.powerSpectrumEstimate[0][ii] > highFreq) {
1914 high = ii - 1;
1915 if (high < 0) {
1916 System.out.println("plotPowerLog - highFreq out of range - reset to highest value");
1917 high = n;
1918 }
1919 test = false;
1920 } else {
1921 ii++;
1922 if (ii >= n) {
1923 high = n;
1924 System.out.println("plotPowerSpectrum - highFreq out of range - reset to highest value");
1925 test = false;
1926 }
1927 }
1928 }
1929 this.plotPowerSpectrumLog(low, high, graphTitle);
1930 }
1931 }
1932
1933 // Display a log plot of the power spectrum
1934 // no graph title provided
1935 public void plotPowerLog() {
1936 if (!this.powSpecDone)
1937 this.powerSpectrum();
1938
1939 String graphTitle = "Estimation of Power Spectrum Density";
1940 this.plotPowerSpectrumLog(0, this.powerSpectrumEstimate[0].length - 1, graphTitle);
1941 }
1942
1943 // Display a log plot of the power spectrum
1944 public void plotPowerLog(String graphTitle) {
1945 if (!this.powSpecDone)
1946 this.powerSpectrum();
1947
1948 this.plotPowerSpectrumLog(0, this.powerSpectrumEstimate[0].length - 1, graphTitle);
1949 }
1950
1951 // Prepare a plot of the power spectrum (log)
1952 private void plotPowerSpectrumLog(int low, int high, String graphTitle) {
1953
1954 int nData = this.powerSpectrumEstimate[0].length;
1955 int nNew = high - low + 1;
1956 double[][] spectrum = new double[2][nNew];
1957 for (int i = 0; i < nNew; i++) {
1958 spectrum[0][i] = this.powerSpectrumEstimate[0][i + low];
1959 spectrum[1][i] = this.powerSpectrumEstimate[1][i + low];
1960 }
1961
1962 // Find minimum of amplitudes that is not zero
1963 // find first non-zero value
1964 boolean test = true;
1965 int ii = 0;
1966 double minimum = 0.0D;
1967 while (test) {
1968 if (spectrum[1][ii] > 0.0D) {
1969 minimum = spectrum[1][ii];
1970 test = false;
1971 } else {
1972 ii++;
1973 if (ii >= nNew) {
1974 test = false;
1975 System.out.println("plotPowerSpectrumLog: no non-zero amplitudes");
1976 System.exit(0);
1977 }
1978 }
1979 }
1980
1981 // Find minimum
1982 for (int i = ii + 1; i < nNew; i++)
1983 if (spectrum[1][i] < minimum)
1984 minimum = spectrum[1][i];
1985
1986 // Replace zeros with minimum
1987 for (int i = 0; i < nNew; i++)
1988 if (spectrum[1][i] <= 0.0D)
1989 spectrum[1][i] = minimum;
1990
1991 // Take log to base 10
1992 for (int i = 0; i < nNew; i++)
1993 spectrum[1][i] = Fmath.log10(spectrum[1][i]);
1994
1995 // call display method
1996 String yLegend = "Log10(Mean Square Amplitude)";
1997 plotPowerDisplay(spectrum, low, high, graphTitle, yLegend);
1998 }
1999
2000 // Display a plot of the power spectrum
2001 private void plotPowerDisplay(double[][] spectrum, int low, int high, String graphTitle, String yLegend) {
2002
2003 PlotGraph pg = new PlotGraph(spectrum);
2004 graphTitle = graphTitle + " [plot between points " + low + " and " + high + "]";
2005 pg.setGraphTitle(graphTitle);
2006 String graphTitle2 = "Window: " + this.windowNames[this.windowOption];
2007 if (this.windowOption == 6)
2008 graphTitle2 += " - alpha = " + this.kaiserAlpha;
2009 if (this.windowOption == 7)
2010 graphTitle2 += " - alpha = " + this.gaussianAlpha;
2011 graphTitle2 += ", " + this.segmentNumber + " segment/s of length " + this.segmentLength;
2012 if (this.segmentNumber > 1) {
2013 if (this.overlap) {
2014 graphTitle2 += ", segments overlap by 50%";
2015 } else {
2016 graphTitle2 += ", segments do not overlap";
2017 }
2018 }
2019
2020 pg.setGraphTitle2(graphTitle2);
2021 pg.setXaxisLegend("Frequency");
2022 if (this.deltaTset) {
2023 pg.setXaxisUnitsName("cycles per unit time");
2024 } else {
2025 pg.setXaxisUnitsName("cycles per grid point");
2026 }
2027 pg.setYaxisLegend(yLegend);
2028
2029 switch (this.plotLineOption) {
2030 case 0:
2031 pg.setLine(3);
2032 break;
2033 case 1:
2034 pg.setLine(1);
2035 break;
2036 case 2:
2037 pg.setLine(2);
2038 break;
2039 default:
2040 pg.setLine(3);
2041 }
2042
2043 switch (this.plotPointOption) {
2044 case 0:
2045 pg.setPoint(0);
2046 break;
2047 case 1:
2048 pg.setPoint(4);
2049 break;
2050 default:
2051 pg.setPoint(0);
2052 }
2053
2054 pg.plot();
2055
2056 }
2057
2058 // Set the line option in plotting the power spectrum or correlation
2059 // = 0 join points with straight lines
2060 // = 1 cubic spline interpolation
2061 // = 3 no line - only points
2062 public void setPlotLineOption(int lineOpt) {
2063 this.plotLineOption = lineOpt;
2064 }
2065
2066 // Get the line option in ploting the power spectrum or correlation
2067 // = 0 join points with straight lines
2068 // = 1 cubic spline interpolation
2069 // = 3 no line - only points
2070 public int getPlotLineOption() {
2071 return this.plotLineOption;
2072 }
2073
2074 // Set the point option in plotting the power spectrum or correlation
2075 // = 0 no point symbol
2076 // = 1 filled circles
2077 public void setPlotPointOption(int pointOpt) {
2078 this.plotPointOption = pointOpt;
2079 }
2080
2081 // Get the point option in plotting the power spectrum or correlation
2082 // = 0 no point symbol
2083 // = 1 filled circles
2084 public int getPlotPointOption() {
2085 return this.plotPointOption;
2086 }
2087
2088 // Return correlation of data already entered with data passed as this
2089 // method's argument
2090 // data must be real
2091 public double[][] correlate(double[] data) {
2092 int nLen = data.length;
2093 if (!this.fftDataSet)
2094 throw new IllegalArgumentException("No data has been previously entered");
2095 if (nLen != this.originalDataLength)
2096 throw new IllegalArgumentException("The two data sets to be correlated are of different length");
2097 if (!FourierTransform.checkPowerOfTwo(nLen))
2098 throw new IllegalArgumentException(
2099 "The length of the correlation data sets is not equal to an integer power of two");
2100
2101 this.complexCorr = Complex.oneDarray(nLen);
2102 for (int i = 0; i < nLen; i++) {
2103 this.complexCorr[i].setReal(data[i]);
2104 this.complexCorr[i].setImag(0.0D);
2105 }
2106
2107 this.fftCorr = new double[2 * nLen];
2108 int j = -1;
2109 for (int i = 0; i < nLen; i++) {
2110 this.fftCorr[++j] = data[i];
2111 this.fftCorr[++j] = 0.0D;
2112 }
2113
2114 return correlation(nLen);
2115 }
2116
2117 // Return correlation of data1 and data2 passed as this method's arguments
2118 // data must be real
2119 public double[][] correlate(double[] data1, double[] data2) {
2120 int nLen = data1.length;
2121 int nLen2 = data2.length;
2122 if (nLen != nLen2)
2123 throw new IllegalArgumentException("The two data sets to be correlated are of different length");
2124 if (!FourierTransform.checkPowerOfTwo(nLen))
2125 throw new IllegalArgumentException(
2126 "The length of the correlation data sets is not equal to an integer power of two");
2127
2128 this.fftDataLength = nLen;
2129 this.complexData = Complex.oneDarray(this.fftDataLength);
2130 for (int i = 0; i < this.fftDataLength; i++) {
2131 this.complexData[i].setReal(data1[i]);
2132 this.complexData[i].setImag(0.0D);
2133 }
2134
2135 this.fftData = new double[2 * this.fftDataLength];
2136 int j = 0;
2137 for (int i = 0; i < this.fftDataLength; i++) {
2138 this.fftData[j] = data1[i];
2139 j++;
2140 this.fftData[j] = 0.0D;
2141 j++;
2142 }
2143 this.fftDataSet = true;
2144
2145 this.fftDataWindow = new double[2 * this.fftDataLength];
2146 this.weights = new double[this.fftDataLength];
2147 this.sumOfSquaredWeights = windowData(this.fftData, this.fftDataWindow, this.weights);
2148
2149 this.transformedDataFft = new double[2 * this.fftDataLength];
2150 this.transformedDataComplex = Complex.oneDarray(this.fftDataLength);
2151
2152 this.complexCorr = Complex.oneDarray(nLen);
2153 for (int i = 0; i < nLen; i++) {
2154 this.complexCorr[i].setReal(data2[i]);
2155 this.complexCorr[i].setImag(0.0D);
2156 }
2157
2158 this.fftCorr = new double[2 * nLen];
2159 j = -1;
2160 for (int i = 0; i < nLen; i++) {
2161 this.fftCorr[++j] = data2[i];
2162 this.fftCorr[++j] = 0.0D;
2163 }
2164
2165 return correlation(nLen);
2166 }
2167
2168 // Returns the correlation of the data in fftData and fftCorr
2169 private double[][] correlation(int nLen) {
2170
2171 this.fftDataWindow = new double[2 * nLen];
2172 this.fftCorrWindow = new double[2 * nLen];
2173 this.weights = new double[nLen];
2174
2175 this.sumOfSquaredWeights = windowData(this.fftData, this.fftDataWindow, this.weights);
2176 windowData(this.fftCorr, this.fftCorrWindow, this.weights);
2177
2178 // Perform fft on first set of stored data
2179 int isign = 1;
2180 double[] hold1 = new double[2 * nLen];
2181 for (int i = 0; i < nLen * 2; i++)
2182 hold1[i] = this.fftDataWindow[i];
2183 basicFft(hold1, nLen, isign);
2184
2185 // Perform fft on second set of stored data
2186 isign = 1;
2187 double[] hold2 = new double[2 * nLen];
2188 for (int i = 0; i < nLen * 2; i++)
2189 hold2[i] = this.fftCorrWindow[i];
2190 basicFft(hold2, nLen, isign);
2191
2192 // multiply hold1 by complex congugate of hold2
2193 double[] hold3 = new double[2 * nLen];
2194 int j = 0;
2195 for (int i = 0; i < nLen; i++) {
2196 hold3[j] = (hold1[j] * hold2[j] + hold1[j + 1] * hold2[j + 1]) / nLen;
2197 hold3[j + 1] = (-hold1[j] * hold2[j + 1] + hold1[j + 1] * hold2[j]) / nLen;
2198 j += 2;
2199 }
2200
2201 // Inverse transform -> correlation
2202 isign = -1;
2203 basicFft(hold3, nLen, isign);
2204
2205 // fill correlation array
2206 for (int i = 0; i < 2 * nLen; i++)
2207 this.transformedDataFft[i] = hold3[i];
2208 this.correlationArray = new double[2][nLen];
2209 j = 0;
2210 int k = nLen;
2211 for (int i = nLen / 2 + 1; i < nLen; i++) {
2212 this.correlationArray[1][j] = hold3[k] / nLen;
2213 j++;
2214 k += 2;
2215 }
2216 k = 0;
2217 for (int i = 0; i < nLen / 2; i++) {
2218 this.correlationArray[1][j] = hold3[k] / nLen;
2219 j++;
2220 k += 2;
2221 }
2222
2223 // calculate time lags
2224 this.correlationArray[0][0] = -(double) (nLen / 2) * this.deltaT;
2225 for (int i = 1; i < nLen; i++) {
2226 this.correlationArray[0][i] = this.correlationArray[0][i - 1] + this.deltaT;
2227 }
2228
2229 this.correlateDone = true;
2230 return this.correlationArray;
2231 }
2232
2233 // Get the correlation
2234 public double[][] getCorrelation() {
2235 if (!this.correlateDone) {
2236 System.out.println("getCorrelation - correlation has not been called - no correlation returned");
2237 }
2238 return this.correlationArray;
2239 }
2240
2241 // Print the correlation to a text file
2242 // default file name
2243 public void printCorrelation() {
2244 String filename = "Correlation.txt";
2245 printCorrelation(filename);
2246 }
2247
2248 // Print the correlation to a text file
2249 public void printCorrelation(String filename) {
2250 if (!this.correlateDone) {
2251 System.out.println("printCorrelation - correlate has not been called - no file printed");
2252 } else {
2253 FileOutput fout = new FileOutput(filename);
2254 fout.println("Correlation Output File from FourierTransform");
2255 fout.dateAndTimeln(filename);
2256 String title = "Window: " + this.windowNames[this.windowOption];
2257 if (this.windowOption == 6)
2258 title += ", alpha = " + this.kaiserAlpha;
2259 if (this.windowOption == 7)
2260 title += ", alpha = " + this.gaussianAlpha;
2261 fout.println(title);
2262 fout.printtab("Data length = ");
2263 fout.println(this.fftDataLength);
2264 fout.println();
2265
2266 fout.printtab("Time lag");
2267 fout.println("Correlation");
2268 if (this.deltaTset) {
2269 fout.printtab("/unit time");
2270 } else {
2271 fout.printtab("/grid interval)");
2272 }
2273 fout.println("Coefficient");
2274
2275 int n = this.correlationArray[0].length;
2276 for (int i = 0; i < n; i++) {
2277 fout.printtab(Fmath.truncate(this.correlationArray[0][i], 4));
2278 fout.println(Fmath.truncate(this.correlationArray[1][i], 4));
2279 }
2280 fout.close();
2281 }
2282 }
2283
2284 // Display a plot of the correlation
2285 // no graph title provided
2286 public void plotCorrelation() {
2287 if (!this.correlateDone) {
2288 System.out.println("plotCorrelation - correlation has not been called - no plot displayed");
2289 } else {
2290 String graphTitle = "Correlation Plot";
2291 plotCorrelation(graphTitle);
2292 }
2293 }
2294
2295 // Display a plot of the correlation
2296 public void plotCorrelation(String graphTitle) {
2297 if (!this.correlateDone) {
2298 System.out.println("plotCorrelation - correlate has not been called - no plot displayed");
2299 } else {
2300
2301 PlotGraph pg = new PlotGraph(this.correlationArray);
2302 pg.setGraphTitle(graphTitle);
2303 String graphTitle2 = "Window: " + this.windowNames[this.windowOption];
2304 if (this.windowOption == 6)
2305 graphTitle2 += " - alpha = " + this.kaiserAlpha;
2306 if (this.windowOption == 7)
2307 graphTitle2 += " - alpha = " + this.gaussianAlpha;
2308
2309 pg.setGraphTitle2(graphTitle2);
2310 pg.setXaxisLegend("Correlation Lag");
2311 if (this.deltaTset) {
2312 pg.setXaxisUnitsName("unit time");
2313 } else {
2314 pg.setXaxisUnitsName("grid interval");
2315 }
2316 pg.setYaxisLegend("Correlation coefficient");
2317
2318 switch (this.plotLineOption) {
2319 case 0:
2320 pg.setLine(3);
2321 break;
2322 case 1:
2323 pg.setLine(1);
2324 break;
2325 case 2:
2326 pg.setLine(2);
2327 break;
2328 default:
2329 pg.setLine(3);
2330 }
2331
2332 switch (this.plotPointOption) {
2333 case 0:
2334 pg.setPoint(0);
2335 break;
2336 case 1:
2337 pg.setPoint(4);
2338 break;
2339 default:
2340 pg.setPoint(0);
2341 }
2342
2343 pg.plot();
2344 }
2345 }
2346
2347 // Performs a fft power spectrum density (psd) estimation
2348 // on a moving window throughout the original data set
2349 // returning the results as a frequency time matrix
2350 // windowLength is the length of the window in time units
2351 public double[][] shortTime(double windowTime) {
2352 int windowLength = (int) Math.round(windowTime / this.deltaT);
2353 if (!this.checkPowerOfTwo(windowLength)) {
2354 int low = this.lastPowerOfTwo(windowLength);
2355 int high = this.nextPowerOfTwo(windowLength);
2356
2357 if ((windowLength - low) <= (high - windowLength)) {
2358 windowLength = low;
2359 if (low == 0)
2360 windowLength = high;
2361 } else {
2362 windowLength = high;
2363 }
2364 System.out.println("Method - shortTime");
2365 System.out.println("Window length, provided as time, " + windowTime
2366 + ", did not convert to an integer power of two data points");
2367 System.out.println("A value of " + ((windowLength - 1) * this.deltaT) + " was substituted");
2368 }
2369
2370 return shortTime(windowLength);
2371 }
2372
2373 // Performs a fft power spectrum density (psd) estimation
2374 // on a moving window throughout the original data set
2375 // returning the results as a frequency time matrix
2376 // windowLength is the number of points in the window
2377 public double[][] shortTime(int windowLength) {
2378
2379 if (!FourierTransform.checkPowerOfTwo(windowLength))
2380 throw new IllegalArgumentException(
2381 "Moving window data length ," + windowLength + ", is not an integer power of two");
2382 if (!this.fftDataSet)
2383 throw new IllegalArgumentException("No data has been entered for the Fast Fourier Transform");
2384 if (windowLength > this.originalDataLength)
2385 throw new IllegalArgumentException("The window length, " + windowLength
2386 + ", is greater than the data length, " + this.originalDataLength + ".");
2387
2388 // if no window option has been set - default = Gaussian with alpha =
2389 // 2.5
2390 if (this.windowOption == 0)
2391 this.setGaussian();
2392 // set up time-frequency matrix
2393 // first row = blank cell followed by time vector
2394 // first column = blank cell followed by frequency vector
2395 // each cell is then the mean square amplitude at that frequency and
2396 // time
2397 this.numShortTimes = this.originalDataLength - windowLength + 1;
2398 this.numShortFreq = windowLength / 2;
2399 this.timeFrequency = new double[this.numShortFreq + 1][this.numShortTimes + 1];
2400 this.timeFrequency[0][0] = 0.0D;
2401 this.timeFrequency[0][1] = (double) (windowLength - 1) * this.deltaT / 2.0D;
2402 for (int i = 2; i <= this.numShortTimes; i++) {
2403 this.timeFrequency[0][i] = this.timeFrequency[0][i - 1] + this.deltaT;
2404 }
2405 for (int i = 0; i < this.numShortFreq; i++) {
2406 this.timeFrequency[i + 1][0] = (double) i / ((double) windowLength * this.deltaT);
2407 }
2408
2409 // set up window details
2410 this.segmentLength = windowLength;
2411 int windowStartIndex = 0;
2412 double[] data = new double[2 * windowLength]; // holds data and
2413 // transformed data for
2414 // working window
2415 double[] winPSD = new double[this.numShortFreq]; // holds psd for
2416 // working window
2417 int isign = 1;
2418
2419 // loop through time shifts
2420 for (int i = 1; i <= this.numShortTimes; i++) {
2421
2422 // collect window data
2423 for (int j = 0; j < 2 * windowLength; j++)
2424 data[j] = this.fftData[windowStartIndex + j];
2425
2426 // window data
2427 if (i == 1) {
2428 this.sumOfSquaredWeights = this.windowData(data, data, this.weights);
2429 } else {
2430 int k = 0;
2431 for (int j = 0; j < this.segmentLength; j++) {
2432 data[k] = data[k] * this.weights[j];
2433 data[++k] = data[k] * this.weights[j];
2434 ++k;
2435 }
2436 }
2437
2438 // perform fft on windowed segment
2439 basicFft(data, windowLength, isign);
2440
2441 // obtain weighted mean square amplitudes
2442 winPSD[0] = Fmath.square(data[0]) + Fmath.square(data[1]);
2443 for (int j = 1; j < this.numShortFreq; j++) {
2444 winPSD[j] = Fmath.square(data[2 * j]) + Fmath.square(data[2 * j + 1])
2445 + Fmath.square(data[2 * windowLength - 2 * j])
2446 + Fmath.square(data[2 * windowLength - 2 * j + 1]);
2447 }
2448
2449 // Normalise and place in time-frequency matrix
2450 for (int j = 0; j < this.numShortFreq; j++) {
2451 timeFrequency[j + 1][i] = 2.0D * winPSD[j] / (windowLength * this.sumOfSquaredWeights);
2452 }
2453
2454 // increment segment start index
2455 windowStartIndex += 2;
2456 }
2457
2458 this.shortTimeDone = true;
2459 return this.timeFrequency;
2460 }
2461
2462 // Return time frequency matrix
2463 public double[][] getTimeFrequencyMatrix() {
2464 if (!this.shortTimeDone)
2465 throw new IllegalArgumentException("No short time Fourier transform has been performed");
2466 return this.timeFrequency;
2467 }
2468
2469 // Return number of times in short time Fourier transform
2470 public int getShortTimeNumberOfTimes() {
2471 if (!this.shortTimeDone)
2472 throw new IllegalArgumentException("No short time Fourier transform has been performed");
2473 return this.numShortTimes;
2474 }
2475
2476 // Return number of frequencies in short time Fourier transform
2477 public int getShortTimeNumberOfFrequencies() {
2478 if (!this.shortTimeDone)
2479 throw new IllegalArgumentException("No short time Fourier transform has been performed");
2480 return this.numShortFreq;
2481 }
2482
2483 // Return number of points in short time Fourier transform window
2484 public int getShortTimeWindowLength() {
2485 if (!this.shortTimeDone)
2486 throw new IllegalArgumentException("No short time Fourier transform has been performed");
2487 return this.segmentLength;
2488 }
2489
2490 // Print the short time Fourier transform to a text file
2491 // default file name
2492 public void printShortTime() {
2493 String filename = "ShortTime.txt";
2494 printShortTime(filename);
2495 }
2496
2497 // Print the short time Fourier transform to a text file
2498 public void printShortTime(String filename) {
2499 if (!this.shortTimeDone) {
2500 System.out.println("printShortTime- shortTime has not been called - no file printed");
2501 } else {
2502 FileOutput fout = new FileOutput(filename);
2503 fout.println("Short Time Fourier Transform Output File from FourierTransform");
2504 fout.dateAndTimeln(filename);
2505 String title = "Window: " + this.windowNames[this.windowOption];
2506 if (this.windowOption == 6)
2507 title += ", alpha = " + this.kaiserAlpha;
2508 if (this.windowOption == 7)
2509 title += ", alpha = " + this.gaussianAlpha;
2510 fout.println(title);
2511 fout.printtab("Data length = ");
2512 fout.println(this.originalDataLength);
2513 fout.printtab("Delta T = ");
2514 fout.println(this.deltaT);
2515 fout.printtab("Window length (points) = ");
2516 fout.println(this.segmentLength);
2517 fout.printtab("Window length (time units) = ");
2518 fout.println((this.segmentLength - 1) * this.deltaT);
2519 fout.printtab("Number of frequency points = ");
2520 fout.println(this.numShortFreq);
2521 fout.printtab("Number of time points = ");
2522 fout.println(this.numShortTimes);
2523
2524 // Average points if output would be greater than a text file line
2525 // length
2526 boolean checkAve = false;
2527 int newTp = this.numShortTimes;
2528 int maxN = 100;
2529 int nAve = this.numShortTimes / maxN;
2530 int nLast = this.numShortTimes % maxN;
2531 if (this.numShortTimes > 127) {
2532 checkAve = true;
2533 if (nLast > 0) {
2534 nAve += 1;
2535 newTp = maxN;
2536 nLast = this.numShortTimes - nAve * (newTp - 1);
2537 } else {
2538 newTp = maxN;
2539 nLast = nAve;
2540 }
2541 if (nLast != nAve) {
2542 fout.println("In the output below, each of the first " + (newTp - 2)
2543 + " magnitude points, along the time axis, is the average of " + nAve
2544 + " calculated points");
2545 fout.println("The last point is the average of " + nLast + " calculated points");
2546 } else {
2547 fout.println("In the output below, each magnitude point is the average of " + nAve
2548 + " calculated points");
2549 }
2550 fout.println("The data, without averaging, may be accessed using the method getTimeFrequencyMatrix()");
2551 }
2552 fout.println();
2553
2554 fout.println("first row = times");
2555 fout.println("first column = frequencies");
2556 fout.println("all other cells = mean square amplitudes at the corresponding time and frequency");
2557 if (checkAve) {
2558 double sum = 0.0D;
2559 int start = 1;
2560 int workingAve = nAve;
2561 for (int i = 0; i <= this.numShortFreq; i++) {
2562 fout.printtab(Fmath.truncate(this.timeFrequency[i][0], 4));
2563 start = 1;
2564 for (int j = 1; j <= newTp; j++) {
2565 workingAve = nAve;
2566 if (j == newTp)
2567 workingAve = nLast;
2568 sum = 0.0D;
2569 for (int k = start; k <= (start + workingAve - 1); k++) {
2570 sum += this.timeFrequency[i][k];
2571 }
2572 sum /= workingAve;
2573 fout.printtab(Fmath.truncate(sum, 4));
2574 start += workingAve;
2575 }
2576 fout.println();
2577 }
2578 } else {
2579 for (int i = 0; i <= this.numShortFreq; i++) {
2580 for (int j = 0; j <= newTp; j++) {
2581 fout.printtab(Fmath.truncate(this.timeFrequency[i][j], 4));
2582 }
2583 fout.println();
2584 }
2585 }
2586 fout.close();
2587 }
2588 }
2589
2590 // The paint method to draw the graph for plotShortTime.
2591 public void paint(Graphics g) {
2592
2593 // Call graphing method
2594 graph(g);
2595 }
2596
2597 // Set up the window and show graph for short time Fourier transform
2598 // user provided graph title
2599 public void plotShortTime(String title) {
2600 this.shortTitle = title;
2601 plotShortTime();
2602 }
2603
2604 // Set up the window and show graph for short time Fourier transform
2605 // No user provided graph title
2606 public void plotShortTime() {
2607 // Create the window object
2608 JFrame window = new JFrame("Michael T Flanagan's plotting program - FourierTransform.plotShortTime");
2609
2610 // Set the initial size of the graph window
2611 setSize(800, 600);
2612
2613 // Set background colour
2614 window.getContentPane().setBackground(Color.white);
2615
2616 // Choose close box
2617 window.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
2618
2619 // Add graph canvas
2620 window.getContentPane().add("Center", this);
2621
2622 // Set the window up
2623 window.pack();
2624 window.setResizable(true);
2625 window.toFront();
2626
2627 // Show the window
2628 window.setVisible(true);
2629 }
2630
2631 // graph method for plotShortTime short time Fourier Transform as a contour
2632 // plot
2633 public void graph(Graphics g) {
2634
2635 // graph axes positions
2636 int xLen = 512;
2637 int yLen = 256;
2638 int yTop = 100;
2639 int xBot = 100;
2640 int numBands = 18;
2641 // colours for contour map
2642 Color[] color = new Color[numBands + 1];
2643 color[18] = Color.black;
2644 color[17] = Color.darkGray;
2645 color[16] = Color.gray;
2646 color[15] = Color.lightGray;
2647 color[14] = Color.red.darker();
2648 color[13] = Color.red;
2649 color[12] = Color.magenta.darker();
2650 color[11] = Color.magenta;
2651 color[10] = Color.pink;
2652 color[9] = Color.pink.darker();
2653 color[8] = Color.orange.darker();
2654 color[7] = Color.orange;
2655 color[6] = Color.yellow;
2656 color[5] = Color.green;
2657 color[4] = Color.green.darker();
2658 color[3] = Color.cyan;
2659 color[2] = Color.cyan.darker();
2660 color[1] = Color.blue;
2661 color[0] = Color.blue.darker();
2662
2663 // Check and set parameters in case need to average or expand to match
2664 // fixed x-axis pixels
2665 int pixelsPerXpoint = 0;
2666 int xTp = 0;
2667 int xAve = 0;
2668 int xLast = 0;
2669 boolean xCheck = true;
2670 if (this.numShortTimes <= xLen) {
2671 pixelsPerXpoint = xLen / this.numShortTimes;
2672 xLen = pixelsPerXpoint * this.numShortTimes;
2673 xTp = this.numShortTimes;
2674 } else {
2675 xCheck = false;
2676 pixelsPerXpoint = 1;
2677 xTp = this.numShortTimes;
2678 xAve = this.numShortTimes / xLen;
2679 xLast = this.numShortTimes % xLen;
2680 if (xLast > 0) {
2681 xAve += 1;
2682 xTp = this.numShortTimes / xAve + 1;
2683 xLast = this.numShortTimes - xAve * (xTp - 1);
2684 } else {
2685 xTp = this.numShortTimes / xAve;
2686 xLast = xAve;
2687 }
2688 xLen = xTp;
2689 }
2690
2691 // Check and set parameters in case need to average or expand to match
2692 // fixed y-axis pixels
2693 int pixelsPerYpoint = 0;
2694 int yTp = 0;
2695 int yAve = 0;
2696 int yLast = 0;
2697 boolean yCheck = true;
2698
2699 if (this.numShortFreq <= yLen) {
2700 pixelsPerYpoint = yLen / this.numShortFreq;
2701 yLen = pixelsPerYpoint * this.numShortFreq;
2702 yTp = this.numShortFreq;
2703 } else {
2704 yCheck = false;
2705 pixelsPerYpoint = 1;
2706 yTp = this.numShortFreq;
2707 yAve = this.numShortFreq / yLen;
2708 yLast = this.numShortFreq % yLen;
2709 if (yLast > 0) {
2710 yAve += 1;
2711 yTp = this.numShortFreq / yAve + 1;
2712 yLast = this.numShortFreq - yAve * (yTp - 1);
2713 } else {
2714 yTp = this.numShortFreq / yAve;
2715 yLast = yAve;
2716 }
2717 yLen = yTp;
2718 }
2719
2720 // Complete axes positions
2721 int yBot = yTop + yLen;
2722 int xTop = xBot + xLen;
2723
2724 // declare contour map arrays
2725 double[][] averages = new double[yTp][xTp];
2726 int[][] pixels = new int[yTp][xTp];
2727 double[] times = new double[xTp];
2728 int[] timesPixels = new int[xTp];
2729 double[] freqs = new double[yTp];
2730 int[] freqPixels = new int[yTp];
2731
2732 double[][] hold = new double[this.numShortFreq][xTp];
2733
2734 // If necessary average or expand to match fixed y-axis pixels
2735 if (xCheck) {
2736 for (int i = 0; i <= this.numShortFreq; i++) {
2737 for (int j = 1; j <= this.numShortTimes; j++) {
2738 if (i == 0) {
2739 times[j - 1] = this.timeFrequency[0][j];
2740 } else {
2741 hold[i - 1][j - 1] = this.timeFrequency[i][j];
2742 }
2743 }
2744 }
2745 } else {
2746 double sum = 0.0D;
2747 int start = 1;
2748 int workingAve = xAve;
2749 for (int i = 0; i <= this.numShortFreq; i++) {
2750 start = 1;
2751 for (int j = 1; j <= xTp; j++) {
2752 workingAve = xAve;
2753 if (j == xTp)
2754 workingAve = xLast;
2755 sum = 0.0D;
2756 for (int k = start; k <= (start + workingAve - 1); k++) {
2757 sum += this.timeFrequency[i][k];
2758 }
2759 if (i == 0) {
2760 times[j - 1] = sum / workingAve;
2761 } else {
2762 hold[i - 1][j - 1] = sum / workingAve;
2763 }
2764 start += workingAve;
2765 }
2766 }
2767 }
2768
2769 // If necessary average or expand to match fixed x-axis pixels
2770 if (yCheck) {
2771 for (int i = 0; i < this.numShortFreq; i++) {
2772 freqs[i] = this.timeFrequency[i + 1][0];
2773 for (int j = 0; j < xTp; j++) {
2774 averages[i][j] = hold[i][j];
2775 }
2776 }
2777 } else {
2778 double sum = 0.0D;
2779 double sFreq = 0.0D;
2780 int start = 0;
2781 int workingAve = yAve;
2782 for (int i = 0; i < xTp; i++) {
2783 start = 0;
2784 for (int j = 0; j < yTp; j++) {
2785 workingAve = yAve;
2786 if (j == yTp - 1)
2787 workingAve = yLast;
2788 sum = 0.0D;
2789 sFreq = 0.0D;
2790 for (int k = start; k <= (start + workingAve - 1); k++) {
2791 sum += hold[k][i];
2792 sFreq += this.timeFrequency[k + 1][0];
2793 }
2794 averages[j][i] = sum;
2795 freqs[j] = sFreq / workingAve;
2796 start += workingAve;
2797 }
2798 }
2799 }
2800
2801 // Calculate contour bands
2802 double max = averages[0][0];
2803 double min = max;
2804 for (int i = 0; i < yTp; i++) {
2805 for (int j = 0; j < xTp; j++) {
2806 if (averages[i][j] > max)
2807 max = averages[i][j];
2808 if (averages[i][j] < min)
2809 min = averages[i][j];
2810 }
2811 }
2812
2813 double bandZero = 0.0D;
2814 if (min > 0.1D * max)
2815 bandZero = 0.99D * min;
2816 double bandWidth = (1.01D * max - 0.99D * min) / numBands;
2817 double[] band = new double[numBands];
2818 band[0] = bandZero + bandWidth;
2819 for (int i = 1; i < numBands; i++) {
2820 band[i] = band[i - 1] + bandWidth;
2821 }
2822 boolean test = true;
2823 for (int i = 0; i < yTp; i++) {
2824 for (int j = 0; j < xTp; j++) {
2825 test = true;
2826 int k = 0;
2827 while (test) {
2828 if (averages[i][j] <= band[k]) {
2829 pixels[i][j] = k;
2830 test = false;
2831 } else {
2832 k++;
2833 }
2834 }
2835 }
2836 }
2837
2838 // Plot contour coloured bands
2839 int yPixels = 0;
2840 int xPixels = 0;
2841 int yInner = 0;
2842 int xInner = 0;
2843 int xx = xBot;
2844 int yy = yTop;
2845 for (int i = 0; i < yTp; i++) {
2846 for (int j = 0; j < xTp; j++) {
2847 yInner = 0;
2848 for (int k = 0; k < pixelsPerYpoint; k++) {
2849 xInner = 0;
2850 for (int l = 0; l < pixelsPerXpoint; l++) {
2851 g.setColor(color[pixels[i][j]]);
2852 xx = xBot + (xPixels + xInner);
2853 yy = yBot - (yPixels + yInner);
2854 g.drawLine(xx, yy, xx, yy);
2855 xInner++;
2856 }
2857 yInner++;
2858 }
2859 xPixels += xInner;
2860 }
2861 yPixels += yInner;
2862 xPixels = 0;
2863 }
2864
2865 // draw axes
2866 g.setColor(color[numBands]);
2867 g.drawLine(xBot, yBot, xBot, yTop);
2868 g.drawLine(xTop, yBot, xTop, yTop);
2869 g.drawLine(xBot, yBot, xTop, yBot);
2870 g.drawLine(xBot, yTop, xTop, yTop);
2871
2872 // calculate axis legends and units
2873 int yInc = yLen / 4;
2874 int yScale = this.numShortFreq / 4;
2875 double yUnits = yInc * (freqs[1] - freqs[0]) / (pixelsPerYpoint * yScale);
2876 String[] yArray = new String[5];
2877 int yArr = 0;
2878 yArray[0] = "0 ";
2879 for (int i = 1; i < 5; i++) {
2880 yArr += yScale;
2881 yArray[i] = yArr + " ";
2882 }
2883 xx = xBot;
2884 yy = yBot;
2885 int yWord = 6 * (yArray[4].length() + 1);
2886 for (int i = 0; i < 5; i++) {
2887 g.drawLine(xx - 5, yy, xx, yy);
2888 g.drawString(yArray[i], xx - yWord, yy + 4);
2889 yy -= yInc;
2890 }
2891
2892 int xInc = xLen / 8;
2893 int xScale = this.numShortTimes / 8;
2894 double xUnits = xInc * (times[1] - times[0]) / (pixelsPerXpoint * xScale);
2895 String[] xArray = new String[9];
2896 int xArr = 0;
2897 xArray[0] = "0 ";
2898 for (int i = 1; i < 9; i++) {
2899 xArr += xScale;
2900 xArray[i] = xArr + " ";
2901 }
2902 xx = xBot;
2903 yy = yBot;
2904 for (int i = 0; i < 9; i++) {
2905 g.drawLine(xx, yy, xx, yy + 5);
2906 g.drawString(xArray[i], xx - 4, yy + 20);
2907 xx += xInc;
2908 }
2909
2910 // write graph and axis legends and units
2911 g.drawString("Short Time Fourier Transfer Time-Frequency Plot", xBot - 80, yTop - 80);
2912 g.drawString(this.shortTitle, xBot - 80, yTop - 60);
2913
2914 String yAxis = "Frequency / (" + Fmath.truncate(yUnits, 3) + " cycles per time unit)";
2915 g.drawString(yAxis, xBot - 60, yTop - 20);
2916 String xAxis = "Time / (" + Fmath.truncate(xUnits, 3) + " time units)";
2917 g.drawString(xAxis, xBot, yBot + 40);
2918 String totalTime = "Total time = " + (Fmath.truncate((xLen * (times[1] - times[0])) / pixelsPerXpoint, 3))
2919 + " time units";
2920 g.drawString(totalTime, xBot, yBot + 80);
2921
2922 String totalFreq = "Frequecy range = 0 to "
2923 + (Fmath.truncate((yLen * (freqs[1] - freqs[0])) / pixelsPerYpoint, 3)) + " cycles per time unit";
2924 g.drawString(totalFreq, xBot, yBot + 100);
2925
2926 g.drawString("Widow length = " + Fmath.truncate((this.segmentLength - 1) * this.deltaT, 3) + " time units",
2927 xBot, yBot + 120);
2928 String filter = "Window filter = " + this.windowNames[this.windowOption];
2929 if (this.windowOption == 6)
2930 filter += ", alpha = " + this.kaiserAlpha;
2931 if (this.windowOption == 7)
2932 filter += ", alpha = " + this.gaussianAlpha;
2933 g.drawString(filter, xBot, yBot + 140);
2934
2935 // draw contour key
2936 yy = yBot + 100;
2937 xx = xTop + 40;
2938 double ss = Fmath.truncate(bandZero, 3);
2939 for (int i = 0; i < numBands; i++) {
2940 double ff = Fmath.truncate(band[i], 3);
2941 g.setColor(color[numBands]);
2942 g.drawString(ss + " - " + ff, xx + 25, yy);
2943 ss = ff;
2944 g.setColor(color[i]);
2945 for (int j = 0; j < 20; j++) {
2946 yy = yy - 1;
2947 g.drawLine(xx, yy, xx + 20, yy);
2948 }
2949 }
2950 g.setColor(Color.black);
2951 g.drawString("Mean square", xx + 25, yy - 25);
2952 g.drawString("amplitudes ", xx + 25, yy - 10);
2953
2954 }
2955
2956 // returns nearest power of two that is equal to or lower than argument
2957 // length
2958 public static int lastPowerOfTwo(int len) {
2959
2960 boolean test0 = true;
2961 while (test0) {
2962 if (FourierTransform.checkPowerOfTwo(len)) {
2963 test0 = false;
2964 } else {
2965 len--;
2966 }
2967 }
2968 return len;
2969 }
2970
2971 // returns nearest power of two that is equal to or higher than argument
2972 // length
2973 public static int nextPowerOfTwo(int len) {
2974
2975 boolean test0 = true;
2976 while (test0) {
2977 if (FourierTransform.checkPowerOfTwo(len)) {
2978 test0 = false;
2979 } else {
2980 len++;
2981 }
2982 }
2983 return len;
2984 }
2985
2986 // Checks whether the argument n is a power of 2
2987 public static boolean checkPowerOfTwo(int n) {
2988 boolean test = true;
2989 int m = n;
2990 while (test && m > 1) {
2991 if ((m % 2) != 0) {
2992 test = false;
2993 } else {
2994 m /= 2;
2995 }
2996 }
2997 return test;
2998 }
2999
3000 // Checks whether the argument n is an integer times a integer power of 2
3001 // returns integer multiplier if true
3002 // returns zero if false
3003 public static int checkIntegerTimesPowerOfTwo(int n) {
3004 boolean testOuter1 = true;
3005 boolean testInner1 = true;
3006 boolean testInner2 = true;
3007 boolean testReturn = true;
3008
3009 int m = n;
3010 int j = 1;
3011 int mult = 0;
3012
3013 while (testOuter1) {
3014 testInner1 = FourierTransform.checkPowerOfTwo(m);
3015 if (testInner1) {
3016 testReturn = true;
3017 testOuter1 = false;
3018 } else {
3019 testInner2 = true;
3020 while (testInner2) {
3021 m /= ++j;
3022 if (m < 1) {
3023 testInner2 = false;
3024 testInner1 = false;
3025 testOuter1 = false;
3026 testReturn = false;
3027 } else {
3028 if ((m % 2) == 0)
3029 testInner2 = false;
3030 }
3031 }
3032 }
3033 }
3034 if (testReturn)
3035 mult = j;
3036 return mult;
3037 }
3038
3039 // Return the serial version unique identifier
3040 public static long getSerialVersionUID() {
3041 return FourierTransform.serialVersionUID;
3042 }
3043
3044}
Note: See TracBrowser for help on using the repository browser.