001/******************************************************************************* 002 * This software is provided as a supplement to the authors' textbooks on digital 003 * image processing published by Springer-Verlag in various languages and editions. 004 * Permission to use and distribute this software is granted under the BSD 2-Clause 005 * "Simplified" License (see http://opensource.org/licenses/BSD-2-Clause). 006 * Copyright (c) 2006-2023 Wilhelm Burger, Mark J. Burge. All rights reserved. 007 * Visit https://imagingbook.com for additional details. 008 ******************************************************************************/ 009package imagingbook.spectral.dft; 010 011import org.jtransforms.fft.DoubleFFT_1D; 012import org.jtransforms.fft.FloatFFT_1D; 013 014/** 015 * <p> 016 * FFT (fast) implementation of the DFT, based on the JTransforms package by 017 * Piotr Wendykier (see <a href="https://github.com/wendykierp/JTransforms"> 018 * https://github.com/wendykierp/JTransforms</a>). Note that this class has no 019 * public constructor - instantiate sub-class {@link Dft1dFast.Float} or 020 * {@link Dft1dFast.Double} instead, as shown below. See Sec. 18.4.2 of [1] for 021 * additional details. 022 * </p> 023 * <p> 024 * Usage example (for {@code float} data): 025 * </p> 026 * <pre> 027 * float[] re = {1, 2, 3, 4, 5, 3, 0}; 028 * float[] im = {0, 0, 0, 0, 0, 0, 0}; 029 * Dft1d.Float dft = new Dft1dFast.Float(re.length); 030 * dct.forward(re, im); // re/im now is the DFT spectrum 031 * dct.inverse(re, im); // re/im now is the original signal 032 * ... 033 * </pre> 034 * <p> 035 * [1] W. Burger, M.J. Burge, <em>Digital Image Processing – An Algorithmic 036 * Introduction</em>, 3rd ed, Springer (2022). 037 * </p> 038 */ 039public abstract class Dft1dFast extends Dft1dImp { 040 041 private Dft1dFast(int size, ScalingMode sm) { 042 super(size, sm); 043 } 044 045 // ---------------------------------------------------------------------- 046 047 /** 048 * One-dimensional FFT implementation using {@code float} data. 049 */ 050 public static class Float extends Dft1dFast implements Dft1d.Float { 051 052 private final float[] A; // temporary array for FFT composed of re/im values 053 private final FloatFFT_1D fft; 054 055 /** 056 * Constructor using a specific scaling mode. 057 * @param M the size of the data vectors 058 * @param sm the scaling mode 059 */ 060 public Float(int M, ScalingMode sm) { 061 super(M, sm); 062 this.A = new float[2 * M]; 063 this.fft = new FloatFFT_1D(M); 064 } 065 066 /** 067 * Constructor using the default scaling mode. 068 * @param M the size of the data vectors 069 */ 070 public Float(int M) { 071 this(M, ScalingMode.DEFAULT); 072 } 073 074 @Override 075 public void forward(float[] gRe, float[] gIm) { 076 transform(gRe, gIm, true); 077 } 078 079 @Override 080 public void inverse(float[] GRe, float[] GIm) { 081 transform(GRe, GIm, false); 082 } 083 084 @Override 085 public void transform(float[] inRe, float[] inIm, boolean forward) { 086 checkSize(inRe, inIm); 087 final float scale = (float) sm.getScale(M, forward); 088 composeA(inRe, inIm, A); 089 if (forward) 090 fft.complexForward(A); 091 else 092 fft.complexInverse(A, false); 093 decomposeA(A, inRe, inIm, scale); 094 } 095 096 // (re, im) -> A 097 private void composeA(float[] re, float[] im, float[] A) { 098 for (int i = 0; i < M; i++) { 099 A[2*i] = re[i]; 100 A[2*i + 1] = im[i]; 101 } 102 } 103 104 // A -> (re, im) 105 private void decomposeA(float[] A, float[] re, float[] im, float scale) { 106 for (int i = 0; i < M; i++) { 107 re[i] = A[2*i] * scale; 108 im[i] = A[2*i + 1] * scale; 109 } 110 } 111 } 112 113 // ---------------------------------------------------------------------- 114 115 /** 116 * One-dimensional FFT implementation using {@code double} data. 117 */ 118 public static class Double extends Dft1dFast implements Dft1d.Double { 119 120 private final double[] A; // temporary array for FFT composed of re/im values 121 private final DoubleFFT_1D fft; 122 123 /** 124 * Constructor using a specific scaling mode. 125 * @param M the size of the data vectors 126 * @param sm the scaling mode 127 */ 128 public Double(int M, ScalingMode sm) { 129 super(M, sm); 130 this.A = new double[2 * M]; 131 this.fft = new DoubleFFT_1D(M); 132 } 133 134 /** 135 * Constructor using the default scaling mode. 136 * @param M the size of the data vectors 137 */ 138 public Double(int M) { 139 this(M, ScalingMode.DEFAULT); 140 } 141 142 @Override 143 public void forward(double[] gRe, double[] gIm) { 144 transform(gRe, gIm, true); 145 } 146 147 @Override 148 public void inverse(double[] GRe, double[] GIm) { 149 transform(GRe, GIm, false); 150 } 151 152 @Override 153 public void transform(double[] inRe, double[] inIm, boolean forward) { 154 checkSize(inRe, inIm); 155 final double scale = sm.getScale(M, forward); 156 composeA(inRe, inIm, A); 157 if (forward) 158 fft.complexForward(A); 159 else 160 fft.complexInverse(A, false); 161 decomposeA(A, inRe, inIm, scale); 162 } 163 164 // (re, im) -> A 165 private void composeA(double[] re, double[] im, double[] A) { 166 for (int i = 0; i < M; i++) { 167 A[2*i] = re[i]; 168 A[2*i + 1] = im[i]; 169 } 170 } 171 172 // A -> (re, im) 173 private void decomposeA(double[] A, double[] re, double[] im, double scale) { 174 for (int i = 0; i < M; i++) { 175 re[i] = A[2*i] * scale; 176 im[i] = A[2*i + 1] * scale; 177 } 178 } 179 } 180 181 // ---------------------------------------------------------------------- 182 183 /* 184 * Direct implementation of the one-dimensional DFT for arbitrary signal lengths. 185 * This DFT uses the same definition as Mathematica. Example: 186 * Fourier[{1, 2, 3, 4, 5, 6, 7, 8}, FourierParameters -> {0, -1}]: 187 {12.7279 + 0. i, 188 -1.41421 + 3.41421 i, 189 -1.41421 + 1.41421 i, 190 -1.41421 + 0.585786 i, 191 -1.41421 + 0. i, 192 -1.41421 - 0.585786 i, 193 -1.41421 - 1.41421 i, 194 -1.41421 - 3.41421 i} 195 */ 196 197 //test example 198// public static void main(String[] args) { 199// 200// System.out.println("******************** Float test (FFT) ********************"); 201// { 202// float[] re = { 1, 2, 3, 4, 5, 6, 7, 8 }; 203// float[] im = new float[re.length]; 204// 205// System.out.println("original signal:"); 206// System.out.println("gRe = " + Matrix.toString(re)); 207// System.out.println("gIm = " + Matrix.toString(im)); 208// 209// Dft1d.Float dft = new Dft1dFast.Float(re.length); 210// dft.forward(re, im); 211// 212// System.out.println("DFT spectrum:"); 213// System.out.println("GRe = " + Matrix.toString(re)); 214// System.out.println("GIm = " + Matrix.toString(im)); 215// 216// dft.inverse(re, im); 217// 218// System.out.println("reconstructed signal:"); 219// System.out.println("gRe' = " + Matrix.toString(re)); 220// System.out.println("gIm' = " + Matrix.toString(im)); 221// } 222// 223// System.out.println(); 224// System.out.println("******************** Double test (FFT) ********************"); 225// 226// { 227// double[] re = { 1, 2, 3, 4, 5, 6, 7, 8 }; 228// double[] im = new double[re.length]; 229// 230// System.out.println("original signal:"); 231// System.out.println("gRe = " + Matrix.toString(re)); 232// System.out.println("gIm = " + Matrix.toString(im)); 233// 234// Dft1d.Double dft = new Dft1dFast.Double(re.length); 235// dft.forward(re, im); 236// 237// System.out.println("DFT spectrum:"); 238// System.out.println("GRe = " + Matrix.toString(re)); 239// System.out.println("GIm = " + Matrix.toString(im)); 240// 241// dft.inverse(re, im); 242// 243// System.out.println("reconstructed signal:"); 244// System.out.println("gRe' = " + Matrix.toString(re)); 245// System.out.println("gIm' = " + Matrix.toString(im)); 246// } 247// } 248 249}