Maintenance scheduled for Thursday, September 24th at 15:00 MDT. Expected downtime <1 hour.

...
 
Commits (4)
......@@ -27,6 +27,9 @@ import org.jfree.data.xy.XYSeries;
*/
public class FFTResult {
// 5% single side taper is default for FFT calculations
final public static double DEFAULT_TAPER_WIDTH = 0.05;
final private Complex[] transform; // the FFT data
final private double[] freqs; // array of frequencies matching the fft data
......@@ -152,7 +155,8 @@ public class FFTResult {
public static FFTResult crossPower(double[] timeseries1, double[] timeSeries2, long interval,
Complex[] response1, Complex[] response2, int range, int slider) {
FFTResult psd = FFTResult.spectralCalc(timeseries1, timeSeries2, interval, range, slider);
FFTResult psd = FFTResult.spectralCalc(timeseries1, timeSeries2, interval,
range, slider, DEFAULT_TAPER_WIDTH);
return crossPower(psd.getFFT(), psd.getFreqs(), response1, response2);
}
......@@ -487,13 +491,15 @@ public class FFTResult {
* including a tapering operation on the data to remove noise at the ends of the series
* @param toFFT timeseries data to produce FFT over
* @param padding power of 2 to pad the data out to
* @param taperWidth width of cosine taper to apply to the data (per-side) as a decimal;
* use {@link #DEFAULT_TAPER_WIDTH} for a default choice of a 5% taper
* @return pair of values, first representing the frequency-space representation
* of the positive-half of the FFT and second representing the taper power correction
*/
static Pair<Complex[], Double> getSpectralWindow(double[] toFFT, int padding) {
static Pair<Complex[], Double> getSpectralWindow(double[] toFFT, int padding, double taperWidth) {
// demean and detrend work in-place on the list
NumericUtils.demeanInPlace(toFFT);
Double wss = cosineTaper(toFFT, 0.05);
Double wss = cosineTaper(toFFT, taperWidth); // defaults to 5% width per side
// presumably we only need the last value of wss
toFFT = Arrays.copyOfRange(toFFT, 0, padding);
......@@ -749,7 +755,38 @@ public class FFTResult {
int range = list1.length / 4;
int slider = range / 4;
return spectralCalc(list1, list2, interval, range, slider);
return spectralCalc(list1, list2, interval, range, slider, DEFAULT_TAPER_WIDTH);
}
/**
* Helper function to calculate power spectral density / crosspower.
* Takes in two time series data and produces the windowed FFT over each.
* The first is multiplied by the complex conjugate of the second.
* If the two series are the same, this is the PSD of that series. If they
* are different, this result is the crosspower.
* The result is smoothed but does not have the frequency response applied,
* and so does not give a full result -- this is merely a helper function
* for the crossPower function.
*
* @param list1 First list of data to be given as input
* @param list2 Second list of data to be given as input, which can be
* the same as the first (and if so, is ignored)
* @param interval Interval of the data (same for both lists)
* @return FFTResult (FFT values and frequencies as a pair of arrays)
* representing the power-spectral density / crosspower of the input data.
*/
public static FFTResult
spectralCalc(double[] list1, double[] list2, long interval, double taperWidth) {
//Only the same data if the arrays are actually the same objects.
//noinspection ArrayEquals
boolean sameData = list1.equals(list2);
// divide into windows of 1/4, moving up 1/16 of the data at a time
int range = list1.length / 4;
int slider = range / 4;
return spectralCalc(list1, list2, interval, range, slider, taperWidth);
}
......@@ -768,11 +805,12 @@ public class FFTResult {
* the same as the first (and if so, is ignored)
* @param interval Sampling interval of the data, in ms
* @param range Length of data in a single window
* @param slider Length by which
* @return
* @param slider Length by which to adjust the window over the full data
* (i.e., if range is 1/4, then a window of 1/16 has 75% overlap between PSD window calculations)
* @return Crosspower of given timeseries data
*/
public static FFTResult spectralCalc(double[] list1, double[] list2, long interval,
int range, int slider) {
int range, int slider, double taperWidth) {
// As in the preceding method, only the same data if the arrays are actually the same objects.
// noinspection ArrayEquals
......@@ -818,12 +856,12 @@ public class FFTResult {
}
// this function converts the current data to the padded size, applies taper, and gets FFT
Pair<Complex[], Double> windFFTData = getSpectralWindow(toFFT1, padding);
Pair<Complex[], Double> windFFTData = getSpectralWindow(toFFT1, padding, taperWidth);
Complex[] fftResult1 = windFFTData.getFirst(); // actual fft data
wss = windFFTData.getSecond(); // represents some measure of power loss
Complex[] fftResult2 = fftResult1;
if (toFFT2 != null) {
fftResult2 = getSpectralWindow(toFFT2, padding).getFirst();
fftResult2 = getSpectralWindow(toFFT2, padding, taperWidth).getFirst();
}
for (int i = 0; i < singleSide; ++i) {
......
package asl.utils;
import static asl.utils.FFTResult.DEFAULT_TAPER_WIDTH;
import static org.junit.Assert.assertEquals;
import static org.junit.Assert.assertNotEquals;
import static org.junit.Assert.assertTrue;
......@@ -385,7 +386,7 @@ public class FFTResultTest {
// FFTResult.cosineTaper(toFFT, 0.05);
assertEquals(range, toFFT.length);
Pair<Complex[], Double> tempResult = FFTResult
.getSpectralWindow(toFFT, padding);
.getSpectralWindow(toFFT, padding, DEFAULT_TAPER_WIDTH);
double cos = tempResult.getSecond();
assertEquals(20249, cos, 1E-10);
Complex[] result = tempResult.getFirst();
......@@ -457,7 +458,7 @@ public class FFTResultTest {
// FFTResult.cosineTaper(toFFT, 0.05);
assertEquals(range, toFFT.length);
Pair<Complex[], Double> tempResult = FFTResult
.getSpectralWindow(toFFT, padding);
.getSpectralWindow(toFFT, padding, DEFAULT_TAPER_WIDTH);
double cos = tempResult.getSecond();
assertEquals(404999., cos, 1E-10);
Complex[] result = tempResult.getFirst();
......