jm.audio.math.RealFloatFFT_Radix2 Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of jmusic Show documentation
Show all versions of jmusic Show documentation
JMusic - Java Music Library
The newest version!
package jm.audio.math;
/**
* 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 < nh; i++) {
newdata[2 * i] = data[i0 + stride * i];
newdata[2 * i + 1] = data[i0 + stride * (n - i)];
newdata[2 * (n - i)] = data[i0 + stride * i];
newdata[2 * (n - i) + 1] = -data[i0 + stride * (n - i)];
}
return newdata;
}
protected void bitreverse(float data[], int i0, int stride) {
/* This is the Goldrader bit-reversal algorithm */
for (int i = 0, j = 0; i < n - 1; i++) {
int k = n / 2;
if (i < j) {
float tmp = data[i0 + stride * i];
data[i0 + stride * i] = data[i0 + stride * j];
data[i0 + stride * j] = tmp;
}
while (k <= j) {
j = j - k;
k = k / 2;
}
j += k;
}
}
}