001/* 002 * This file is part of McIDAS-V 003 * 004 * Copyright 2007-2024 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 https://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 Complex input and output data in 051 * arrays x and y Length is n. Inputs : RealArray_Input - Input data array 052 * to perform FFT analysis CmplxArray_Input - Empty on input NumBins - 053 * Number of histogram bins in input array Outputs : RealArray_Input - Real 054 * part of FFT Analysis CmplxArray_Input - Complex part of FFT Analysis 055 * 056 * @return Values {@code <= 0} are errors, while anything {@code > 0} is ok. 057 */ 058 private static int dfft() { 059 /* real values array */ 060 double[] RealArr = new double[FFTBINS + 1]; 061 062 /* complex values array */ 063 double[] CmplxArr = new double[FFTBINS + 1]; 064 065 int LocalA; 066 int LocalB; 067 int LocalC; 068 int LocalD; 069 int LocalE; 070 int LocalX0; 071 int LocalX1; 072 int LocalX2; 073 int LocalX3; 074 int LocalY; 075 int LocalZ; 076 int LocalE1; 077 int LocalE2; 078 int LocalE4; 079 double LocalDblA; 080 double LocalDblB; 081 double LocalDblA3; 082 double LocalDblX1; 083 double LocalDblX3; 084 double LocalDblY1; 085 double LocalDblY3; 086 double LocalDblW1; 087 double LocalDblW2; 088 double LocalDblZ1; 089 double LocalDblZ2; 090 double LocalDblZ3; 091 double LocalDblXt; 092 int i; 093 094 RealArr[0] = 0.0; 095 CmplxArr[0] = 0.0; 096 for (i = 1; i <= FFTBINS; i++) { 097 RealArr[i] = fftReal[i - 1]; 098 CmplxArr[i] = fftComplex[i - 1]; 099 } 100 LocalA = 2; 101 LocalD = 1; 102 while (LocalA < FFTBINS) { 103 LocalA = LocalA + LocalA; 104 LocalD = LocalD + 1; 105 } 106 LocalE = LocalA; 107 108 if (LocalE != FFTBINS) { 109 for (LocalA = FFTBINS + 1; LocalA <= LocalE; LocalA++) { 110 RealArr[LocalA] = 0.0; 111 CmplxArr[LocalA] = 0.0; 112 } 113 } 114 115 LocalE2 = LocalE + LocalE; 116 for (LocalC = 1; LocalC <= LocalD - 1; LocalC++) { 117 LocalE2 = LocalE2 / 2; 118 LocalE4 = LocalE2 / 4; 119 LocalDblB = 2.0 * PI / LocalE2; 120 LocalDblA = 0.0; 121 for (LocalB = 1; LocalB <= LocalE4; LocalB++) { 122 LocalDblA3 = 3.0 * LocalDblA; 123 LocalDblX1 = cos(LocalDblA); 124 LocalDblY1 = sin(LocalDblA); 125 LocalDblX3 = cos(LocalDblA3); 126 LocalDblY3 = sin(LocalDblA3); 127 LocalDblA = ((double) LocalB) * LocalDblB; 128 LocalY = LocalB; 129 LocalZ = 2 * LocalE2; 130 while (LocalY < LocalE) { 131 132 for (LocalX0 = LocalY; LocalX0 <= LocalE - 1; LocalX0 = LocalX0 + LocalZ) { 133 LocalX1 = LocalX0 + LocalE4; 134 LocalX2 = LocalX1 + LocalE4; 135 LocalX3 = LocalX2 + LocalE4; 136 LocalDblW1 = RealArr[LocalX0] - RealArr[LocalX2]; 137 138 RealArr[LocalX0] = RealArr[LocalX0] + RealArr[LocalX2]; 139 LocalDblW2 = RealArr[LocalX1] - RealArr[LocalX3]; 140 RealArr[LocalX1] = RealArr[LocalX1] + RealArr[LocalX3]; 141 LocalDblZ1 = CmplxArr[LocalX0] - CmplxArr[LocalX2]; 142 CmplxArr[LocalX0] = CmplxArr[LocalX0] + CmplxArr[LocalX2]; 143 LocalDblZ2 = CmplxArr[LocalX1] - CmplxArr[LocalX3]; 144 CmplxArr[LocalX1] = CmplxArr[LocalX1] + CmplxArr[LocalX3]; 145 LocalDblZ3 = LocalDblW1 - LocalDblZ2; 146 LocalDblW1 = LocalDblW1 + LocalDblZ2; 147 LocalDblZ2 = LocalDblW2 - LocalDblZ1; 148 LocalDblW2 = LocalDblW2 + LocalDblZ1; 149 RealArr[LocalX2] = LocalDblW1 * LocalDblX1 - LocalDblZ2 * LocalDblY1; 150 CmplxArr[LocalX2] = -LocalDblZ2 * LocalDblX1 - LocalDblW1 * LocalDblY1; 151 RealArr[LocalX3] = LocalDblZ3 * LocalDblX3 + LocalDblW2 * LocalDblY3; 152 CmplxArr[LocalX3] = LocalDblW2 * LocalDblX3 - LocalDblZ3 * LocalDblY3; 153 } 154 LocalY = 2 * LocalZ - LocalE2 + LocalB; 155 LocalZ = 4 * LocalZ; 156 } 157 } 158 } 159 160 /* 161 * ---------------------Last stage, length=2 162 * butterfly--------------------- 163 */ 164 LocalY = 1; 165 LocalZ = 4; 166 while (LocalY < LocalE) { 167 for (LocalX0 = LocalY; LocalX0 <= LocalE; LocalX0 = LocalX0 + LocalZ) { 168 LocalX1 = LocalX0 + 1; 169 LocalDblW1 = RealArr[LocalX0]; 170 RealArr[LocalX0] = LocalDblW1 + RealArr[LocalX1]; 171 RealArr[LocalX1] = LocalDblW1 - RealArr[LocalX1]; 172 LocalDblW1 = CmplxArr[LocalX0]; 173 CmplxArr[LocalX0] = LocalDblW1 + CmplxArr[LocalX1]; 174 CmplxArr[LocalX1] = LocalDblW1 - CmplxArr[LocalX1]; 175 } 176 LocalY = 2 * LocalZ - 1; 177 LocalZ = 4 * LocalZ; 178 } 179 180 /* 181 * c--------------------------Bit reverse counter 182 */ 183 LocalB = 1; 184 LocalE1 = LocalE - 1; 185 for (LocalA = 1; LocalA <= LocalE1; LocalA++) { 186 if (LocalA < LocalB) { 187 LocalDblXt = RealArr[LocalB]; 188 RealArr[LocalB] = RealArr[LocalA]; 189 RealArr[LocalA] = LocalDblXt; 190 LocalDblXt = CmplxArr[LocalB]; 191 CmplxArr[LocalB] = CmplxArr[LocalA]; 192 CmplxArr[LocalA] = LocalDblXt; 193 } 194 LocalC = LocalE / 2; 195 while (LocalC < LocalB) { 196 LocalB = LocalB - LocalC; 197 LocalC = LocalC / 2; 198 } 199 LocalB = LocalB + LocalC; 200 } 201 202 /* write Real/CmplxArr back to FFT_Real/Comples arrays */ 203 for (i = 1; i <= FFTBINS; i++) { 204 fftReal[i - 1] = RealArr[i]; 205 fftComplex[i - 1] = CmplxArr[i]; 206 } 207 208 return LocalE; 209 210 } 211 212 private static double complexAbs(double realValue, double imaginaryValue) { 213 double storageValue; 214 215 if (realValue < 0.0) { 216 realValue = -realValue; 217 } 218 if (imaginaryValue < 0.0) { 219 imaginaryValue = -imaginaryValue; 220 } 221 if (imaginaryValue > realValue) { 222 storageValue = realValue; 223 realValue = imaginaryValue; 224 imaginaryValue = storageValue; 225 } 226 227 final double complexAbs; 228 if ((realValue + imaginaryValue) == realValue) { 229 complexAbs = realValue; 230 } else { 231 storageValue = imaginaryValue / realValue; 232 storageValue = realValue * sqrt(1.0 + (storageValue * storageValue)); 233 complexAbs = storageValue; 234 } 235 return complexAbs; 236 } 237 238 public static int calculateFFT(double[] inputArray) { 239 240 int fftValue = -99; 241 242 for (int i = 0; i < FFTBINS; i++) { 243 fftReal[i] = inputArray[i]; 244 fftComplex[i] = 0.0; 245 fftMagnitude[i] = 0.0; 246 } 247 248 int retErr = FFT.dfft(); 249 if (retErr <= 0) { 250 /* throw exception */ 251 } else { 252 int harmonicCounter = 0; 253 254 for (int i = 0; i < FFTBINS; i++) { 255 fftMagnitude[i] = FFT.complexAbs(fftReal[i], fftComplex[i]); 256 /* 257 * System.out.printf( 258 * "arrayinc=%d FFT real=%f cmplx=%f magnitude=%f\n" 259 * ,i,fftReal[i],fftComplex[i],fftMagnitude[i]); 260 */ 261 } 262 263 double fftTotalAllBins = 0.0; 264 for (int i = 2; i <= 31; i++) { 265 double fftBinM2 = fftMagnitude[i - 2]; 266 double fftBinM1 = fftMagnitude[i - 1]; 267 fftTotalAllBins = fftTotalAllBins + (fftBinM1 + fftBinM2) / 2.0; 268 if ((fftMagnitude[i - 1] > fftMagnitude[i - 2]) 269 && (fftMagnitude[i - 1] > fftMagnitude[i])) { 270 ++harmonicCounter; 271 /* 272 * System.out.printf("i=%d magnitude=%f counter=%d\n",i, 273 * fftMagnitude[i],harmonicCounter); 274 */ 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}