jnt.FFT.RealFloatFFT_Radix2 Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of sspace-wordsi Show documentation
Show all versions of sspace-wordsi Show documentation
The S-Space Package is a collection of algorithms for building
Semantic Spaces as well as a highly-scalable library for designing new
distributional semantics algorithms. Distributional algorithms process text
corpora and represent the semantic for words as high dimensional feature
vectors. This package also includes matrices, vectors, and numerous
clustering algorithms. These approaches are known by many names, such as
word spaces, semantic spaces, or distributed semantics and rest upon the
Distributional Hypothesis: words that appear in similar contexts have
similar meanings.
The newest version!
package jnt.FFT;
/** Computes FFT's of real, single 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 RealFloatFFT_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 RealFloatFFT_Radix2 extends RealFloatFFT {
private int logn;
/** Create an FFT for transforming n points of real, single precision data.
* n must be an integral power of 2. */
public RealFloatFFT_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 (float 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++) {
float t0_real = data[i0+stride*b*p] + data[i0+stride*(b*p + p_1)];
float 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 */
{
float w_real = 1.0f;
float w_imag = 0.0f;
double theta = - 2.0 * Math.PI / p;
float s = (float) Math.sin(theta);
float t = (float) Math.sin(theta / 2.0);
float s2 = 2.0f * t * t;
for (a = 1; a < (p_1)/2; a++) {
/* trignometric recurrence for w-> exp(i theta) w */
{
float tmp_real = w_real - s * w_imag - s2 * w_real;
float tmp_imag = w_imag + s * w_real - s2 * w_imag;
w_real = tmp_real;
w_imag = tmp_imag;
}
for (b = 0; b < q; b++) {
float z0_real = data[i0+stride*(b*p + a)];
float z0_imag = data[i0+stride*(b*p + p_1 - a)];
float z1_real = data[i0+stride*(b*p + p_1 + a)];
float 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 (float 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++) {
float z0 = data[i0+stride*b*p];
float 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 */
{
float w_real = 1.0f;
float w_imag = 0.0f;
double theta = 2.0 * Math.PI / p;
float s = (float) Math.sin(theta);
float t = (float) Math.sin(theta / 2.0);
float s2 = 2.0f * t * t;
for (a = 1; a < (p_1)/2; a++) {
/* trignometric recurrence for w-> exp(i theta) w */
float tmp_real = w_real - s * w_imag - s2 * w_real;
float tmp_imag = w_imag + s * w_real - s2 * w_imag;
w_real = tmp_real;
w_imag = tmp_imag;
for (b = 0; b < q; b++) {
float z0_real = data[i0+stride*(b*p + a)];
float z0_imag = data[i0+stride*(b*p + p - a)];
float z1_real = data[i0+stride*(b*p + p_1 - a)];
float 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) */
float t1_real = z0_real - z1_real;
float 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 float[] toWraparoundOrder(float 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 float[] toWraparoundOrder(float data[], int i0, int stride) {
checkData(data,i0,stride);
float newdata[] = new float[2*n];
int nh = n/2;
newdata[0] = data[i0];
newdata[1] = 0.0f;
newdata[n] = data[i0+stride*nh];
newdata[n+1] = 0.0f;
for(int i=1; i