gov.fnal.controls.tools.math
Class FourierTransforms

java.lang.Object
  |
  +--gov.fnal.controls.tools.math.FourierTransforms

public class FourierTransforms

extends Object

Java implementation of FFT package from Tokyo University

Version:

1 This class is actually a set of static functions and static data ...

 Original files can be found at http://momonga.t.u-tokyo.ac.jp/~ooura/fft.html
 Actually I want to implement only
      fftsg.c    : FFT Package in C       - Fast Version   III (Split-Radix)
 As is written in comments, this file is optimized for computers with multylevel cashe
 I actually don't know, what kind of optimization JVM have and how it would work
 in Java. I would place benchmark's here ... I can't benchmark the original C code
 on this computer. At least now ...
 All original mathematical comments:
 Fast Fourier/Cosine/Sine Transform
  dimension   :one
  data length :power of 2
  decimation  :frequency
  radix       :split-radix
  data        :inplace
  table       :use
 functions
   cdft: Complex Discrete Fourier Transform
   _cdft: Complex Discrete Fourier Transform with working area parameters
   rdft: Real Discrete Fourier Transform
   _rdft: Real Discrete Fourier Transform with working area parameters
   ddct: Discrete Cosine Transform
   _ddct: Discrete Cosine Transform with working area parameters
   ddst: Discrete Sine Transform
   _ddst: Discrete Sine Transform with working area parameters
   dfct: Cosine Transform of RDFT (Real Symmetric DFT)
   _dfct: Cosine Transform of RDFT (Real Symmetric DFT) with working area parameters
   dfst: Sine Transform of RDFT (Real Anti-symmetric DFT)
  _dfst: Sine Transform of RDFT (Real Anti-symmetric DFT) with working area parameters
 Appendix :
   The cos/sin table are recalculated when the larger table required.
   Its are static class variables;

Author:

Timofei Bolshakov

Field Summary

protected static TreeMap	c_tables Map Integer n -> double [] c table for transform
static int	CDFT_RECURSIVE_N
protected static int[]	ip IP working area
protected static double[]	t Working area for dfct & dfst
protected static TreeMap	w_tables Map Integer n -> double [] w table for transform

Constructor Summary

FourierTransforms()

Method Summary

static void	_cdft(int n, int isgn, double[] a, int[] ip) The same, as cdft but with working area parameter Can be used in multy-threaded application (each thread have to have own working area) int [] ip - length sgrt(n/2)+1; (@see cdft)
static void	_ddct(int n, int isgn, double[] a, int[] ip) The same, as ddct but with working area parameter Can be used in multy-threaded application (each thread have to have own working area) int [] ip - length sgrt(n/2)+1; (@see ddct)
static void	_ddst(int n, int isgn, double[] a, int[] ip) The same, as ddst but with working area parameter Can be used in multy-threaded application (each thread have to have own working area) int [] ip - length sgrt(n/2)+1; (@see ddst)
static void	_dfct(int n, double[] a, int[] ip, double[] t) The same, as dfct but with working area parameter Can be used in multy-threaded application (each thread have to have own working area) int [] ip - length sgrt(n/2)+1; double [] t - length n/2+1 (@see dfct)
static void	_dfst(int n, double[] a, int[] ip, double[] t) The same, as dfst but with working area parameter Can be used in multy-threaded application (each thread have to have own working area) int [] ip - length sgrt(n/2)+1; double [] t - length n/2+1 (@see dfst)
static void	_rdft(int n, int isgn, double[] a, int[] ip) The same, as rdft but with working area parameter Can be used in multy-threaded application (each thread have to have own working area) int [] ip - length sgrt(n/2)+1; (@see rdft)
static void	_sym_cdft(int n, int isgn, double[] a, int[] ip) The same, as sym_cdft but with working area parameter Can be used in multy-threaded application (each thread have to have own working area) int [] ip - length sgrt(n/2)+1; (@see cdft)
static void	_sym_rdft(int n, int isgn, double[] a, int[] ip) The same, as sym_rdft but with working area parameter Can be used in multy-threaded application (each thread have to have own working area) int [] ip - length sgrt(n/2)+1; (@see rdft)
static void	cdft(int n, int isgn, double[] a) cdft: Complex Discrete Fourier Transform
static void	compute_w_table(int n)
static void	ddct(int n, int isgn, double[] a) ddct: Discrete Cosine Transform
static void	ddst(int n, int isgn, double[] a) ddst: Discrete Sine Transform
static void	dfct(int n, double[] a) dfct: Cosine Transform of RDFT (Real Symmetric DFT)
static void	dfst(int n, double[] a) dfst: Sine Transform of RDFT (Real Anti-symmetric DFT)
static void	main(String[] args) Check all transforms doing direct and inverse transform, check and print the error
static void	rdft(int n, int isgn, double[] a) rdft: Real Discrete Fourier Transform
static void	sym_cdft(int n, int isgn, double[] a) The same, as cdft but normalized on koefficient sqrt(2/n), So, for sym_cdft(n, 1, a); inverse is sym_cdft(n, -1, a); (@see cdft)
static void	sym_rdft(int n, int isgn, double[] a) The same, as rdft but normalized on koefficient sqrt(2/n), So, for sym_rdft(n, 1, a); inverse is sym_rdft(n, -1, a); (@see rdft)

Methods inherited from class java.lang.Object

clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait

Field Detail

c_tables

protected static TreeMap c_tables

Map Integer n -> double [] c table for transform

w_tables

protected static TreeMap w_tables

Map Integer n -> double [] w table for transform

ip

protected static int[] ip

IP working area

t

protected static double[] t

Working area for dfct & dfst

CDFT_RECURSIVE_N

public static int CDFT_RECURSIVE_N

Constructor Detail

FourierTransforms

public FourierTransforms()

Method Detail

compute_w_table

public static void compute_w_table(int n)

dfst

public static void dfst(int n,
                        double[] a)

dfst: Sine Transform of RDFT (Real Anti-symmetric DFT)

-------- Sine Transform of RDFT (Real Anti-symmetric DFT) --------
[definition]
S[k] = sum_j=1^n-1 a[j]*sin(pi*j*k/n), 0= 2, n = power of 2
a[0...n-1]     :input/output data (double *)
output data
a[k] = S[k], 0

_dfst

public static void _dfst(int n,
                         double[] a,
                         int[] ip,
                         double[] t)

The same, as dfst but with working area parameter Can be used in multy-threaded application (each thread have to have own working area) int [] ip - length sgrt(n/2)+1; double [] t - length n/2+1 (@see dfst)

dfct

public static void dfct(int n,
                        double[] a)

dfct: Cosine Transform of RDFT (Real Symmetric DFT)

-------- Cosine Transform of RDFT (Real Symmetric DFT) --------
[definition]
C[k] = sum_j=0^n a[j]*cos(pi*j*k/n), 0<=k<=n
[usage]
dfct(n, a);
[parameters]
n              :data length - 1 (int)
n >= 2, n = power of 2
a[0...n]       :input/output data (double *)
output data
a[k] = C[k], 0<=k<=n
[remark]
Inverse of
a[0] *= 0.5;
a[n] *= 0.5;
dfct(n, a);
is
a[0] *= 0.5;
a[n] *= 0.5;
dfct(n, a);
for (j = 0; j <= n; j++) {
a[j] *= 2.0 / n;
}

_dfct

public static void _dfct(int n,
                         double[] a,
                         int[] ip,
                         double[] t)

The same, as dfct but with working area parameter Can be used in multy-threaded application (each thread have to have own working area) int [] ip - length sgrt(n/2)+1; double [] t - length n/2+1 (@see dfct)

ddst

public static void ddst(int n,
                        int isgn,
                        double[] a)

ddst: Discrete Sine Transform

-------- DST (Discrete Sine Transform) / Inverse of DST --------
[definition]
 IDST (excluding scale)
S[k] = sum_j=1^n A[j]*sin(pi*j*(k+1/2)/n), 0<=k DST
S[k] = sum_j=0^n-1 a[j]*sin(pi*(j+1/2)*k/n), 0
ddst(n, 1, a);

ddst(n, -1, a);
[parameters]
n              :data length (int)
n >= 2, n = power of 2
a[0...n-1]     :input/output data (double *)

input data
a[j] = A[j], 0
output data
a[k] = S[k], 0

_ddst

public static void _ddst(int n,
                         int isgn,
                         double[] a,
                         int[] ip)

The same, as ddst but with working area parameter Can be used in multy-threaded application (each thread have to have own working area) int [] ip - length sgrt(n/2)+1; (@see ddst)

ddct

public static void ddct(int n,
                        int isgn,
                        double[] a)

ddct: Discrete Cosine Transform

-------- DCT (Discrete Cosine Transform) / Inverse of DCT --------
[definition]
 IDCT (excluding scale)
C[k] = sum_j=0^n-1 a[j]*cos(pi*j*(k+1/2)/n), 0<=k DCT
C[k] = sum_j=0^n-1 a[j]*cos(pi*(j+1/2)*k/n), 0<=k
ddct(n, 1, a);

ddct(n, -1, a);
[parameters]
n              :data length (int)
n >= 2, n = power of 2
a[0...n-1]     :input/output data (double *)
output data
a[k] = C[k], 0<=k

_ddct

public static void _ddct(int n,
                         int isgn,
                         double[] a,
                         int[] ip)

The same, as ddct but with working area parameter Can be used in multy-threaded application (each thread have to have own working area) int [] ip - length sgrt(n/2)+1; (@see ddct)

sym_rdft

public static void sym_rdft(int n,
                            int isgn,
                            double[] a)

The same, as rdft but normalized on koefficient sqrt(2/n), So, for sym_rdft(n, 1, a); inverse is sym_rdft(n, -1, a); (@see rdft)

_sym_rdft

public static void _sym_rdft(int n,
                             int isgn,
                             double[] a,
                             int[] ip)

The same, as sym_rdft but with working area parameter Can be used in multy-threaded application (each thread have to have own working area) int [] ip - length sgrt(n/2)+1; (@see rdft)

rdft

public static void rdft(int n,
                        int isgn,
                        double[] a)

rdft: Real Discrete Fourier Transform

-------- Real DFT / Inverse of Real DFT --------
[definition]
 RDFT
R[k] = sum_j=0^n-1 a[j]*cos(2*pi*j*k/n), 0<=k<=n/2
I[k] = sum_j=0^n-1 a[j]*sin(2*pi*j*k/n), 0 IRDFT (excluding scale)
a[k] = (R[0] + R[n/2]*cos(pi*k))/2 +
sum_j=1^n/2-1 R[j]*cos(2*pi*j*k/n) +
sum_j=1^n/2-1 I[j]*sin(2*pi*j*k/n), 0<=k
rdft(n, 1, a);

rdft(n, -1, a);
[parameters]
n              :data length (int)
n >= 2, n = power of 2
a[0...n-1]     :input/output data (double *)

output data
a[2*k] = R[k], 0<=k
input data
a[2*j] = R[j], 0<=j

_rdft

public static void _rdft(int n,
                         int isgn,
                         double[] a,
                         int[] ip)

The same, as rdft but with working area parameter Can be used in multy-threaded application (each thread have to have own working area) int [] ip - length sgrt(n/2)+1; (@see rdft)

sym_cdft

public static void sym_cdft(int n,
                            int isgn,
                            double[] a)

The same, as cdft but normalized on koefficient sqrt(2/n), So, for sym_cdft(n, 1, a); inverse is sym_cdft(n, -1, a); (@see cdft)

_sym_cdft

public static void _sym_cdft(int n,
                             int isgn,
                             double[] a,
                             int[] ip)

The same, as sym_cdft but with working area parameter Can be used in multy-threaded application (each thread have to have own working area) int [] ip - length sgrt(n/2)+1; (@see cdft)

cdft

public static void cdft(int n,
                        int isgn,
                        double[] a)

cdft: Complex Discrete Fourier Transform

[definition]

X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k
X[k] = sum_j=0^n-1 x[j]*exp(-2*pi*i*j*k/n), 0<=k
cdft(2*n, 1, a);

cdft(2*n, -1, a);
[parameters]
2*n            :data length (int)
n >= 1, n = power of 2
a[0...2*n-1]   :input/output data (double *)
input data
a[2*j] = Re(x[j]),
a[2*j+1] = Im(x[j]), 0<=j

_cdft

public static void _cdft(int n,
                         int isgn,
                         double[] a,
                         int[] ip)

The same, as cdft but with working area parameter Can be used in multy-threaded application (each thread have to have own working area) int [] ip - length sgrt(n/2)+1; (@see cdft)

main

public static void main(String[] args)

Check all transforms doing direct and inverse transform, check and print the error