archive-gh-me/MeteoInfo
1
0
Fork 0
mirror of https://github.com/meteoinfo/MeteoInfo.git synced 2025-12-08 20:36:05 +00:00
Code Issues Projects Releases Wiki Activity

fft support Bluestein chirp-z transform algorithm

Browse Source
This commit is contained in:
wyq 2024-10-11 10:20:58 +08:00
parent b8fac5eb11
commit 76612e03da
7 changed files with 267 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
meteoinfo-common/src/main/java/org/meteoinfo/common/util/GlobalUtil.java
View File

@ -67,7 +67,7 @@ import java.util.zip.ZipInputStream;
public static String getVersion() {
String version = GlobalUtil.class.getPackage().getImplementationVersion();
if (version == null || version.equals("")) {
version = "3.9.5";
version = "3.9.6";
}
return version;
}

10
meteoinfo-lab/milconfig.xml
View File

@ -21,16 +21,16 @@
<OpenedFiles>
<OpenedFile File="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\slice\slice_2d.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\geoshow\layer_1.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\common_math\fft\fft_3.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\common_math\fft\fft_3-1.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\common_math\fft\fft_5.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\common_math\fft\fft2_1.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\common_math\fft\ifft2_1.py"/>
</OpenedFiles>
<RecentFiles>
<RecentFile File="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\slice\slice_2d.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\geoshow\layer_1.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\common_math\fft\fft_3.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\common_math\fft\fft_3-1.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\common_math\fft\fft_5.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\common_math\fft\fft2_1.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\common_math\fft\ifft2_1.py"/>
</RecentFiles>
</File>
<Font>
@ -38,5 +38,5 @@
</Font>
<LookFeel DockWindowDecorated="true" LafDecorated="true" Name="FlatDarkLaf"/>
<Figure DoubleBuffering="true"/>
<Startup MainFormLocation="-7,0" MainFormSize="1373,822"/>
<Startup MainFormLocation="-7,0" MainFormSize="1327,827"/>
</MeteoInfo>

163
meteoinfo-math/src/main/java/org/meteoinfo/math/transform/FastFourierTransform.java
View File

@ -2,9 +2,15 @@ package org.meteoinfo.math.transform;
import org.apache.commons.numbers.core.ArithmeticUtils;
import org.meteoinfo.ndarray.Array;
import org.meteoinfo.ndarray.ArrayComplex;
import org.meteoinfo.ndarray.Complex;
import org.meteoinfo.ndarray.DataType;
import org.meteoinfo.ndarray.math.ArrayMath;
import org.meteoinfo.ndarray.math.ArrayUtil;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
import java.util.function.DoubleUnaryOperator;
/**
@ -28,6 +34,9 @@ public class FastFourierTransform implements ComplexTransform {
/** Number of array slots: 1 for "real" parts 1 for "imaginary" parts. */
private static final int NUM_PARTS = 2;
/** The IEEE 754 machine epsilon from Cephes: {@code (2^-53)} */
private static final double MACH_EPS = 1.11022302462515654042e-16;
static final double TOL = 5.0 * MACH_EPS;
/**
* {@code W_SUB_N_R[i]} is the real part of
* {@code exp(- 2 * i * pi / n)}:
@ -318,14 +327,151 @@ public class FastFourierTransform implements ComplexTransform {
@Override
public Array apply(Array f) {
if (!ArithmeticUtils.isPowerOfTwo(f.getSize())) {
f = pad(f);
//f = pad(f);
if (inverse) {
return ifftBluestein(f);
} else {
return fftBluestein(f);
}
}
final double[][] dataRI = TransformUtils.createRealImaginary(f);
transformInPlace(dataRI);
tolTransformedData(dataRI);
return TransformUtils.createComplexArray(dataRI);
}
/**
* Fast Fourier Transform using Bluestein method
* @param a Input data array
* @return FFT data array
*/
public Array fftBluestein(Array a) {
// find a power of 2 convolution length m such that m >= n * 2 + 1
int n = (int) a.getSize();
if (n >= 0x20000000) {
throw new IllegalArgumentException("array too large: " + n);
}
double[] cos = new double[n];
double[] sin = new double[n];
for (int i = 0; i < n; ++i) {
int j = (int) ((long) i * i % (n * 2));
double angle = Math.PI * j / n;
cos[i] = Math.cos(angle);
sin[i] = Math.sin(angle);
}
int m = Integer.highestOneBit(n) * 4;
// temporary arrays
double[] a_re = new double[m];
double[] a_im = new double[m];
double[] b_re = new double[m];
double[] b_im = new double[m];
b_re[0] = cos[0];
b_im[0] = sin[0];
Complex c;
for (int i = 0; i < n; ++i) {
double sin_i = sin[i];
double cos_i = cos[i];
c = a.getComplex(i);
double re_i = c.real();
double im_i = c.imag();
a_re[i] = re_i * cos_i + im_i * sin_i;
a_im[i] = -re_i * sin_i + im_i * cos_i;
if (i != 0) {
b_re[i] = b_re[m - i] = cos_i;
b_im[i] = b_im[m - i] = sin_i;
}
}
// convolution
Array conv = convolve(ArrayUtil.arrayComplex(a_re, a_im), ArrayUtil.arrayComplex(b_re, b_im));
List<double[]> re_im = ((ArrayComplex) conv).getRealImage();
double[] c_re = re_im.get(0);
double[] c_im = re_im.get(1);
// result
double[] re = new double[n];
double[] im = new double[n];
// postprocessing
Array r = Array.factory(DataType.COMPLEX, new int[]{n});
for (int i = 0; i < n; ++i) {
double sin_i = sin[i];
double cos_i = cos[i];
double c_re_i = c_re[i];
double c_im_i = c_im[i];
double re_i = c_re_i * cos_i + c_im_i * sin_i;
double im_i = -c_re_i * sin_i + c_im_i * cos_i;
c = new Complex((Math.abs(re_i) <= TOL) ? 0.0 : re_i,
(Math.abs(im_i) <= TOL) ? 0.0 : im_i);
r.setComplex(i, c);
}
return r;
}
/**
* Fast Fourier Inverse Transform using Bluestein method
* @param a Input data array
* @return Inverse FFT data array
*/
public Array ifftBluestein(Array freqs) {
Array inv = fftBluestein(freqs);
List<double[]> re_im = ((ArrayComplex) inv).getRealImage();
double[] re = re_im.get(0);
double[] im = re_im.get(1);
final int n = re.length;
for (int i = 0; i < n; ++i) {
double re_i = re[i] / n;
double im_i = im[i] / n;
re[i] = (Math.abs(re_i) <= TOL) ? 0.0 : re_i;
im[i] = (Math.abs(im_i) <= TOL) ? 0.0 : im_i;
}
for (int i = 1; i <= n / 2; ++i) {
double re_tmp = re[n - i];
double im_tmp = im[n - i];
re[n - i] = re[i];
re[i] = re_tmp;
im[n - i] = im[i];
im[i] = im_tmp;
}
//return inv;
return new ArrayComplex(re, im);
}
private Array convolve(Array x, Array y) {
FastFourierTransform fft = new FastFourierTransform(this.normalization, false);
x = fft.apply(x);
y = fft.apply(y);
Array x_re = ArrayMath.getReal(x);
Array x_im = ArrayMath.getImage(x);
Array y_re = ArrayMath.getReal(y);
Array y_im = ArrayMath.getImage(y);
Array temp = Array.factory(DataType.COMPLEX, x.getShape());
Complex c;
for (int i = 0; i < x_re.getSize(); ++i) {
double x_re_i = x_re.getDouble(i);
double y_re_i = y_re.getDouble(i);
double x_im_i = x_im.getDouble(i);
double y_im_i = y_im.getDouble(i);
c = new Complex(x_re_i * y_re_i - x_im_i * y_im_i,
x_im_i * y_re_i + x_re_i * y_im_i);
temp.setComplex(i, c);
}
fft = new FastFourierTransform(this.normalization, true);
return fft.apply(temp);
}
/**
* Applies normalization to the transformed data.
*
@ -366,6 +512,21 @@ public class FastFourierTransform implements ComplexTransform {
}
}
private void tolTransformedData(final double[][] dataRI) {
final double[] dataR = dataRI[0];
final double[] dataI = dataRI[1];
final int n = dataR.length;
for (int i = 0; i < n; i++) {
if (Math.abs(dataR[i]) <= TOL) {
dataR[i] = 0;
}
if (Math.abs(dataI[i]) <= TOL) {
dataI[i] = 0;
}
}
}
/**
* Performs identical index bit reversal shuffles on two arrays of
* identical size.

11
meteoinfo-math/src/main/java/org/meteoinfo/math/transform/FastFourierTransform2D.java
View File

@ -12,6 +12,13 @@ import java.util.List;
public class FastFourierTransform2D extends FastFourierTransform {
/**
* Constructor
*/
public FastFourierTransform2D() {
super();
}
/**
* Constructor
* @param inverse Whether is inverse transform
@ -60,8 +67,4 @@ public class FastFourierTransform2D extends FastFourierTransform {
return r;
}
@Override
public Array apply(double[] f) {
return null;
}
}

5
meteoinfo-math/src/main/java/org/meteoinfo/math/transform/Fourier.java Normal file
View File

@ -0,0 +1,5 @@
package org.meteoinfo.math.transform;
public class Fourier {
}

71
meteoinfo-ndarray/src/main/java/org/meteoinfo/ndarray/ArrayComplex.java
View File

@ -35,6 +35,8 @@ package org.meteoinfo.ndarray;
import java.nio.ByteBuffer;
import java.nio.DoubleBuffer;
import java.time.LocalDateTime;
import java.util.Arrays;
import java.util.List;
/**
* Concrete implementation of Array specialized for complex. Data storage is
@ -95,6 +97,19 @@ public class ArrayComplex extends Array {
storage = new Complex[(int) indexCalc.getSize()];
}
/**
* Constructor
* @param real Real data array
* @param image Image data array
*/
public ArrayComplex(double[] real, double[] image) {
super(DataType.COMPLEX, new int[]{real.length});
storage = new Complex[real.length];
for (int i = 0; i < real.length; i++) {
storage[i] = new Complex(real[i], image[1]);
}
}
/**
* create new Array with given indexImpl and the same backing store
* @return New Array view
@ -387,6 +402,62 @@ public class ArrayComplex extends Array {
storage[index] = (Complex)value;
}
/**
* Get real data array
* @return Real data array
*/
public double[] getReal() {
double[] real = new double[(int) this.getSize()];
IndexIterator iter = this.getIndexIterator();
int i = 0;
Complex c;
while (iter.hasNext()) {
c = iter.getComplexNext();
real[i] = c.real();
i += 1;
}
return real;
}
/**
* Get image data array
* @return Image data array
*/
public double[] getImage() {
double[] image = new double[(int) this.getSize()];
IndexIterator iter = this.getIndexIterator();
int i = 0;
Complex c;
while (iter.hasNext()) {
c = iter.getComplexNext();
image[i] = c.imag();
i += 1;
}
return image;
}
/**
* Get real and image data array
* @return Real and image data array
*/
public List<double[]> getRealImage() {
double[] real = new double[(int) this.getSize()];
double[] image = new double[(int) this.getSize()];
IndexIterator iter = this.getIndexIterator();
int i = 0;
Complex c;
while (iter.hasNext()) {
c = iter.getComplexNext();
real[i] = c.real();
image[i] = c.imag();
i += 1;
}
return Arrays.asList(real, image);
}
/**
* Concrete implementation of Array specialized for doubles, rank 0.
*/

16
meteoinfo-ndarray/src/main/java/org/meteoinfo/ndarray/math/ArrayUtil.java
View File

@ -752,6 +752,22 @@ public class ArrayUtil {
}
}
/**
* Create complex array from real and image data array
* @param real Real data array
* @param image Image data array
* @return Complex array
*/
public static Array arrayComplex(double[] real, double[] image) {
int n = real.length;
Array r = Array.factory(DataType.COMPLEX, new int[]{n});
for (int i = 0; i < n; i++) {
r.setComplex(i, new Complex(real[i], image[i]));
}
return r;
}
/**
* Create an array
*