jnt.FFT.RealDoubleFFT_Radix2 Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of sspace Show documentation
Show all versions of sspace Show documentation
The S-Space Package is a Natural Language Processing library for
distributional semantics representations. Distributional semantics
representations model the meaning of words, phrases, and sentences as high
dimensional vectors or probability distributions. The library includes common
algorithms such as Latent Semantic Analysis, Random Indexing, and Latent
Dirichlet Allocation. The S-Space package also includes software libraries
for matrices, vectors, graphs, and numerous clustering
algorithms.
The newest version!
package jnt.FFT;
/** Computes FFT's of real, double precision data where n is an integral power of 2.
* The physical layout of the mathematical data d[i] in the array data is as follows:
*
* d[i] = data[i0 + stride*i]
*
* The FFT (D[i]) of real data (d[i]) is complex, but restricted by symmetry:
*
* D[n-i] = conj(D[i])
*
* It turns out that there are still n `independent' values, so the transformation
* can still be carried out in-place.
* For RealDoubleFFT_Radix2, the correspondence is as follows:
*
*Logical Physical Re(D[0]) = data[0] Im(D[0]) = 0 Re(D[1]) = data[1] Im(D[1]) = data[n-1] ... ... Re(D[k]) = data[k] Im(D[k]) = data[n-k] ... ... Re(D[n/2]) = data[n/2] Im(D[n/2]) = 0 ... ... Re(D[n-k]) = data[k] Im(D[n-k]) = -data[n-k] ... ... Re(D[n-1]) = data[1] Im(D[n-1]) = -data[n-1]
*
* @author Bruce R. Miller [email protected]
* @author Contribution of the National Institute of Standards and Technology,
* @author not subject to copyright.
* @author Derived from GSL (Gnu Scientific Library)
* @author GSL's FFT Code by Brian Gough [email protected]
* @author Since GSL is released under
* @author GPL,
* @author this package must also be.
*/
public class RealDoubleFFT_Radix2 extends RealDoubleFFT {
private int logn;
/** Create an FFT for transforming n points of real, double precision data.
* n must be an integral power of 2. */
public RealDoubleFFT_Radix2(int n){
/* make sure that n is a power of 2 */
super(n);
logn = Factorize.log2(n);
if (logn < 0)
throw new IllegalArgumentException(n+" is not a power of 2");
}
/** Compute the Fast Fourier Transform of data leaving the result in data.
* See {@link Radix2 Transform Layout} for description of
* the resulting data layout.*/
public void transform (double data[], int i0, int stride) {
checkData(data,i0,stride);
int p, p_1, q;
if (n == 1) return; /* identity operation */
/* bit reverse the ordering of input data for decimation in time algorithm */
bitreverse(data, i0, stride);
/* apply fft recursion */
p = 1; q = n ;
for (int i = 1; i <= logn; i++) {
int a, b;
p_1 = p ;
p = 2 * p ;
q = q / 2 ;
/* a = 0 */
for (b = 0; b < q; b++) {
double t0_real = data[i0+stride*b*p] + data[i0+stride*(b*p + p_1)];
double t1_real = data[i0+stride*b*p] - data[i0+stride*(b*p + p_1)];
data[i0+stride*b*p] = t0_real;
data[i0+stride*(b*p + p_1)] = t1_real;
}
/* a = 1 ... p_{i-1}/2 - 1 */
{
double w_real = 1.0;
double w_imag = 0.0;
double theta = - 2.0 * Math.PI / p;
double s = Math.sin(theta);
double t = Math.sin(theta / 2.0);
double s2 = 2.0 * t * t;
for (a = 1; a < (p_1)/2; a++) {
/* trignometric recurrence for w-> exp(i theta) w */
{
double tmp_real = w_real - s * w_imag - s2 * w_real;
double tmp_imag = w_imag + s * w_real - s2 * w_imag;
w_real = tmp_real;
w_imag = tmp_imag;
}
for (b = 0; b < q; b++) {
double z0_real = data[i0+stride*(b*p + a)];
double z0_imag = data[i0+stride*(b*p + p_1 - a)];
double z1_real = data[i0+stride*(b*p + p_1 + a)];
double z1_imag = data[i0+stride*(b*p + p - a)];
/* t0 = z0 + w * z1 */
data[i0+stride*(b*p + a)] =z0_real + w_real * z1_real - w_imag * z1_imag;
data[i0+stride*(b*p + p - a)]=z0_imag + w_real * z1_imag + w_imag * z1_real;
/* t1 = -(z0 - w * z1) */
data[i0+stride*(b*p + p_1 - a)]=z0_real - w_real * z1_real + w_imag * z1_imag;
data[i0+stride*(b*p + p_1 + a)]=-(z0_imag - w_real * z1_imag - w_imag * z1_real);
}
}
}
if (p_1 > 1) {
for (b = 0; b < q; b++) {
/* a = p_{i-1}/2 */
data[i0+stride*(b*p + p - p_1/2)] *= -1 ;
}}
}
}
/** Compute the (unnomalized) inverse FFT of data, leaving it in place.
* The data must be in the same arrangement as that produced by
{@link #transform transform}. */
public void backtransform (double data[], int i0, int stride) {
checkData(data,i0,stride);
int p, p_1, q;
if (n == 1) return; /* identity operation */
/* apply fft recursion */
p = n; q = 1 ; p_1 = n/2 ;
for (int i = 1; i <= logn; i++) {
int a, b;
/* a = 0 */
for (b = 0; b < q; b++) {
double z0 = data[i0+stride*b*p];
double z1 = data[i0+stride*(b*p + p_1)];
data[i0+stride*b*p] = z0 + z1 ;
data[i0+stride*(b*p + p_1)] = z0 - z1 ; }
/* a = 1 ... p_{i-1}/2 - 1 */
{
double w_real = 1.0;
double w_imag = 0.0;
double theta = 2.0 * Math.PI / p;
double s = Math.sin(theta);
double t = Math.sin(theta / 2.0);
double s2 = 2.0 * t * t;
for (a = 1; a < (p_1)/2; a++) {
/* trignometric recurrence for w-> exp(i theta) w */
double tmp_real = w_real - s * w_imag - s2 * w_real;
double tmp_imag = w_imag + s * w_real - s2 * w_imag;
w_real = tmp_real;
w_imag = tmp_imag;
for (b = 0; b < q; b++) {
double z0_real = data[i0+stride*(b*p + a)];
double z0_imag = data[i0+stride*(b*p + p - a)];
double z1_real = data[i0+stride*(b*p + p_1 - a)];
double z1_imag = -data[i0+stride*(b*p + p_1 + a)];
/* t0 = z0 + z1 */
data[i0+stride*(b*p + a)] = z0_real + z1_real;
data[i0+stride*(b*p + p_1 - a)] = z0_imag + z1_imag;
/* t1 = (z0 - z1) */
double t1_real = z0_real - z1_real;
double t1_imag = z0_imag - z1_imag;
data[i0+stride*(b*p + p_1 + a)] = (w_real * t1_real - w_imag * t1_imag) ;
data[i0+stride*(b*p + p - a)] = (w_real * t1_imag + w_imag * t1_real) ;
}
}
}
if (p_1 > 1) {
for (b = 0; b < q; b++) {
data[i0+stride*(b*p + p_1/2)] *= 2 ;
data[i0+stride*(b*p + p_1 + p_1/2)] *= -2 ;
}
}
p_1 = p_1 / 2 ;
p = p / 2 ;
q = q * 2 ;
}
/* bit reverse the ordering of output data for decimation in
frequency algorithm */
bitreverse(data, i0, stride);
}
/** Return data in wraparound order.
* @see wraparound format */
public double[] toWraparoundOrder(double data[]){
return toWraparoundOrder(data,0,1); }
/** Return data in wraparound order.
* i0 and stride are used to traverse data; the new array is in
* packed (i0=0, stride=1) format.
* @see wraparound format */
public double[] toWraparoundOrder(double data[], int i0, int stride) {
checkData(data,i0,stride);
double newdata[] = new double[2*n];
int nh = n/2;
newdata[0] = data[i0];
newdata[1] = 0.0;
newdata[n] = data[i0+stride*nh];
newdata[n+1] = 0.0;
for(int i=1; i