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 */ 045@SuppressWarnings("javadoc") 046public class AccuracyCheckDoubleFFT { 047 048 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, 049 1024, 1056, 2048, 8192, 10158, 16384, 32768, 65530, 65536, 131072 }; 050 051 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, 052 512, 1024 }; 053 054 private static int[] sizes3D = { 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 16, 32, 64, 100, 128 }; 055 056 private static int[] sizes2D2 = { 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024 }; 057 058 private static int[] sizes3D2 = { 2, 4, 8, 16, 32, 64, 128 }; 059 060 private static double eps = Math.pow(2, -52); 061 062 private AccuracyCheckDoubleFFT() { 063 064 } 065 066 public static void checkAccuracyComplexFFT_1D() { 067 System.out.println("Checking accuracy of 1D complex FFT..."); 068 for (int i = 0; i < sizes1D.length; i++) { 069 DoubleFFT_1D fft = new DoubleFFT_1D(sizes1D[i]); 070 double err = 0.0; 071 double[] a = new double[2 * sizes1D[i]]; 072 IOUtils.fillMatrix_1D(2 * sizes1D[i], a); 073 double[] b = new double[2 * sizes1D[i]]; 074 IOUtils.fillMatrix_1D(2 * sizes1D[i], b); 075 fft.complexForward(a); 076 fft.complexInverse(a, true); 077 err = computeRMSE(a, b); 078 if (err > eps) { 079 System.err.println("\tsize = " + sizes1D[i] + ";\terror = " + err); 080 } else { 081 System.out.println("\tsize = " + sizes1D[i] + ";\terror = " + err); 082 } 083 a = null; 084 b = null; 085 fft = null; 086 System.gc(); 087 } 088 089 } 090 091 public static void checkAccuracyComplexFFT_2D() { 092 System.out.println("Checking accuracy of 2D complex FFT (double[] input)..."); 093 for (int i = 0; i < sizes2D.length; i++) { 094 DoubleFFT_2D fft2 = new DoubleFFT_2D(sizes2D[i], sizes2D[i]); 095 double err = 0.0; 096 double[] a = new double[2 * sizes2D[i] * sizes2D[i]]; 097 IOUtils.fillMatrix_2D(sizes2D[i], 2 * sizes2D[i], a); 098 double[] b = new double[2 * sizes2D[i] * sizes2D[i]]; 099 IOUtils.fillMatrix_2D(sizes2D[i], 2 * sizes2D[i], b); 100 fft2.complexForward(a); 101 fft2.complexInverse(a, true); 102 err = computeRMSE(a, b); 103 if (err > eps) { 104 System.err.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err); 105 } else { 106 System.out.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err); 107 } 108 a = null; 109 b = null; 110 fft2 = null; 111 System.gc(); 112 } 113 System.out.println("Checking accuracy of 2D complex FFT (double[][] input)..."); 114 for (int i = 0; i < sizes2D.length; i++) { 115 DoubleFFT_2D fft2 = new DoubleFFT_2D(sizes2D[i], sizes2D[i]); 116 double err = 0.0; 117 double[][] a = new double[sizes2D[i]][2 * sizes2D[i]]; 118 IOUtils.fillMatrix_2D(sizes2D[i], 2 * sizes2D[i], a); 119 double[][] b = new double[sizes2D[i]][2 * sizes2D[i]]; 120 IOUtils.fillMatrix_2D(sizes2D[i], 2 * sizes2D[i], b); 121 fft2.complexForward(a); 122 fft2.complexInverse(a, true); 123 err = computeRMSE(a, b); 124 if (err > eps) { 125 System.err.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err); 126 } else { 127 System.out.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err); 128 } 129 a = null; 130 b = null; 131 fft2 = null; 132 System.gc(); 133 } 134 } 135 136 public static void checkAccuracyComplexFFT_3D() { 137 System.out.println("Checking accuracy of 3D complex FFT (double[] input)..."); 138 for (int i = 0; i < sizes3D.length; i++) { 139 DoubleFFT_3D fft3 = new DoubleFFT_3D(sizes3D[i], sizes3D[i], sizes3D[i]); 140 double err = 0.0; 141 double[] a = new double[2 * sizes3D[i] * sizes3D[i] * sizes3D[i]]; 142 IOUtils.fillMatrix_3D(sizes3D[i], sizes3D[i], 2 * sizes3D[i], a); 143 double[] b = new double[2 * sizes3D[i] * sizes3D[i] * sizes3D[i]]; 144 IOUtils.fillMatrix_3D(sizes3D[i], sizes3D[i], 2 * sizes3D[i], b); 145 fft3.complexForward(a); 146 fft3.complexInverse(a, true); 147 err = computeRMSE(a, b); 148 if (err > eps) { 149 System.err.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " 150 + err); 151 } else { 152 System.out.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " 153 + err); 154 } 155 a = null; 156 b = null; 157 fft3 = null; 158 System.gc(); 159 } 160 161 System.out.println("Checking accuracy of 3D complex FFT (double[][][] input)..."); 162 for (int i = 0; i < sizes3D.length; i++) { 163 DoubleFFT_3D fft3 = new DoubleFFT_3D(sizes3D[i], sizes3D[i], sizes3D[i]); 164 double err = 0.0; 165 double[][][] a = new double[sizes3D[i]][sizes3D[i]][2 * sizes3D[i]]; 166 IOUtils.fillMatrix_3D(sizes3D[i], sizes3D[i], 2 * sizes3D[i], a); 167 double[][][] b = new double[sizes3D[i]][sizes3D[i]][2 * sizes3D[i]]; 168 IOUtils.fillMatrix_3D(sizes3D[i], sizes3D[i], 2 * sizes3D[i], b); 169 fft3.complexForward(a); 170 fft3.complexInverse(a, true); 171 err = computeRMSE(a, b); 172 if (err > eps) { 173 System.err.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " 174 + err); 175 } else { 176 System.out.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " 177 + err); 178 } 179 a = null; 180 b = null; 181 fft3 = null; 182 System.gc(); 183 } 184 } 185 186 public static void checkAccuracyRealFFT_1D() { 187 System.out.println("Checking accuracy of 1D real FFT..."); 188 for (int i = 0; i < sizes1D.length; i++) { 189 DoubleFFT_1D fft = new DoubleFFT_1D(sizes1D[i]); 190 double err = 0.0; 191 double[] a = new double[sizes1D[i]]; 192 IOUtils.fillMatrix_1D(sizes1D[i], a); 193 double[] b = new double[sizes1D[i]]; 194 IOUtils.fillMatrix_1D(sizes1D[i], b); 195 fft.realForward(a); 196 fft.realInverse(a, true); 197 err = computeRMSE(a, b); 198 if (err > eps) { 199 System.err.println("\tsize = " + sizes1D[i] + ";\terror = " + err); 200 } else { 201 System.out.println("\tsize = " + sizes1D[i] + ";\terror = " + err); 202 } 203 a = null; 204 b = null; 205 fft = null; 206 System.gc(); 207 } 208 System.out.println("Checking accuracy of on 1D real forward full FFT..."); 209 for (int i = 0; i < sizes1D.length; i++) { 210 DoubleFFT_1D fft = new DoubleFFT_1D(sizes1D[i]); 211 double err = 0.0; 212 double[] a = new double[2 * sizes1D[i]]; 213 IOUtils.fillMatrix_1D(sizes1D[i], a); 214 double[] b = new double[2 * sizes1D[i]]; 215 for (int j = 0; j < sizes1D[i]; j++) { 216 b[2 * j] = a[j]; 217 } 218 fft.realForwardFull(a); 219 fft.complexInverse(a, true); 220 err = computeRMSE(a, b); 221 if (err > eps) { 222 System.err.println("\tsize = " + sizes1D[i] + ";\terror = " + err); 223 } else { 224 System.out.println("\tsize = " + sizes1D[i] + ";\terror = " + err); 225 } 226 a = null; 227 b = null; 228 fft = null; 229 System.gc(); 230 } 231 System.out.println("Checking accuracy of 1D real inverse full FFT..."); 232 for (int i = 0; i < sizes1D.length; i++) { 233 DoubleFFT_1D fft = new DoubleFFT_1D(sizes1D[i]); 234 double err = 0.0; 235 double[] a = new double[2 * sizes1D[i]]; 236 IOUtils.fillMatrix_1D(sizes1D[i], a); 237 double[] b = new double[2 * sizes1D[i]]; 238 for (int j = 0; j < sizes1D[i]; j++) { 239 b[2 * j] = a[j]; 240 } 241 fft.realInverseFull(a, true); 242 fft.complexForward(a); 243 err = computeRMSE(a, b); 244 if (err > eps) { 245 System.err.println("\tsize = " + sizes1D[i] + ";\terror = " + err); 246 } else { 247 System.out.println("\tsize = " + sizes1D[i] + ";\terror = " + err); 248 } 249 a = null; 250 b = null; 251 fft = null; 252 System.gc(); 253 } 254 255 } 256 257 public static void checkAccuracyRealFFT_2D() { 258 System.out.println("Checking accuracy of 2D real FFT (double[] input)..."); 259 for (int i = 0; i < sizes2D2.length; i++) { 260 DoubleFFT_2D fft2 = new DoubleFFT_2D(sizes2D2[i], sizes2D2[i]); 261 double err = 0.0; 262 double[] a = new double[sizes2D2[i] * sizes2D2[i]]; 263 IOUtils.fillMatrix_2D(sizes2D2[i], sizes2D2[i], a); 264 double[] b = new double[sizes2D2[i] * sizes2D2[i]]; 265 IOUtils.fillMatrix_2D(sizes2D2[i], sizes2D2[i], b); 266 fft2.realForward(a); 267 fft2.realInverse(a, true); 268 err = computeRMSE(a, b); 269 if (err > eps) { 270 System.err.println("\tsize = " + sizes2D2[i] + " x " + sizes2D2[i] + ";\terror = " + err); 271 } else { 272 System.out.println("\tsize = " + sizes2D2[i] + " x " + sizes2D2[i] + ";\terror = " + err); 273 } 274 a = null; 275 b = null; 276 fft2 = null; 277 System.gc(); 278 } 279 System.out.println("Checking accuracy of 2D real FFT (double[][] input)..."); 280 for (int i = 0; i < sizes2D2.length; i++) { 281 DoubleFFT_2D fft2 = new DoubleFFT_2D(sizes2D2[i], sizes2D2[i]); 282 double err = 0.0; 283 double[][] a = new double[sizes2D2[i]][sizes2D2[i]]; 284 IOUtils.fillMatrix_2D(sizes2D2[i], sizes2D2[i], a); 285 double[][] b = new double[sizes2D2[i]][sizes2D2[i]]; 286 IOUtils.fillMatrix_2D(sizes2D2[i], sizes2D2[i], b); 287 fft2.realForward(a); 288 fft2.realInverse(a, true); 289 err = computeRMSE(a, b); 290 if (err > eps) { 291 System.err.println("\tsize = " + sizes2D2[i] + " x " + sizes2D2[i] + ";\terror = " + err); 292 } else { 293 System.out.println("\tsize = " + sizes2D2[i] + " x " + sizes2D2[i] + ";\terror = " + err); 294 } 295 a = null; 296 b = null; 297 fft2 = null; 298 System.gc(); 299 } 300 System.out.println("Checking accuracy of 2D real forward full FFT (double[] input)..."); 301 for (int i = 0; i < sizes2D.length; i++) { 302 DoubleFFT_2D fft2 = new DoubleFFT_2D(sizes2D[i], sizes2D[i]); 303 double err = 0.0; 304 double[] a = new double[2 * sizes2D[i] * sizes2D[i]]; 305 IOUtils.fillMatrix_2D(sizes2D[i], sizes2D[i], a); 306 double[] b = new double[2 * sizes2D[i] * sizes2D[i]]; 307 for (int r = 0; r < sizes2D[i]; r++) { 308 for (int c = 0; c < sizes2D[i]; c++) { 309 b[r * 2 * sizes2D[i] + 2 * c] = a[r * sizes2D[i] + c]; 310 } 311 } 312 fft2.realForwardFull(a); 313 fft2.complexInverse(a, true); 314 err = computeRMSE(a, b); 315 if (err > eps) { 316 System.err.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err); 317 } else { 318 System.out.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err); 319 } 320 a = null; 321 b = null; 322 fft2 = null; 323 System.gc(); 324 } 325 System.out.println("Checking accuracy of 2D real forward full FFT (double[][] input)..."); 326 for (int i = 0; i < sizes2D.length; i++) { 327 DoubleFFT_2D fft2 = new DoubleFFT_2D(sizes2D[i], sizes2D[i]); 328 double err = 0.0; 329 double[][] a = new double[sizes2D[i]][2 * sizes2D[i]]; 330 IOUtils.fillMatrix_2D(sizes2D[i], sizes2D[i], a); 331 double[][] b = new double[sizes2D[i]][2 * sizes2D[i]]; 332 for (int r = 0; r < sizes2D[i]; r++) { 333 for (int c = 0; c < sizes2D[i]; c++) { 334 b[r][2 * c] = a[r][c]; 335 } 336 } 337 fft2.realForwardFull(a); 338 fft2.complexInverse(a, true); 339 err = computeRMSE(a, b); 340 if (err > eps) { 341 System.err.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err); 342 } else { 343 System.out.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err); 344 } 345 a = null; 346 b = null; 347 fft2 = null; 348 System.gc(); 349 } 350 System.out.println("Checking accuracy of 2D real inverse full FFT (double[] input)..."); 351 for (int i = 0; i < sizes2D.length; i++) { 352 DoubleFFT_2D fft2 = new DoubleFFT_2D(sizes2D[i], sizes2D[i]); 353 double err = 0.0; 354 double[] a = new double[2 * sizes2D[i] * sizes2D[i]]; 355 IOUtils.fillMatrix_2D(sizes2D[i], sizes2D[i], a); 356 double[] b = new double[2 * sizes2D[i] * sizes2D[i]]; 357 for (int r = 0; r < sizes2D[i]; r++) { 358 for (int c = 0; c < sizes2D[i]; c++) { 359 b[r * 2 * sizes2D[i] + 2 * c] = a[r * sizes2D[i] + c]; 360 } 361 } 362 fft2.realInverseFull(a, true); 363 fft2.complexForward(a); 364 err = computeRMSE(a, b); 365 if (err > eps) { 366 System.err.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err); 367 } else { 368 System.out.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err); 369 } 370 a = null; 371 b = null; 372 fft2 = null; 373 System.gc(); 374 } 375 System.out.println("Checking accuracy of 2D real inverse full FFT (double[][] input)..."); 376 for (int i = 0; i < sizes2D.length; i++) { 377 DoubleFFT_2D fft2 = new DoubleFFT_2D(sizes2D[i], sizes2D[i]); 378 double err = 0.0; 379 double[][] a = new double[sizes2D[i]][2 * sizes2D[i]]; 380 IOUtils.fillMatrix_2D(sizes2D[i], sizes2D[i], a); 381 double[][] b = new double[sizes2D[i]][2 * sizes2D[i]]; 382 for (int r = 0; r < sizes2D[i]; r++) { 383 for (int c = 0; c < sizes2D[i]; c++) { 384 b[r][2 * c] = a[r][c]; 385 } 386 } 387 fft2.realInverseFull(a, true); 388 fft2.complexForward(a); 389 err = computeRMSE(a, b); 390 if (err > eps) { 391 System.err.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err); 392 } else { 393 System.out.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err); 394 } 395 a = null; 396 b = null; 397 fft2 = null; 398 System.gc(); 399 } 400 401 } 402 403 public static void checkAccuracyRealFFT_3D() { 404 System.out.println("Checking accuracy of 3D real FFT (double[] input)..."); 405 for (int i = 0; i < sizes3D2.length; i++) { 406 DoubleFFT_3D fft3 = new DoubleFFT_3D(sizes3D2[i], sizes3D2[i], sizes3D2[i]); 407 double err = 0.0; 408 double[] a = new double[sizes3D2[i] * sizes3D2[i] * sizes3D2[i]]; 409 IOUtils.fillMatrix_3D(sizes3D2[i], sizes3D2[i], sizes3D2[i], a); 410 double[] b = new double[sizes3D2[i] * sizes3D2[i] * sizes3D2[i]]; 411 IOUtils.fillMatrix_3D(sizes3D2[i], sizes3D2[i], sizes3D2[i], b); 412 fft3.realForward(b); 413 fft3.realInverse(b, true); 414 err = computeRMSE(a, b); 415 if (err > eps) { 416 System.err.println("\tsize = " + sizes3D2[i] + " x " + sizes3D2[i] + " x " + sizes3D2[i] 417 + ";\t\terror = " + err); 418 } else { 419 System.out.println("\tsize = " + sizes3D2[i] + " x " + sizes3D2[i] + " x " + sizes3D2[i] 420 + ";\t\terror = " + err); 421 } 422 a = null; 423 b = null; 424 fft3 = null; 425 System.gc(); 426 } 427 System.out.println("Checking accuracy of 3D real FFT (double[][][] input)..."); 428 for (int i = 0; i < sizes3D2.length; i++) { 429 DoubleFFT_3D fft3 = new DoubleFFT_3D(sizes3D2[i], sizes3D2[i], sizes3D2[i]); 430 double err = 0.0; 431 double[][][] a = new double[sizes3D2[i]][sizes3D2[i]][sizes3D2[i]]; 432 IOUtils.fillMatrix_3D(sizes3D2[i], sizes3D2[i], sizes3D2[i], a); 433 double[][][] b = new double[sizes3D2[i]][sizes3D2[i]][sizes3D2[i]]; 434 IOUtils.fillMatrix_3D(sizes3D2[i], sizes3D2[i], sizes3D2[i], b); 435 fft3.realForward(b); 436 fft3.realInverse(b, true); 437 err = computeRMSE(a, b); 438 if (err > eps) { 439 System.err.println("\tsize = " + sizes3D2[i] + " x " + sizes3D2[i] + " x " + sizes3D2[i] 440 + ";\t\terror = " + err); 441 } else { 442 System.out.println("\tsize = " + sizes3D2[i] + " x " + sizes3D2[i] + " x " + sizes3D2[i] 443 + ";\t\terror = " + err); 444 } 445 a = null; 446 b = null; 447 fft3 = null; 448 System.gc(); 449 } 450 System.out.println("Checking accuracy of 3D real forward full FFT (double[] input)..."); 451 for (int i = 0; i < sizes3D.length; i++) { 452 DoubleFFT_3D fft3 = new DoubleFFT_3D(sizes3D[i], sizes3D[i], sizes3D[i]); 453 double err = 0.0; 454 double[] a = new double[2 * sizes3D[i] * sizes3D[i] * sizes3D[i]]; 455 IOUtils.fillMatrix_3D(sizes3D[i], sizes3D[i], sizes3D[i], a); 456 double[] b = new double[2 * sizes3D[i] * sizes3D[i] * sizes3D[i]]; 457 for (int s = 0; s < sizes3D[i]; s++) { 458 for (int r = 0; r < sizes3D[i]; r++) { 459 for (int c = 0; c < sizes3D[i]; c++) { 460 b[s * 2 * sizes3D[i] * sizes3D[i] + r * 2 * sizes3D[i] + 2 * c] = a[s * sizes3D[i] * sizes3D[i] 461 + r * sizes3D[i] + c]; 462 } 463 } 464 } 465 fft3.realForwardFull(a); 466 fft3.complexInverse(a, true); 467 err = computeRMSE(a, b); 468 if (err > eps) { 469 System.err.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " 470 + err); 471 } else { 472 System.out.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " 473 + err); 474 } 475 a = null; 476 b = null; 477 fft3 = null; 478 System.gc(); 479 } 480 481 System.out.println("Checking accuracy of 3D real forward full FFT (double[][][] input)..."); 482 for (int i = 0; i < sizes3D.length; i++) { 483 DoubleFFT_3D fft3 = new DoubleFFT_3D(sizes3D[i], sizes3D[i], sizes3D[i]); 484 double err = 0.0; 485 double[][][] a = new double[sizes3D[i]][sizes3D[i]][2 * sizes3D[i]]; 486 IOUtils.fillMatrix_3D(sizes3D[i], sizes3D[i], sizes3D[i], a); 487 double[][][] b = new double[sizes3D[i]][sizes3D[i]][2 * sizes3D[i]]; 488 for (int s = 0; s < sizes3D[i]; s++) { 489 for (int r = 0; r < sizes3D[i]; r++) { 490 for (int c = 0; c < sizes3D[i]; c++) { 491 b[s][r][2 * c] = a[s][r][c]; 492 } 493 } 494 } 495 fft3.realForwardFull(a); 496 fft3.complexInverse(a, true); 497 err = computeRMSE(a, b); 498 if (err > eps) { 499 System.err.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " 500 + err); 501 } else { 502 System.out.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " 503 + err); 504 } 505 a = null; 506 b = null; 507 fft3 = null; 508 System.gc(); 509 } 510 511 System.out.println("Checking accuracy of 3D real inverse full FFT (double[] input)..."); 512 for (int i = 0; i < sizes3D.length; i++) { 513 DoubleFFT_3D fft3 = new DoubleFFT_3D(sizes3D[i], sizes3D[i], sizes3D[i]); 514 double err = 0.0; 515 double[] a = new double[2 * sizes3D[i] * sizes3D[i] * sizes3D[i]]; 516 IOUtils.fillMatrix_3D(sizes3D[i], sizes3D[i], sizes3D[i], a); 517 double[] b = new double[2 * sizes3D[i] * sizes3D[i] * sizes3D[i]]; 518 for (int s = 0; s < sizes3D[i]; s++) { 519 for (int r = 0; r < sizes3D[i]; r++) { 520 for (int c = 0; c < sizes3D[i]; c++) { 521 b[s * 2 * sizes3D[i] * sizes3D[i] + r * 2 * sizes3D[i] + 2 * c] = a[s * sizes3D[i] * sizes3D[i] 522 + r * sizes3D[i] + c]; 523 } 524 } 525 } 526 fft3.realInverseFull(a, true); 527 fft3.complexForward(a); 528 err = computeRMSE(a, b); 529 if (err > eps) { 530 System.err.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " 531 + err); 532 } else { 533 System.out.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " 534 + err); 535 } 536 a = null; 537 b = null; 538 fft3 = null; 539 System.gc(); 540 } 541 System.out.println("Checking accuracy of 3D real inverse full FFT (double[][][] input)..."); 542 for (int i = 0; i < sizes3D.length; i++) { 543 DoubleFFT_3D fft3 = new DoubleFFT_3D(sizes3D[i], sizes3D[i], sizes3D[i]); 544 double err = 0.0; 545 double[][][] a = new double[sizes3D[i]][sizes3D[i]][2 * sizes3D[i]]; 546 IOUtils.fillMatrix_3D(sizes3D[i], sizes3D[i], sizes3D[i], a); 547 double[][][] b = new double[sizes3D[i]][sizes3D[i]][2 * sizes3D[i]]; 548 for (int s = 0; s < sizes3D[i]; s++) { 549 for (int r = 0; r < sizes3D[i]; r++) { 550 for (int c = 0; c < sizes3D[i]; c++) { 551 b[s][r][2 * c] = a[s][r][c]; 552 } 553 } 554 } 555 fft3.realInverseFull(a, true); 556 fft3.complexForward(a); 557 err = computeRMSE(a, b); 558 if (err > eps) { 559 System.err.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " 560 + err); 561 } else { 562 System.out.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " 563 + err); 564 } 565 a = null; 566 b = null; 567 fft3 = null; 568 System.gc(); 569 } 570 } 571 572 private static double computeRMSE(double[] a, double[] b) { 573 if (a.length != b.length) { 574 throw new IllegalArgumentException("Arrays are not the same size."); 575 } 576 double rms = 0; 577 double tmp; 578 for (int i = 0; i < a.length; i++) { 579 tmp = (a[i] - b[i]); 580 rms += tmp * tmp; 581 } 582 return Math.sqrt(rms / a.length); 583 } 584 585 private static double computeRMSE(double[][] a, double[][] b) { 586 if (a.length != b.length || a[0].length != b[0].length) { 587 throw new IllegalArgumentException("Arrays are not the same size."); 588 } 589 double rms = 0; 590 double tmp; 591 for (int r = 0; r < a.length; r++) { 592 for (int c = 0; c < a[0].length; c++) { 593 tmp = (a[r][c] - b[r][c]); 594 rms += tmp * tmp; 595 } 596 } 597 return Math.sqrt(rms / (a.length * a[0].length)); 598 } 599 600 private static double computeRMSE(double[][][] a, double[][][] b) { 601 if (a.length != b.length || a[0].length != b[0].length || a[0][0].length != b[0][0].length) { 602 throw new IllegalArgumentException("Arrays are not the same size."); 603 } 604 double rms = 0; 605 double tmp; 606 for (int s = 0; s < a.length; s++) { 607 for (int r = 0; r < a[0].length; r++) { 608 for (int c = 0; c < a[0][0].length; c++) { 609 tmp = (a[s][r][c] - b[s][r][c]); 610 rms += tmp * tmp; 611 } 612 } 613 } 614 return Math.sqrt(rms / (a.length * a[0].length * a[0][0].length)); 615 } 616 617 public static void main(String[] args) { 618 checkAccuracyComplexFFT_1D(); 619 checkAccuracyRealFFT_1D(); 620 checkAccuracyComplexFFT_2D(); 621 checkAccuracyRealFFT_2D(); 622 checkAccuracyComplexFFT_3D(); 623 checkAccuracyRealFFT_3D(); 624 System.exit(0); 625 } 626}