001/* 002 * This file is part of McIDAS-V 003 * 004 * Copyright 2007-2017 005 * Space Science and Engineering Center (SSEC) 006 * University of Wisconsin - Madison 007 * 1225 W. Dayton Street, Madison, WI 53706, USA 008 * https://www.ssec.wisc.edu/mcidas 009 * 010 * All Rights Reserved 011 * 012 * McIDAS-V is built on Unidata's IDV and SSEC's VisAD libraries, and 013 * some McIDAS-V source code is based on IDV and VisAD source code. 014 * 015 * McIDAS-V is free software; you can redistribute it and/or modify 016 * it under the terms of the GNU Lesser Public License as published by 017 * the Free Software Foundation; either version 3 of the License, or 018 * (at your option) any later version. 019 * 020 * McIDAS-V is distributed in the hope that it will be useful, 021 * but WITHOUT ANY WARRANTY; without even the implied warranty of 022 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 023 * GNU Lesser Public License for more details. 024 * 025 * You should have received a copy of the GNU Lesser Public License 026 * along with this program. If not, see http://www.gnu.org/licenses. 027 */ 028 029package edu.wisc.ssec.mcidasv.adt; 030 031import static java.lang.Math.PI; 032import static java.lang.Math.cos; 033import static java.lang.Math.sin; 034import static java.lang.Math.sqrt; 035 036public class FFT { 037 038 private static int FFTBINS = 64; 039 040 private static double fftReal[] = new double[FFTBINS]; 041 private static double fftComplex[] = new double[FFTBINS]; 042 private static double fftMagnitude[] = new double[FFTBINS]; 043 044 private FFT() { 045 } 046 047 /** 048 * A Duhamel-Hollman split-radix dif FFT. 049 * 050 * Ref: Electronics Letters, Jan. 5, 1984 051 * Complex input and output data in arrays x and y 052 * Length is n. 053 * Inputs : RealArray_Input - Input data array to perform FFT analysis 054 * CmplxArray_Input - Empty on input 055 * NumBins - Number of histogram bins in input array 056 * Outputs : RealArray_Input - Real part of FFT Analysis 057 * CmplxArray_Input - Complex part of FFT Analysis 058 * 059 * @return Values {@code <= 0} are errors, while anything {@code > 0} is ok. 060 */ 061 private static int dfft() { 062 /* real values array */ 063 double[] RealArr = new double[FFTBINS + 1]; 064 065 /* complex values array */ 066 double[] CmplxArr = new double[FFTBINS + 1]; 067 068 int LocalA; 069 int LocalB; 070 int LocalC; 071 int LocalD; 072 int LocalE; 073 int LocalX0; 074 int LocalX1; 075 int LocalX2; 076 int LocalX3; 077 int LocalY; 078 int LocalZ; 079 int LocalE1; 080 int LocalE2; 081 int LocalE4; 082 double LocalDblA; 083 double LocalDblB; 084 double LocalDblA3; 085 double LocalDblX1; 086 double LocalDblX3; 087 double LocalDblY1; 088 double LocalDblY3; 089 double LocalDblW1; 090 double LocalDblW2; 091 double LocalDblZ1; 092 double LocalDblZ2; 093 double LocalDblZ3; 094 double LocalDblXt; 095 int i; 096 097 RealArr[0] = 0.0; 098 CmplxArr[0] = 0.0; 099 for (i = 1; i <= FFTBINS; i++ ) { 100 RealArr[i] = fftReal[i-1]; 101 CmplxArr[i] = fftComplex[i-1]; 102 } 103 LocalA = 2; 104 LocalD = 1; 105 while (LocalA < FFTBINS) { 106 LocalA = LocalA+LocalA; 107 LocalD = LocalD+1; 108 } 109 LocalE = LocalA; 110 111 if (LocalE != FFTBINS) { 112 for (LocalA = FFTBINS + 1; LocalA <= LocalE; LocalA++) { 113 RealArr[LocalA] = 0.0; 114 CmplxArr[LocalA] = 0.0; 115 } 116 } 117 118 LocalE2 = LocalE+LocalE; 119 for (LocalC = 1; LocalC <= LocalD-1; LocalC++ ) { 120 LocalE2 = LocalE2 / 2; 121 LocalE4 = LocalE2 / 4; 122 LocalDblB = 2.0 * PI / LocalE2; 123 LocalDblA = 0.0; 124 for (LocalB = 1; LocalB<= LocalE4 ; LocalB++) { 125 LocalDblA3 = 3.0*LocalDblA; 126 LocalDblX1 = cos(LocalDblA); 127 LocalDblY1 = sin(LocalDblA); 128 LocalDblX3 = cos(LocalDblA3); 129 LocalDblY3 = sin(LocalDblA3); 130 LocalDblA = ((double)LocalB)*LocalDblB; 131 LocalY = LocalB; 132 LocalZ = 2*LocalE2; 133 while ( LocalY < LocalE ) { 134 135 for (LocalX0 = LocalY; LocalX0 <= LocalE-1; 136 LocalX0 = LocalX0 + LocalZ) { 137 LocalX1 = LocalX0 + LocalE4; 138 LocalX2 = LocalX1 + LocalE4; 139 LocalX3 = LocalX2 + LocalE4; 140 LocalDblW1 = RealArr[LocalX0] - RealArr[LocalX2]; 141 142 RealArr[LocalX0] = RealArr[LocalX0] + RealArr[LocalX2]; 143 LocalDblW2 = RealArr[LocalX1] - RealArr[LocalX3]; 144 RealArr[LocalX1] = RealArr[LocalX1] + RealArr[LocalX3]; 145 LocalDblZ1 = CmplxArr[LocalX0] - CmplxArr[LocalX2]; 146 CmplxArr[LocalX0] = CmplxArr[LocalX0] + CmplxArr[LocalX2]; 147 LocalDblZ2 = CmplxArr[LocalX1] - CmplxArr[LocalX3]; 148 CmplxArr[LocalX1] = CmplxArr[LocalX1] + CmplxArr[LocalX3]; 149 LocalDblZ3 = LocalDblW1 - LocalDblZ2; 150 LocalDblW1 = LocalDblW1 + LocalDblZ2; 151 LocalDblZ2 = LocalDblW2 - LocalDblZ1; 152 LocalDblW2 = LocalDblW2 + LocalDblZ1; 153 RealArr[LocalX2] = LocalDblW1*LocalDblX1 - 154 LocalDblZ2*LocalDblY1; 155 CmplxArr[LocalX2] = -LocalDblZ2*LocalDblX1 - 156 LocalDblW1*LocalDblY1; 157 RealArr[LocalX3] = LocalDblZ3*LocalDblX3 + 158 LocalDblW2*LocalDblY3; 159 CmplxArr[LocalX3] = LocalDblW2*LocalDblX3 - 160 LocalDblZ3*LocalDblY3; 161 } 162 LocalY = 2*LocalZ - LocalE2 + LocalB; 163 LocalZ = 4*LocalZ; 164 } 165 } 166 } 167 168 /* 169 ---------------------Last stage, length=2 butterfly--------------------- 170 */ 171 LocalY = 1; 172 LocalZ = 4; 173 while ( LocalY < LocalE) { 174 for (LocalX0 = LocalY; LocalX0 <= LocalE; LocalX0 = LocalX0 + LocalZ) { 175 LocalX1 = LocalX0 + 1; LocalDblW1 = RealArr[LocalX0]; 176 RealArr[LocalX0] = LocalDblW1 + RealArr[LocalX1]; 177 RealArr[LocalX1] = LocalDblW1 - RealArr[LocalX1]; 178 LocalDblW1 = CmplxArr[LocalX0]; 179 CmplxArr[LocalX0] = LocalDblW1 + CmplxArr[LocalX1]; 180 CmplxArr[LocalX1] = LocalDblW1 - CmplxArr[LocalX1]; 181 } 182 LocalY = 2 * LocalZ - 1; 183 LocalZ = 4 * LocalZ; 184 } 185 186 /* 187 c--------------------------Bit reverse counter 188 */ 189 LocalB = 1; 190 LocalE1 = LocalE - 1; 191 for (LocalA = 1; LocalA <= LocalE1; LocalA++) { 192 if (LocalA < LocalB) { 193 LocalDblXt = RealArr[LocalB]; 194 RealArr[LocalB] = RealArr[LocalA]; 195 RealArr[LocalA] = LocalDblXt; 196 LocalDblXt = CmplxArr[LocalB]; 197 CmplxArr[LocalB] = CmplxArr[LocalA]; 198 CmplxArr[LocalA] = LocalDblXt; 199 } 200 LocalC = LocalE / 2; 201 while (LocalC < LocalB) { 202 LocalB = LocalB - LocalC; 203 LocalC = LocalC / 2; 204 } 205 LocalB = LocalB + LocalC; 206 } 207 208 /* write Real/CmplxArr back to FFT_Real/Comples arrays */ 209 for (i = 1; i <= FFTBINS; i++ ) { 210 fftReal[i-1] = RealArr[i]; 211 fftComplex[i-1] = CmplxArr[i]; 212 } 213 214 return LocalE; 215 216 } 217 218 private static double complexAbs(double realValue, double imaginaryValue) { 219 double storageValue; 220 221 if (realValue < 0.0) { 222 realValue = -realValue; 223 } 224 if (imaginaryValue < 0.0) { 225 imaginaryValue = -imaginaryValue; 226 } 227 if (imaginaryValue > realValue){ 228 storageValue = realValue; 229 realValue = imaginaryValue; 230 imaginaryValue = storageValue; 231 } 232 233 final double complexAbs; 234 if ((realValue + imaginaryValue) == realValue) { 235 complexAbs = realValue; 236 } else { 237 storageValue = imaginaryValue / realValue; 238 storageValue = realValue * sqrt(1.0 + (storageValue * storageValue)); 239 complexAbs = storageValue; 240 } 241 return complexAbs; 242 } 243 244 public static int calculateFFT(double[] inputArray) { 245 246 int fftValue = -99; 247 248 for (int i = 0; i < FFTBINS; i++ ) { 249 fftReal[i] = inputArray[i]; 250 fftComplex[i] = 0.0; 251 fftMagnitude[i] = 0.0; 252 } 253 254 int retErr = FFT.dfft(); 255 if (retErr <= 0) { 256 /* throw exception */ 257 } else { 258 int harmonicCounter = 0; 259 260 for (int i = 0; i < FFTBINS; i++ ) { 261 fftMagnitude[i] = FFT.complexAbs(fftReal[i], fftComplex[i]); 262 /* System.out.printf("arrayinc=%d FFT real=%f cmplx=%f magnitude=%f\n",i,fftReal[i],fftComplex[i],fftMagnitude[i]); */ 263 } 264 265 double fftTotalAllBins = 0.0; 266 for (int i = 2; i <= 31; i++ ) { 267 double fftBinM2 = fftMagnitude[i - 2]; 268 double fftBinM1 = fftMagnitude[i - 1]; 269 fftTotalAllBins = fftTotalAllBins + (fftBinM1+fftBinM2) / 2.0; 270 if ((fftMagnitude[i - 1] > fftMagnitude[i - 2]) && 271 (fftMagnitude[i - 1] > fftMagnitude[i])) { 272 ++harmonicCounter; 273 /* System.out.printf("i=%d magnitude=%f counter=%d\n",i, 274 * fftMagnitude[i],harmonicCounter); */ 275 } 276 } 277 if (fftMagnitude[0] == 0) { 278 /* throw exception */ 279 } else { 280 fftValue = harmonicCounter; 281 } 282 } 283 284 /* System.out.printf("Amplitude=%f\n",Amplitude); */ 285 return fftValue; 286 } 287 288}