001/* ***** BEGIN LICENSE BLOCK ***** 002 * Version: MPL 1.1/GPL 2.0/LGPL 2.1 003 * 004 * The contents of this file are subject to the Mozilla Public License Version 005 * 1.1 (the "License"); you may not use this file except in compliance with 006 * the License. You may obtain a copy of the License at 007 * http://www.mozilla.org/MPL/ 008 * 009 * Software distributed under the License is distributed on an "AS IS" basis, 010 * WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License 011 * for the specific language governing rights and limitations under the 012 * License. 013 * 014 * The Original Code is JTransforms. 015 * 016 * The Initial Developer of the Original Code is 017 * Piotr Wendykier, Emory University. 018 * Portions created by the Initial Developer are Copyright (C) 2007-2009 019 * the Initial Developer. All Rights Reserved. 020 * 021 * Alternatively, the contents of this file may be used under the terms of 022 * either the GNU General Public License Version 2 or later (the "GPL"), or 023 * the GNU Lesser General Public License Version 2.1 or later (the "LGPL"), 024 * in which case the provisions of the GPL or the LGPL are applicable instead 025 * of those above. If you wish to allow use of your version of this file only 026 * under the terms of either the GPL or the LGPL, and not to allow others to 027 * use your version of this file under the terms of the MPL, indicate your 028 * decision by deleting the provisions above and replace them with the notice 029 * and other provisions required by the GPL or the LGPL. If you do not delete 030 * the provisions above, a recipient may use your version of this file under 031 * the terms of any one of the MPL, the GPL or the LGPL. 032 * 033 * ***** END LICENSE BLOCK ***** */ 034 035package edu.emory.mathcs.jtransforms.fft; 036 037import edu.emory.mathcs.utils.IOUtils; 038 039/** 040 * Accuracy check of double precision FFT's 041 * 042 * @author Piotr Wendykier (piotr.wendykier@gmail.com) 043 * 044 */ 045public class AccuracyCheckDoubleFFT { 046 047 private static int[] sizes1D = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 16, 32, 64, 100, 120, 128, 256, 310, 512, 1024, 1056, 2048, 8192, 10158, 16384, 32768, 65530, 65536, 131072 }; 048 049 private static int[] sizes2D = { 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 16, 32, 64, 100, 120, 128, 256, 310, 511, 512, 1024 }; 050 051 private static int[] sizes3D = { 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 16, 32, 64, 100, 128 }; 052 053 private static int[] sizes2D2 = { 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024 }; 054 055 private static int[] sizes3D2 = { 2, 4, 8, 16, 32, 64, 128 }; 056 057 private static double eps = Math.pow(2, -52); 058 059 private AccuracyCheckDoubleFFT() { 060 061 } 062 063 public static void checkAccuracyComplexFFT_1D() { 064 System.out.println("Checking accuracy of 1D complex FFT..."); 065 for (int i = 0; i < sizes1D.length; i++) { 066 DoubleFFT_1D fft = new DoubleFFT_1D(sizes1D[i]); 067 double err = 0.0; 068 double[] a = new double[2 * sizes1D[i]]; 069 IOUtils.fillMatrix_1D(2 * sizes1D[i], a); 070 double[] b = new double[2 * sizes1D[i]]; 071 IOUtils.fillMatrix_1D(2 * sizes1D[i], b); 072 fft.complexForward(a); 073 fft.complexInverse(a, true); 074 err = computeRMSE(a, b); 075 if (err > eps) { 076 System.err.println("\tsize = " + sizes1D[i] + ";\terror = " + err); 077 } else { 078 System.out.println("\tsize = " + sizes1D[i] + ";\terror = " + err); 079 } 080 a = null; 081 b = null; 082 fft = null; 083 System.gc(); 084 } 085 086 } 087 088 public static void checkAccuracyComplexFFT_2D() { 089 System.out.println("Checking accuracy of 2D complex FFT (double[] input)..."); 090 for (int i = 0; i < sizes2D.length; i++) { 091 DoubleFFT_2D fft2 = new DoubleFFT_2D(sizes2D[i], sizes2D[i]); 092 double err = 0.0; 093 double[] a = new double[2 * sizes2D[i] * sizes2D[i]]; 094 IOUtils.fillMatrix_2D(sizes2D[i], 2 * sizes2D[i], a); 095 double[] b = new double[2 * sizes2D[i] * sizes2D[i]]; 096 IOUtils.fillMatrix_2D(sizes2D[i], 2 * sizes2D[i], b); 097 fft2.complexForward(a); 098 fft2.complexInverse(a, true); 099 err = computeRMSE(a, b); 100 if (err > eps) { 101 System.err.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err); 102 } else { 103 System.out.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err); 104 } 105 a = null; 106 b = null; 107 fft2 = null; 108 System.gc(); 109 } 110 System.out.println("Checking accuracy of 2D complex FFT (double[][] input)..."); 111 for (int i = 0; i < sizes2D.length; i++) { 112 DoubleFFT_2D fft2 = new DoubleFFT_2D(sizes2D[i], sizes2D[i]); 113 double err = 0.0; 114 double[][] a = new double[sizes2D[i]][2 * sizes2D[i]]; 115 IOUtils.fillMatrix_2D(sizes2D[i], 2 * sizes2D[i], a); 116 double[][] b = new double[sizes2D[i]][2 * sizes2D[i]]; 117 IOUtils.fillMatrix_2D(sizes2D[i], 2 * sizes2D[i], b); 118 fft2.complexForward(a); 119 fft2.complexInverse(a, true); 120 err = computeRMSE(a, b); 121 if (err > eps) { 122 System.err.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err); 123 } else { 124 System.out.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err); 125 } 126 a = null; 127 b = null; 128 fft2 = null; 129 System.gc(); 130 } 131 } 132 133 public static void checkAccuracyComplexFFT_3D() { 134 System.out.println("Checking accuracy of 3D complex FFT (double[] input)..."); 135 for (int i = 0; i < sizes3D.length; i++) { 136 DoubleFFT_3D fft3 = new DoubleFFT_3D(sizes3D[i], sizes3D[i], sizes3D[i]); 137 double err = 0.0; 138 double[] a = new double[2 * sizes3D[i] * sizes3D[i] * sizes3D[i]]; 139 IOUtils.fillMatrix_3D(sizes3D[i], sizes3D[i], 2 * sizes3D[i], a); 140 double[] b = new double[2 * sizes3D[i] * sizes3D[i] * sizes3D[i]]; 141 IOUtils.fillMatrix_3D(sizes3D[i], sizes3D[i], 2 * sizes3D[i], b); 142 fft3.complexForward(a); 143 fft3.complexInverse(a, true); 144 err = computeRMSE(a, b); 145 if (err > eps) { 146 System.err.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " + err); 147 } else { 148 System.out.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " + err); 149 } 150 a = null; 151 b = null; 152 fft3 = null; 153 System.gc(); 154 } 155 156 System.out.println("Checking accuracy of 3D complex FFT (double[][][] input)..."); 157 for (int i = 0; i < sizes3D.length; i++) { 158 DoubleFFT_3D fft3 = new DoubleFFT_3D(sizes3D[i], sizes3D[i], sizes3D[i]); 159 double err = 0.0; 160 double[][][] a = new double[sizes3D[i]][sizes3D[i]][2 * sizes3D[i]]; 161 IOUtils.fillMatrix_3D(sizes3D[i], sizes3D[i], 2 * sizes3D[i], a); 162 double[][][] b = new double[sizes3D[i]][sizes3D[i]][2 * sizes3D[i]]; 163 IOUtils.fillMatrix_3D(sizes3D[i], sizes3D[i], 2 * sizes3D[i], b); 164 fft3.complexForward(a); 165 fft3.complexInverse(a, true); 166 err = computeRMSE(a, b); 167 if (err > eps) { 168 System.err.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " + err); 169 } else { 170 System.out.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " + err); 171 } 172 a = null; 173 b = null; 174 fft3 = null; 175 System.gc(); 176 } 177 } 178 179 public static void checkAccuracyRealFFT_1D() { 180 System.out.println("Checking accuracy of 1D real FFT..."); 181 for (int i = 0; i < sizes1D.length; i++) { 182 DoubleFFT_1D fft = new DoubleFFT_1D(sizes1D[i]); 183 double err = 0.0; 184 double[] a = new double[sizes1D[i]]; 185 IOUtils.fillMatrix_1D(sizes1D[i], a); 186 double[] b = new double[sizes1D[i]]; 187 IOUtils.fillMatrix_1D(sizes1D[i], b); 188 fft.realForward(a); 189 fft.realInverse(a, true); 190 err = computeRMSE(a, b); 191 if (err > eps) { 192 System.err.println("\tsize = " + sizes1D[i] + ";\terror = " + err); 193 } else { 194 System.out.println("\tsize = " + sizes1D[i] + ";\terror = " + err); 195 } 196 a = null; 197 b = null; 198 fft = null; 199 System.gc(); 200 } 201 System.out.println("Checking accuracy of on 1D real forward full FFT..."); 202 for (int i = 0; i < sizes1D.length; i++) { 203 DoubleFFT_1D fft = new DoubleFFT_1D(sizes1D[i]); 204 double err = 0.0; 205 double[] a = new double[2 * sizes1D[i]]; 206 IOUtils.fillMatrix_1D(sizes1D[i], a); 207 double[] b = new double[2 * sizes1D[i]]; 208 for (int j = 0; j < sizes1D[i]; j++) { 209 b[2 * j] = a[j]; 210 } 211 fft.realForwardFull(a); 212 fft.complexInverse(a, true); 213 err = computeRMSE(a, b); 214 if (err > eps) { 215 System.err.println("\tsize = " + sizes1D[i] + ";\terror = " + err); 216 } else { 217 System.out.println("\tsize = " + sizes1D[i] + ";\terror = " + err); 218 } 219 a = null; 220 b = null; 221 fft = null; 222 System.gc(); 223 } 224 System.out.println("Checking accuracy of 1D real inverse full FFT..."); 225 for (int i = 0; i < sizes1D.length; i++) { 226 DoubleFFT_1D fft = new DoubleFFT_1D(sizes1D[i]); 227 double err = 0.0; 228 double[] a = new double[2 * sizes1D[i]]; 229 IOUtils.fillMatrix_1D(sizes1D[i], a); 230 double[] b = new double[2 * sizes1D[i]]; 231 for (int j = 0; j < sizes1D[i]; j++) { 232 b[2 * j] = a[j]; 233 } 234 fft.realInverseFull(a, true); 235 fft.complexForward(a); 236 err = computeRMSE(a, b); 237 if (err > eps) { 238 System.err.println("\tsize = " + sizes1D[i] + ";\terror = " + err); 239 } else { 240 System.out.println("\tsize = " + sizes1D[i] + ";\terror = " + err); 241 } 242 a = null; 243 b = null; 244 fft = null; 245 System.gc(); 246 } 247 248 } 249 250 public static void checkAccuracyRealFFT_2D() { 251 System.out.println("Checking accuracy of 2D real FFT (double[] input)..."); 252 for (int i = 0; i < sizes2D2.length; i++) { 253 DoubleFFT_2D fft2 = new DoubleFFT_2D(sizes2D2[i], sizes2D2[i]); 254 double err = 0.0; 255 double[] a = new double[sizes2D2[i] * sizes2D2[i]]; 256 IOUtils.fillMatrix_2D(sizes2D2[i], sizes2D2[i], a); 257 double[] b = new double[sizes2D2[i] * sizes2D2[i]]; 258 IOUtils.fillMatrix_2D(sizes2D2[i], sizes2D2[i], b); 259 fft2.realForward(a); 260 fft2.realInverse(a, true); 261 err = computeRMSE(a, b); 262 if (err > eps) { 263 System.err.println("\tsize = " + sizes2D2[i] + " x " + sizes2D2[i] + ";\terror = " + err); 264 } else { 265 System.out.println("\tsize = " + sizes2D2[i] + " x " + sizes2D2[i] + ";\terror = " + err); 266 } 267 a = null; 268 b = null; 269 fft2 = null; 270 System.gc(); 271 } 272 System.out.println("Checking accuracy of 2D real FFT (double[][] input)..."); 273 for (int i = 0; i < sizes2D2.length; i++) { 274 DoubleFFT_2D fft2 = new DoubleFFT_2D(sizes2D2[i], sizes2D2[i]); 275 double err = 0.0; 276 double[][] a = new double[sizes2D2[i]][sizes2D2[i]]; 277 IOUtils.fillMatrix_2D(sizes2D2[i], sizes2D2[i], a); 278 double[][] b = new double[sizes2D2[i]][sizes2D2[i]]; 279 IOUtils.fillMatrix_2D(sizes2D2[i], sizes2D2[i], b); 280 fft2.realForward(a); 281 fft2.realInverse(a, true); 282 err = computeRMSE(a, b); 283 if (err > eps) { 284 System.err.println("\tsize = " + sizes2D2[i] + " x " + sizes2D2[i] + ";\terror = " + err); 285 } else { 286 System.out.println("\tsize = " + sizes2D2[i] + " x " + sizes2D2[i] + ";\terror = " + err); 287 } 288 a = null; 289 b = null; 290 fft2 = null; 291 System.gc(); 292 } 293 System.out.println("Checking accuracy of 2D real forward full FFT (double[] input)..."); 294 for (int i = 0; i < sizes2D.length; i++) { 295 DoubleFFT_2D fft2 = new DoubleFFT_2D(sizes2D[i], sizes2D[i]); 296 double err = 0.0; 297 double[] a = new double[2 * sizes2D[i] * sizes2D[i]]; 298 IOUtils.fillMatrix_2D(sizes2D[i], sizes2D[i], a); 299 double[] b = new double[2 * sizes2D[i] * sizes2D[i]]; 300 for (int r = 0; r < sizes2D[i]; r++) { 301 for (int c = 0; c < sizes2D[i]; c++) { 302 b[r * 2 * sizes2D[i] + 2 * c] = a[r * sizes2D[i] + c]; 303 } 304 } 305 fft2.realForwardFull(a); 306 fft2.complexInverse(a, true); 307 err = computeRMSE(a, b); 308 if (err > eps) { 309 System.err.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err); 310 } else { 311 System.out.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err); 312 } 313 a = null; 314 b = null; 315 fft2 = null; 316 System.gc(); 317 } 318 System.out.println("Checking accuracy of 2D real forward full FFT (double[][] input)..."); 319 for (int i = 0; i < sizes2D.length; i++) { 320 DoubleFFT_2D fft2 = new DoubleFFT_2D(sizes2D[i], sizes2D[i]); 321 double err = 0.0; 322 double[][] a = new double[sizes2D[i]][2 * sizes2D[i]]; 323 IOUtils.fillMatrix_2D(sizes2D[i], sizes2D[i], a); 324 double[][] b = new double[sizes2D[i]][2 * sizes2D[i]]; 325 for (int r = 0; r < sizes2D[i]; r++) { 326 for (int c = 0; c < sizes2D[i]; c++) { 327 b[r][2 * c] = a[r][c]; 328 } 329 } 330 fft2.realForwardFull(a); 331 fft2.complexInverse(a, true); 332 err = computeRMSE(a, b); 333 if (err > eps) { 334 System.err.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err); 335 } else { 336 System.out.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err); 337 } 338 a = null; 339 b = null; 340 fft2 = null; 341 System.gc(); 342 } 343 System.out.println("Checking accuracy of 2D real inverse full FFT (double[] input)..."); 344 for (int i = 0; i < sizes2D.length; i++) { 345 DoubleFFT_2D fft2 = new DoubleFFT_2D(sizes2D[i], sizes2D[i]); 346 double err = 0.0; 347 double[] a = new double[2 * sizes2D[i] * sizes2D[i]]; 348 IOUtils.fillMatrix_2D(sizes2D[i], sizes2D[i], a); 349 double[] b = new double[2 * sizes2D[i] * sizes2D[i]]; 350 for (int r = 0; r < sizes2D[i]; r++) { 351 for (int c = 0; c < sizes2D[i]; c++) { 352 b[r * 2 * sizes2D[i] + 2 * c] = a[r * sizes2D[i] + c]; 353 } 354 } 355 fft2.realInverseFull(a, true); 356 fft2.complexForward(a); 357 err = computeRMSE(a, b); 358 if (err > eps) { 359 System.err.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err); 360 } else { 361 System.out.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err); 362 } 363 a = null; 364 b = null; 365 fft2 = null; 366 System.gc(); 367 } 368 System.out.println("Checking accuracy of 2D real inverse full FFT (double[][] input)..."); 369 for (int i = 0; i < sizes2D.length; i++) { 370 DoubleFFT_2D fft2 = new DoubleFFT_2D(sizes2D[i], sizes2D[i]); 371 double err = 0.0; 372 double[][] a = new double[sizes2D[i]][2 * sizes2D[i]]; 373 IOUtils.fillMatrix_2D(sizes2D[i], sizes2D[i], a); 374 double[][] b = new double[sizes2D[i]][2 * sizes2D[i]]; 375 for (int r = 0; r < sizes2D[i]; r++) { 376 for (int c = 0; c < sizes2D[i]; c++) { 377 b[r][2 * c] = a[r][c]; 378 } 379 } 380 fft2.realInverseFull(a, true); 381 fft2.complexForward(a); 382 err = computeRMSE(a, b); 383 if (err > eps) { 384 System.err.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err); 385 } else { 386 System.out.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err); 387 } 388 a = null; 389 b = null; 390 fft2 = null; 391 System.gc(); 392 } 393 394 } 395 396 public static void checkAccuracyRealFFT_3D() { 397 System.out.println("Checking accuracy of 3D real FFT (double[] input)..."); 398 for (int i = 0; i < sizes3D2.length; i++) { 399 DoubleFFT_3D fft3 = new DoubleFFT_3D(sizes3D2[i], sizes3D2[i], sizes3D2[i]); 400 double err = 0.0; 401 double[] a = new double[sizes3D2[i] * sizes3D2[i] * sizes3D2[i]]; 402 IOUtils.fillMatrix_3D(sizes3D2[i], sizes3D2[i], sizes3D2[i], a); 403 double[] b = new double[sizes3D2[i] * sizes3D2[i] * sizes3D2[i]]; 404 IOUtils.fillMatrix_3D(sizes3D2[i], sizes3D2[i], sizes3D2[i], b); 405 fft3.realForward(b); 406 fft3.realInverse(b, true); 407 err = computeRMSE(a, b); 408 if (err > eps) { 409 System.err.println("\tsize = " + sizes3D2[i] + " x " + sizes3D2[i] + " x " + sizes3D2[i] + ";\t\terror = " + err); 410 } else { 411 System.out.println("\tsize = " + sizes3D2[i] + " x " + sizes3D2[i] + " x " + sizes3D2[i] + ";\t\terror = " + err); 412 } 413 a = null; 414 b = null; 415 fft3 = null; 416 System.gc(); 417 } 418 System.out.println("Checking accuracy of 3D real FFT (double[][][] input)..."); 419 for (int i = 0; i < sizes3D2.length; i++) { 420 DoubleFFT_3D fft3 = new DoubleFFT_3D(sizes3D2[i], sizes3D2[i], sizes3D2[i]); 421 double err = 0.0; 422 double[][][] a = new double[sizes3D2[i]][sizes3D2[i]][sizes3D2[i]]; 423 IOUtils.fillMatrix_3D(sizes3D2[i], sizes3D2[i], sizes3D2[i], a); 424 double[][][] b = new double[sizes3D2[i]][sizes3D2[i]][sizes3D2[i]]; 425 IOUtils.fillMatrix_3D(sizes3D2[i], sizes3D2[i], sizes3D2[i], b); 426 fft3.realForward(b); 427 fft3.realInverse(b, true); 428 err = computeRMSE(a, b); 429 if (err > eps) { 430 System.err.println("\tsize = " + sizes3D2[i] + " x " + sizes3D2[i] + " x " + sizes3D2[i] + ";\t\terror = " + err); 431 } else { 432 System.out.println("\tsize = " + sizes3D2[i] + " x " + sizes3D2[i] + " x " + sizes3D2[i] + ";\t\terror = " + err); 433 } 434 a = null; 435 b = null; 436 fft3 = null; 437 System.gc(); 438 } 439 System.out.println("Checking accuracy of 3D real forward full FFT (double[] input)..."); 440 for (int i = 0; i < sizes3D.length; i++) { 441 DoubleFFT_3D fft3 = new DoubleFFT_3D(sizes3D[i], sizes3D[i], sizes3D[i]); 442 double err = 0.0; 443 double[] a = new double[2 * sizes3D[i] * sizes3D[i] * sizes3D[i]]; 444 IOUtils.fillMatrix_3D(sizes3D[i], sizes3D[i], sizes3D[i], a); 445 double[] b = new double[2 * sizes3D[i] * sizes3D[i] * sizes3D[i]]; 446 for (int s = 0; s < sizes3D[i]; s++) { 447 for (int r = 0; r < sizes3D[i]; r++) { 448 for (int c = 0; c < sizes3D[i]; c++) { 449 b[s * 2 * sizes3D[i] * sizes3D[i] + r * 2 * sizes3D[i] + 2 * c] = a[s * sizes3D[i] * sizes3D[i] + r * sizes3D[i] + c]; 450 } 451 } 452 } 453 fft3.realForwardFull(a); 454 fft3.complexInverse(a, true); 455 err = computeRMSE(a, b); 456 if (err > eps) { 457 System.err.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " + err); 458 } else { 459 System.out.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " + err); 460 } 461 a = null; 462 b = null; 463 fft3 = null; 464 System.gc(); 465 } 466 467 System.out.println("Checking accuracy of 3D real forward full FFT (double[][][] input)..."); 468 for (int i = 0; i < sizes3D.length; i++) { 469 DoubleFFT_3D fft3 = new DoubleFFT_3D(sizes3D[i], sizes3D[i], sizes3D[i]); 470 double err = 0.0; 471 double[][][] a = new double[sizes3D[i]][sizes3D[i]][2 * sizes3D[i]]; 472 IOUtils.fillMatrix_3D(sizes3D[i], sizes3D[i], sizes3D[i], a); 473 double[][][] b = new double[sizes3D[i]][sizes3D[i]][2 * sizes3D[i]]; 474 for (int s = 0; s < sizes3D[i]; s++) { 475 for (int r = 0; r < sizes3D[i]; r++) { 476 for (int c = 0; c < sizes3D[i]; c++) { 477 b[s][r][2 * c] = a[s][r][c]; 478 } 479 } 480 } 481 fft3.realForwardFull(a); 482 fft3.complexInverse(a, true); 483 err = computeRMSE(a, b); 484 if (err > eps) { 485 System.err.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " + err); 486 } else { 487 System.out.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " + err); 488 } 489 a = null; 490 b = null; 491 fft3 = null; 492 System.gc(); 493 } 494 495 System.out.println("Checking accuracy of 3D real inverse full FFT (double[] input)..."); 496 for (int i = 0; i < sizes3D.length; i++) { 497 DoubleFFT_3D fft3 = new DoubleFFT_3D(sizes3D[i], sizes3D[i], sizes3D[i]); 498 double err = 0.0; 499 double[] a = new double[2 * sizes3D[i] * sizes3D[i] * sizes3D[i]]; 500 IOUtils.fillMatrix_3D(sizes3D[i], sizes3D[i], sizes3D[i], a); 501 double[] b = new double[2 * sizes3D[i] * sizes3D[i] * sizes3D[i]]; 502 for (int s = 0; s < sizes3D[i]; s++) { 503 for (int r = 0; r < sizes3D[i]; r++) { 504 for (int c = 0; c < sizes3D[i]; c++) { 505 b[s * 2 * sizes3D[i] * sizes3D[i] + r * 2 * sizes3D[i] + 2 * c] = a[s * sizes3D[i] * sizes3D[i] + r * sizes3D[i] + c]; 506 } 507 } 508 } 509 fft3.realInverseFull(a, true); 510 fft3.complexForward(a); 511 err = computeRMSE(a, b); 512 if (err > eps) { 513 System.err.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " + err); 514 } else { 515 System.out.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " + err); 516 } 517 a = null; 518 b = null; 519 fft3 = null; 520 System.gc(); 521 } 522 System.out.println("Checking accuracy of 3D real inverse full FFT (double[][][] input)..."); 523 for (int i = 0; i < sizes3D.length; i++) { 524 DoubleFFT_3D fft3 = new DoubleFFT_3D(sizes3D[i], sizes3D[i], sizes3D[i]); 525 double err = 0.0; 526 double[][][] a = new double[sizes3D[i]][sizes3D[i]][2 * sizes3D[i]]; 527 IOUtils.fillMatrix_3D(sizes3D[i], sizes3D[i], sizes3D[i], a); 528 double[][][] b = new double[sizes3D[i]][sizes3D[i]][2 * sizes3D[i]]; 529 for (int s = 0; s < sizes3D[i]; s++) { 530 for (int r = 0; r < sizes3D[i]; r++) { 531 for (int c = 0; c < sizes3D[i]; c++) { 532 b[s][r][2 * c] = a[s][r][c]; 533 } 534 } 535 } 536 fft3.realInverseFull(a, true); 537 fft3.complexForward(a); 538 err = computeRMSE(a, b); 539 if (err > eps) { 540 System.err.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " + err); 541 } else { 542 System.out.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " + err); 543 } 544 a = null; 545 b = null; 546 fft3 = null; 547 System.gc(); 548 } 549 } 550 551 private static double computeRMSE(double[] a, double[] b) { 552 if (a.length != b.length) { 553 throw new IllegalArgumentException("Arrays are not the same size."); 554 } 555 double rms = 0; 556 double tmp; 557 for (int i = 0; i < a.length; i++) { 558 tmp = (a[i] - b[i]); 559 rms += tmp * tmp; 560 } 561 return Math.sqrt(rms / (double) a.length); 562 } 563 564 private static double computeRMSE(double[][] a, double[][] b) { 565 if (a.length != b.length || a[0].length != b[0].length) { 566 throw new IllegalArgumentException("Arrays are not the same size."); 567 } 568 double rms = 0; 569 double tmp; 570 for (int r = 0; r < a.length; r++) { 571 for (int c = 0; c < a[0].length; c++) { 572 tmp = (a[r][c] - b[r][c]); 573 rms += tmp * tmp; 574 } 575 } 576 return Math.sqrt(rms / (a.length * a[0].length)); 577 } 578 579 private static double computeRMSE(double[][][] a, double[][][] b) { 580 if (a.length != b.length || a[0].length != b[0].length || a[0][0].length != b[0][0].length) { 581 throw new IllegalArgumentException("Arrays are not the same size."); 582 } 583 double rms = 0; 584 double tmp; 585 for (int s = 0; s < a.length; s++) { 586 for (int r = 0; r < a[0].length; r++) { 587 for (int c = 0; c < a[0][0].length; c++) { 588 tmp = (a[s][r][c] - b[s][r][c]); 589 rms += tmp * tmp; 590 } 591 } 592 } 593 return Math.sqrt(rms / (a.length * a[0].length * a[0][0].length)); 594 } 595 596 public static void main(String[] args) { 597 checkAccuracyComplexFFT_1D(); 598 checkAccuracyRealFFT_1D(); 599 checkAccuracyComplexFFT_2D(); 600 checkAccuracyRealFFT_2D(); 601 checkAccuracyComplexFFT_3D(); 602 checkAccuracyRealFFT_3D(); 603 System.exit(0); 604 } 605}