diff --git a/meteoinfo-chart/src/main/java/org/meteoinfo/chart/graphic/TriMeshGraphic.java b/meteoinfo-chart/src/main/java/org/meteoinfo/chart/graphic/TriMeshGraphic.java index 7b48e8f8..de3c11ad 100644 --- a/meteoinfo-chart/src/main/java/org/meteoinfo/chart/graphic/TriMeshGraphic.java +++ b/meteoinfo-chart/src/main/java/org/meteoinfo/chart/graphic/TriMeshGraphic.java @@ -2,7 +2,10 @@ package org.meteoinfo.chart.graphic; import org.joml.Vector3f; import org.meteoinfo.chart.jogl.Transform; +import org.meteoinfo.common.Extent; import org.meteoinfo.common.Extent3D; +import org.meteoinfo.geometry.graphic.Graphic; +import org.meteoinfo.geometry.legend.ColorBreak; import org.meteoinfo.geometry.legend.LegendManage; import org.meteoinfo.geometry.colors.TransferFunction; import org.meteoinfo.geometry.graphic.GraphicCollection3D; @@ -17,7 +20,7 @@ import java.awt.*; import java.util.*; import java.util.List; -public class TriMeshGraphic extends GraphicCollection3D { +public class TriMeshGraphic extends Graphic { protected Logger logger = LoggerFactory.getLogger("TriMeshGraphic"); protected float[] vertexPosition; @@ -30,6 +33,10 @@ public class TriMeshGraphic extends GraphicCollection3D { protected boolean edgeInterp; protected boolean mesh; protected boolean normalLoaded = false; + protected Extent extent; + protected LegendScheme legendScheme; + protected boolean singleLegend = true; + protected ColorBreak legendBreak; /** * Constructor @@ -239,6 +246,46 @@ public class TriMeshGraphic extends GraphicCollection3D { updateExtent(); } + /** + * Set triangles + * @param vertexData The triangle vertex array + */ + public void setTriangles(List triangles) { + LinkedHashMap map = new LinkedHashMap(); + int n = triangles.size(); + this.vertexIndices = new int[n]; + Vector3f vector3f; + int idx = 0, vertexIdx = 0, triangleIdx = 0, index, ii; + List idxList = new ArrayList<>(); + for (int i = 0; i < n; i++) { + Triangle3D triangle = triangles.get(i); + for (int j = 0; j < 3; j++) { + vector3f = new Vector3f(triangle.getPoint(j)); + if (map.containsKey(vector3f)) { + index = map.get(vector3f); + vertexIndices[vertexIdx] = index; + } else { + vertexIndices[vertexIdx] = idx; + map.put(vector3f, idx++); + } + + vertexIdx += 1; + } + triangleIdx += 1; + } + + this.vertexPosition = new float[map.size() * 3]; + idx = 0; + for (Map.Entry entry : map.entrySet()) { + vector3f = entry.getKey(); + vertexPosition[idx++] = vector3f.x; + vertexPosition[idx++] = vector3f.y; + vertexPosition[idx++] = vector3f.z; + } + + updateExtent(); + } + /** * Set triangles * @param vertexData The triangle vertex @@ -430,12 +477,77 @@ public class TriMeshGraphic extends GraphicCollection3D { } } + /** + * Get extent + * + * @return The extent + */ @Override + public Extent getExtent() { + return extent; + } + + /** + * Set extent + * + * @param value Extent + */ + @Override + public void setExtent(Extent value) { + this.extent = value; + } + + /** + * Get legend scheme + * @return Legend scheme + */ + public LegendScheme getLegendScheme() { + return this.legendScheme; + } + + /** + * Set legend scheme + * @param ls Legend scheme + */ public void setLegendScheme(LegendScheme ls) { - super.setLegendScheme(ls); + this.legendScheme = ls; updateVertexColor(); } + /** + * Get is single legend or not + * @return Boolean + */ + public boolean isSingleLegend() { + return this.singleLegend; + } + + /** + * Set single legend or not + * @param value Boolean + */ + public void setSingleLegend(boolean value) { + this.singleLegend = value; + } + + /** + * Get legend break + * + * @return Legend break + */ + public ColorBreak getLegendBreak() { + return this.legendBreak; + } + + /** + * Set legend break + * + * @param value Legend break + */ + public void setLegendBreak(ColorBreak value) { + this.legendBreak = value; + } + /** * Set transfer function * @param transferFunction Transfer function diff --git a/meteoinfo-chart/src/main/java/org/meteoinfo/chart/graphic/Triangle3D.java b/meteoinfo-chart/src/main/java/org/meteoinfo/chart/graphic/Triangle3D.java index 291b11da..12e665d6 100644 --- a/meteoinfo-chart/src/main/java/org/meteoinfo/chart/graphic/Triangle3D.java +++ b/meteoinfo-chart/src/main/java/org/meteoinfo/chart/graphic/Triangle3D.java @@ -26,6 +26,24 @@ public class Triangle3D { this.pointC = c; } + /** + * Get point by index + * @param index The index + * @return The point + */ + public Vector3f getPoint(int index) { + switch (index) { + case 0: + return pointA; + case 1: + return pointB; + case 2: + return pointC; + default: + return null; + } + } + /** * Get point a * @return Point a diff --git a/meteoinfo-chart/src/main/java/org/meteoinfo/chart/jogl/GLPlot.java b/meteoinfo-chart/src/main/java/org/meteoinfo/chart/jogl/GLPlot.java index 7a72144f..e509beb5 100644 --- a/meteoinfo-chart/src/main/java/org/meteoinfo/chart/jogl/GLPlot.java +++ b/meteoinfo-chart/src/main/java/org/meteoinfo/chart/jogl/GLPlot.java @@ -2676,6 +2676,18 @@ public class GLPlot extends Plot { } } + protected void drawTriMeshGraphic(GL2 gl, TriMeshGraphic graphic) { + if (!this.renderMap.containsKey(graphic)) { + renderMap.put(graphic, new TriMeshRender(gl, graphic)); + } + TriMeshRender triMeshRender = (TriMeshRender) renderMap.get(graphic); + triMeshRender.setTransform(this.transform, this.alwaysUpdateBuffers); + triMeshRender.setOrthographic(this.orthographic); + triMeshRender.setLighting(this.lighting); + triMeshRender.updateMatrix(); + triMeshRender.draw(); + } + protected void drawGraphics(GL2 gl, Graphic graphic) { boolean lightEnabled = this.lighting.isEnable(); if (graphic instanceof GraphicCollection3D) { @@ -2718,15 +2730,7 @@ public class GLPlot extends Plot { modelRender.setRotateModelView(this.modelViewMatrixR); modelRender.draw(); } else { - if (!this.renderMap.containsKey(graphic)) { - renderMap.put(graphic, new TriMeshRender(gl, (TriMeshGraphic) graphic)); - } - TriMeshRender triMeshRender = (TriMeshRender) renderMap.get(graphic); - triMeshRender.setTransform(this.transform, this.alwaysUpdateBuffers); - triMeshRender.setOrthographic(this.orthographic); - triMeshRender.setLighting(this.lighting); - triMeshRender.updateMatrix(); - triMeshRender.draw(); + drawTriMeshGraphic(gl, (TriMeshGraphic) graphic); } } else if (graphic instanceof VolumeGraphic) { try { @@ -2762,7 +2766,7 @@ public class GLPlot extends Plot { } } if (isDraw) { - switch (graphic.getGraphicN(0).getShape().getShapeType()) { + switch (graphic.getGraphicN(0).getShapeType()) { case POINT_Z: if (!this.renderMap.containsKey(graphic)) { renderMap.put(graphic, new PointRender(gl, (GraphicCollection3D) graphic)); @@ -2815,7 +2819,11 @@ public class GLPlot extends Plot { default: for (int i = 0; i < graphic.getNumGraphics(); i++) { Graphic gg = graphic.getGraphicN(i); - this.drawGraphic(gl, gg); + if (gg instanceof TriMeshGraphic) { + this.drawTriMeshGraphic(gl, (TriMeshGraphic) gg); + } else { + this.drawGraphic(gl, gg); + } } break; } @@ -3339,83 +3347,6 @@ public class GLPlot extends Plot { } } - private void drawPolygon_bak(GL2 gl, PolygonZ aPG, PolygonBreak aPGB) { - PointZ p; - if (aPGB.isDrawFill() && aPGB.getColor().getAlpha() > 0) { - gl.glEnable(GL2.GL_POLYGON_OFFSET_FILL); - gl.glPolygonOffset(1.0f, 1.0f); - - float[] rgba = aPGB.getColor().getRGBComponents(null); - gl.glColor4f(rgba[0], rgba[1], rgba[2], rgba[3]); - - try { - GLUtessellator tobj = glu.gluNewTess(); - //TessCallback tessCallback = new TessCallback(gl, glu); - - glu.gluTessCallback(tobj, GLU.GLU_TESS_VERTEX, tessCallback); - glu.gluTessCallback(tobj, GLU.GLU_TESS_BEGIN, tessCallback); - glu.gluTessCallback(tobj, GLU.GLU_TESS_END, tessCallback); - glu.gluTessCallback(tobj, GLU.GLU_TESS_ERROR, tessCallback); - //glu.gluTessCallback(tobj, GLU.GLU_TESS_COMBINE, tessCallback); - - //gl.glNewList(startList, GL2.GL_COMPILE); - //gl.glShadeModel(GL2.GL_FLAT); - glu.gluTessBeginPolygon(tobj, null); - glu.gluTessBeginContour(tobj); - double[] v; - for (int i = 0; i < aPG.getOutLine().size() - 1; i++) { - p = ((java.util.List) aPG.getOutLine()).get(i); - v = transform.transform(p); - glu.gluTessVertex(tobj, v, 0, v); - } - glu.gluTessEndContour(tobj); - if (aPG.hasHole()) { - for (int i = 0; i < aPG.getHoleLineNumber(); i++) { - glu.gluTessBeginContour(tobj); - for (int j = 0; j < aPG.getHoleLine(i).size() - 1; j++) { - p = ((java.util.List) aPG.getHoleLine(i)).get(j); - v = transform.transform(p); - glu.gluTessVertex(tobj, v, 0, v); - } - glu.gluTessEndContour(tobj); - } - } - glu.gluTessEndPolygon(tobj); - //gl.glEndList(); - glu.gluDeleteTess(tobj); - - //gl.glCallList(startList); - } catch (Exception e) { - e.printStackTrace(); - } - } - if (aPGB.isDrawOutline()) { - float[] rgba = aPGB.getOutlineColor().getRGBComponents(null); - gl.glLineWidth(aPGB.getOutlineSize() * this.dpiScale); - gl.glColor4f(rgba[0], rgba[1], rgba[2], rgba[3]); - gl.glBegin(GL2.GL_LINE_STRIP); - for (int i = 0; i < aPG.getOutLine().size(); i++) { - p = ((java.util.List) aPG.getOutLine()).get(i); - gl.glVertex3f(transform.transform_x((float) p.X), transform.transform_y((float) p.Y), transform.transform_z((float) p.Z)); - } - gl.glEnd(); - - if (aPG.hasHole()) { - java.util.List newPList; - gl.glBegin(GL2.GL_LINE_STRIP); - for (int h = 0; h < aPG.getHoleLines().size(); h++) { - newPList = (java.util.List) aPG.getHoleLines().get(h); - for (int j = 0; j < newPList.size(); j++) { - p = newPList.get(j); - gl.glVertex3f(transform.transform_x((float) p.X), transform.transform_y((float) p.Y), transform.transform_z((float) p.Z)); - } - } - gl.glEnd(); - } - gl.glDisable(GL2.GL_POLYGON_OFFSET_FILL); - } - } - private void drawConvexPolygon(GL2 gl, PolygonZ aPG, PolygonBreak aPGB) { PointZ p; if (aPGB.isDrawFill()) { diff --git a/meteoinfo-chart/src/main/java/org/meteoinfo/chart/jogl/Triangle.java b/meteoinfo-chart/src/main/java/org/meteoinfo/chart/jogl/Triangle.java deleted file mode 100644 index da0a38c7..00000000 --- a/meteoinfo-chart/src/main/java/org/meteoinfo/chart/jogl/Triangle.java +++ /dev/null @@ -1,38 +0,0 @@ -package org.meteoinfo.chart.jogl; - -import org.meteoinfo.geometry.shape.PointZ; - -import java.util.ArrayList; -import java.util.List; - -public class Triangle { - private static final long serialVersionUID = 1; - - public final PointZ p1; - public final PointZ p2; - public final PointZ p3; - - public Triangle(PointZ p1, PointZ p2, PointZ p3) { - super(); - this.p1 = p1; - this.p2 = p2; - this.p3 = p3; - } - - public List getPoints() { - List points = new ArrayList<>(); - points.add(p1); - points.add(p2); - points.add(p3); - return points; - } - - public PointZ[] getPointArray() { - return new PointZ[]{p1, p2, p3}; - } - - @Override - public String toString() { - return "Triangle [p1=" + p1 + ", p2=" + p2 + ", p3=" + p3 + "]"; - } -} diff --git a/meteoinfo-chart/src/main/java/org/meteoinfo/chart/jogl/tessellator/TriangleTessellator.java b/meteoinfo-chart/src/main/java/org/meteoinfo/chart/jogl/tessellator/TriangleTessellator.java index 4f3b9732..7f1f6d36 100644 --- a/meteoinfo-chart/src/main/java/org/meteoinfo/chart/jogl/tessellator/TriangleTessellator.java +++ b/meteoinfo-chart/src/main/java/org/meteoinfo/chart/jogl/tessellator/TriangleTessellator.java @@ -4,9 +4,10 @@ import com.jogamp.opengl.GL2; import com.jogamp.opengl.glu.GLU; import com.jogamp.opengl.glu.GLUtessellator; import com.jogamp.opengl.glu.GLUtessellatorCallbackAdapter; -import org.meteoinfo.chart.jogl.Triangle; +import org.meteoinfo.chart.graphic.Triangle3D; import org.meteoinfo.geometry.shape.PointZ; import org.meteoinfo.geometry.shape.PolygonZ; +import org.joml.Vector3f; import java.util.ArrayList; import java.util.List; @@ -37,7 +38,7 @@ public class TriangleTessellator { * Throws {@link TesselationException} if the tessellation was * unsuccessful, most commonly due to ambiguous shapes */ - public List getTriangles(PolygonZ polygon) + public List getTriangles(PolygonZ polygon) throws TesselationException { makeTriangles(polygon); @@ -83,7 +84,7 @@ public class TriangleTessellator { } public interface TessellatorListener { - public void onTesselationDone(List triangles); + public void onTesselationDone(List triangles); public void onTesselationError(TesselationException err); } @@ -163,9 +164,9 @@ public class TriangleTessellator { */ class TessellationCallback extends GLUtessellatorCallbackAdapter { - protected List triangles = new ArrayList<>(); + protected List triangles = new ArrayList<>(); - private PointZ p1, p2, p3; + private Vector3f p1, p2, p3; private int geometricPrimitiveType; @@ -178,7 +179,8 @@ public class TriangleTessellator { } public void vertex(Object vertexData) { - PointZ thisPoint = new PointZ((double[]) vertexData); + double[] vertex = (double[]) vertexData; + Vector3f thisPoint = new Vector3f((float) vertex[0], (float) vertex[1], (float) vertex[2]); switch (geometricPrimitiveType) { case GL2.GL_TRIANGLE_FAN: @@ -187,7 +189,7 @@ public class TriangleTessellator { } else if (p2 == null) { p2 = thisPoint; } else { - triangles.add(new Triangle(p1, p2, thisPoint)); + triangles.add(new Triangle3D(p1, p2, thisPoint)); p2 = thisPoint; } break; @@ -201,7 +203,7 @@ public class TriangleTessellator { p3 = thisPoint; if (p1 != null) { - triangles.add(new Triangle(p1, p2, p3)); + triangles.add(new Triangle3D(p1, p2, p3)); } } break; diff --git a/meteoinfo-lab/milconfig.xml b/meteoinfo-lab/milconfig.xml index 8afd01e1..3f89eb8b 100644 --- a/meteoinfo-lab/milconfig.xml +++ b/meteoinfo-lab/milconfig.xml @@ -1,32 +1,36 @@ - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + - + + + - + + + @@ -34,5 +38,5 @@
- + diff --git a/meteoinfo-lab/pylib/mipylib/numeric/__init__$py.class b/meteoinfo-lab/pylib/mipylib/numeric/__init__$py.class index c31ace0e..31b67b07 100644 Binary files a/meteoinfo-lab/pylib/mipylib/numeric/__init__$py.class and b/meteoinfo-lab/pylib/mipylib/numeric/__init__$py.class differ diff --git a/meteoinfo-lab/pylib/mipylib/numeric/__init__.py b/meteoinfo-lab/pylib/mipylib/numeric/__init__.py index f8520819..35e7e56f 100644 --- a/meteoinfo-lab/pylib/mipylib/numeric/__init__.py +++ b/meteoinfo-lab/pylib/mipylib/numeric/__init__.py @@ -14,9 +14,10 @@ from . import signal from . import spatial from . import special from . import integrate +from . import fft __all__ = ['linalg', 'fitting', 'random', 'ma', 'stats', 'interpolate', 'optimize', 'signal', 'spatial', - 'special', 'integrate'] + 'special', 'integrate', 'fft'] __all__.extend(['__version__']) __all__.extend(core.__all__) __all__.extend(lib.__all__) diff --git a/meteoinfo-lab/pylib/mipylib/numeric/core/numeric$py.class b/meteoinfo-lab/pylib/mipylib/numeric/core/numeric$py.class index b77d1c3a..ea4a4d35 100644 Binary files a/meteoinfo-lab/pylib/mipylib/numeric/core/numeric$py.class and b/meteoinfo-lab/pylib/mipylib/numeric/core/numeric$py.class differ diff --git a/meteoinfo-lab/pylib/mipylib/numeric/core/numeric.py b/meteoinfo-lab/pylib/mipylib/numeric/core/numeric.py index 5bccf7ef..45d09883 100644 --- a/meteoinfo-lab/pylib/mipylib/numeric/core/numeric.py +++ b/meteoinfo-lab/pylib/mipylib/numeric/core/numeric.py @@ -101,6 +101,19 @@ def array(object, dtype=None, copy=True, order='K', subok=False, ndmin=0): elif isinstance(object[0], (list, tuple)): if miutil.iscomplex(object[0]): a = NDArray(JythonUtil.toComplexArray(object)) + elif isinstance(object[0], NDArray): + shape = [len(object)] + shape.extend(object[0].shape) + if dtype is None: + for o in object: + if dtype is None: + dtype = o.dtype + else: + if o.dtype > dtype: + dtype = o.dtype + a = zeros(shape, dtype=dtype) + for i in range(len(object)): + a[i] = object[i] elif miutil.iscomplex(object): a = NDArray(JythonUtil.toComplexArray(object)) diff --git a/meteoinfo-lab/pylib/mipylib/numeric/fft/__init__.py b/meteoinfo-lab/pylib/mipylib/numeric/fft/__init__.py new file mode 100644 index 00000000..897f8c75 --- /dev/null +++ b/meteoinfo-lab/pylib/mipylib/numeric/fft/__init__.py @@ -0,0 +1,2 @@ +from ._fft import * +from ._helper import * \ No newline at end of file diff --git a/meteoinfo-lab/pylib/mipylib/numeric/fft/_fft.py b/meteoinfo-lab/pylib/mipylib/numeric/fft/_fft.py new file mode 100644 index 00000000..eba2ed8a --- /dev/null +++ b/meteoinfo-lab/pylib/mipylib/numeric/fft/_fft.py @@ -0,0 +1,127 @@ +from org.meteoinfo.math.transform import FastFourierTransform, FastFourierTransform2D + +from .. import core as _nx + +__all__ = ['fft', 'ifft', 'fft2', 'ifft2'] + + +def fft(a): + """ + Compute the one-dimensional discrete Fourier Transform. + + This function computes the one-dimensional *n*-point discrete Fourier + Transform (DFT) with the efficient Fast Fourier Transform (FFT) + algorithm [CT]. + + Parameters + ---------- + a : array_like + Input array, can be complex. + + Returns + ------- + out : complex ndarray + The truncated or zero-padded input, transformed along the axis + indicated by `axis`, or the last one if `axis` is not specified. + + References + ---------- + .. [CT] Cooley, James W., and John W. Tukey, 1965, "An algorithm for the + machine calculation of complex Fourier series," *Math. Comput.* + 19: 297-301. + """ + a = _nx.asarray(a) + jfft = FastFourierTransform() + r = jfft.apply(a._array) + + return _nx.NDArray(r) + +def ifft(a): + """ + Compute the one-dimensional inverse discrete Fourier Transform. + + This function computes the inverse of the one-dimensional *n*-point + discrete Fourier transform computed by `fft`. In other words, + ``ifft(fft(a)) == a`` to within numerical accuracy. + For a general description of the algorithm and definitions, + see `numpy.fft`. + + The input should be ordered in the same way as is returned by `fft`, + i.e., + + Parameters + ---------- + a : array_like + Input array, can be complex. + + Returns + ------- + out : complex ndarray + The truncated or zero-padded input, transformed along the axis + indicated by `axis`, or the last one if `axis` is not specified. + """ + a = _nx.asarray(a) + jfft = FastFourierTransform(True) + r = jfft.apply(a._array) + + return _nx.NDArray(r) + +def fft2(a): + """ + Compute the 2-dimensional discrete Fourier Transform. + + This function computes the *n*-dimensional discrete Fourier Transform + over any axes in an *M*-dimensional array by means of the + Fast Fourier Transform (FFT). By default, the transform is computed over + the last two axes of the input array, i.e., a 2-dimensional FFT. + + Parameters + ---------- + a : array_like + Input array, can be complex + + Returns + ------- + out : complex ndarray + The truncated or zero-padded input, transformed along the axes + indicated by `axes`, or the last two axes if `axes` is not given. + """ + a = _nx.asarray(a) + jfft = FastFourierTransform2D() + r = jfft.apply(a._array) + + return _nx.NDArray(r) + +def ifft2(a, s=None, axes=(-2, -1), norm=None, out=None): + """ + Compute the 2-dimensional inverse discrete Fourier Transform. + + This function computes the inverse of the 2-dimensional discrete Fourier + Transform over any number of axes in an M-dimensional array by means of + the Fast Fourier Transform (FFT). In other words, ``ifft2(fft2(a)) == a`` + to within numerical accuracy. By default, the inverse transform is + computed over the last two axes of the input array. + + The input, analogously to `ifft`, should be ordered in the same way as is + returned by `fft2`, i.e. it should have the term for zero frequency + in the low-order corner of the two axes, the positive frequency terms in + the first half of these axes, the term for the Nyquist frequency in the + middle of the axes and the negative frequency terms in the second half of + both axes, in order of decreasingly negative frequency. + + Parameters + ---------- + a : array_like + Input array, can be complex. + + Returns + ------- + out : complex ndarray + The truncated or zero-padded input, transformed along the axes + indicated by `axes`, or the last two axes if `axes` is not given. + """ + a = _nx.asarray(a) + jfft = FastFourierTransform2D(True) + r = jfft.apply(a._array) + + return _nx.NDArray(r) diff --git a/meteoinfo-lab/pylib/mipylib/numeric/fft/_helper.py b/meteoinfo-lab/pylib/mipylib/numeric/fft/_helper.py new file mode 100644 index 00000000..271acdea --- /dev/null +++ b/meteoinfo-lab/pylib/mipylib/numeric/fft/_helper.py @@ -0,0 +1,216 @@ +""" +Discrete Fourier Transforms - _helper.py + +""" +from ..core import empty, arange, asarray, roll + +# Created by Pearu Peterson, September 2002 + +__all__ = ['fftshift', 'ifftshift', 'fftfreq', 'rfftfreq'] + +integer_types = (int) + + +def _fftshift_dispatcher(x, axes=None): + return (x,) + + +def fftshift(x, axes=None): + """ + Shift the zero-frequency component to the center of the spectrum. + + This function swaps half-spaces for all axes listed (defaults to all). + Note that ``y[0]`` is the Nyquist component only if ``len(x)`` is even. + + Parameters + ---------- + x : array_like + Input array. + axes : int or shape tuple, optional + Axes over which to shift. Default is None, which shifts all axes. + + Returns + ------- + y : ndarray + The shifted array. + + See Also + -------- + ifftshift : The inverse of `fftshift`. + + Examples + -------- + >>> freqs = np.fft.fftfreq(10, 0.1) + >>> freqs + array([ 0., 1., 2., ..., -3., -2., -1.]) + >>> np.fft.fftshift(freqs) + array([-5., -4., -3., -2., -1., 0., 1., 2., 3., 4.]) + + Shift the zero-frequency component only along the second axis: + + >>> freqs = np.fft.fftfreq(9, d=1./9).reshape(3, 3) + >>> freqs + array([[ 0., 1., 2.], + [ 3., 4., -4.], + [-3., -2., -1.]]) + >>> np.fft.fftshift(freqs, axes=(1,)) + array([[ 2., 0., 1.], + [-4., 3., 4.], + [-1., -3., -2.]]) + + """ + x = asarray(x) + if axes is None: + axes = tuple(range(x.ndim)) + shift = [dim // 2 for dim in x.shape] + elif isinstance(axes, integer_types): + shift = x.shape[axes] // 2 + else: + shift = [x.shape[ax] // 2 for ax in axes] + + return roll(x, shift, axes) + + +def ifftshift(x, axes=None): + """ + The inverse of `fftshift`. Although identical for even-length `x`, the + functions differ by one sample for odd-length `x`. + + Parameters + ---------- + x : array_like + Input array. + axes : int or shape tuple, optional + Axes over which to calculate. Defaults to None, which shifts all axes. + + Returns + ------- + y : ndarray + The shifted array. + + See Also + -------- + fftshift : Shift zero-frequency component to the center of the spectrum. + + Examples + -------- + >>> freqs = np.fft.fftfreq(9, d=1./9).reshape(3, 3) + >>> freqs + array([[ 0., 1., 2.], + [ 3., 4., -4.], + [-3., -2., -1.]]) + >>> np.fft.ifftshift(np.fft.fftshift(freqs)) + array([[ 0., 1., 2.], + [ 3., 4., -4.], + [-3., -2., -1.]]) + + """ + x = asarray(x) + if axes is None: + axes = tuple(range(x.ndim)) + shift = [-(dim // 2) for dim in x.shape] + elif isinstance(axes, integer_types): + shift = -(x.shape[axes] // 2) + else: + shift = [-(x.shape[ax] // 2) for ax in axes] + + return roll(x, shift, axes) + + +def fftfreq(n, d=1.0): + """ + Return the Discrete Fourier Transform sample frequencies. + + The returned float array `f` contains the frequency bin centers in cycles + per unit of the sample spacing (with zero at the start). For instance, if + the sample spacing is in seconds, then the frequency unit is cycles/second. + + Given a window length `n` and a sample spacing `d`:: + + f = [0, 1, ..., n/2-1, -n/2, ..., -1] / (d*n) if n is even + f = [0, 1, ..., (n-1)/2, -(n-1)/2, ..., -1] / (d*n) if n is odd + + Parameters + ---------- + n : int + Window length. + d : scalar, optional + Sample spacing (inverse of the sampling rate). Defaults to 1. + + Returns + ------- + f : ndarray + Array of length `n` containing the sample frequencies. + + Examples + -------- + >>> signal = np.array([-2, 8, 6, 4, 1, 0, 3, 5], dtype=float) + >>> fourier = np.fft.fft(signal) + >>> n = signal.size + >>> timestep = 0.1 + >>> freq = np.fft.fftfreq(n, d=timestep) + >>> freq + array([ 0. , 1.25, 2.5 , ..., -3.75, -2.5 , -1.25]) + + """ + if not isinstance(n, integer_types): + raise ValueError("n should be an integer") + val = 1.0 / (n * d) + results = empty(n, dtype='int') + N = (n-1)//2 + 1 + p1 = arange(0, N, dtype='int') + results[:N] = p1 + p2 = arange(-(n//2), 0, dtype='int') + results[N:] = p2 + return results * val + + +def rfftfreq(n, d=1.0): + """ + Return the Discrete Fourier Transform sample frequencies + (for usage with rfft, irfft). + + The returned float array `f` contains the frequency bin centers in cycles + per unit of the sample spacing (with zero at the start). For instance, if + the sample spacing is in seconds, then the frequency unit is cycles/second. + + Given a window length `n` and a sample spacing `d`:: + + f = [0, 1, ..., n/2-1, n/2] / (d*n) if n is even + f = [0, 1, ..., (n-1)/2-1, (n-1)/2] / (d*n) if n is odd + + Unlike `fftfreq` (but like `scipy.fftpack.rfftfreq`) + the Nyquist frequency component is considered to be positive. + + Parameters + ---------- + n : int + Window length. + d : scalar, optional + Sample spacing (inverse of the sampling rate). Defaults to 1. + + Returns + ------- + f : ndarray + Array of length ``n//2 + 1`` containing the sample frequencies. + + Examples + -------- + >>> signal = np.array([-2, 8, 6, 4, 1, 0, 3, 5, -3, 4], dtype=float) + >>> fourier = np.fft.rfft(signal) + >>> n = signal.size + >>> sample_rate = 100 + >>> freq = np.fft.fftfreq(n, d=1./sample_rate) + >>> freq + array([ 0., 10., 20., ..., -30., -20., -10.]) + >>> freq = np.fft.rfftfreq(n, d=1./sample_rate) + >>> freq + array([ 0., 10., 20., 30., 40., 50.]) + + """ + if not isinstance(n, integer_types): + raise ValueError("n should be an integer") + val = 1.0/(n*d) + N = n//2 + 1 + results = arange(0, N, dtype='int') + return results * val \ No newline at end of file diff --git a/meteoinfo-lab/pylib/mipylib/plotlib/miplot$py.class b/meteoinfo-lab/pylib/mipylib/plotlib/miplot$py.class index 08ffe595..eb8d6c84 100644 Binary files a/meteoinfo-lab/pylib/mipylib/plotlib/miplot$py.class and b/meteoinfo-lab/pylib/mipylib/plotlib/miplot$py.class differ diff --git a/meteoinfo-math/src/main/java/org/meteoinfo/math/transform/ComplexTransform.java b/meteoinfo-math/src/main/java/org/meteoinfo/math/transform/ComplexTransform.java new file mode 100644 index 00000000..b2c2bc1c --- /dev/null +++ b/meteoinfo-math/src/main/java/org/meteoinfo/math/transform/ComplexTransform.java @@ -0,0 +1,51 @@ +package org.meteoinfo.math.transform; + +import org.meteoinfo.ndarray.Array; +import org.meteoinfo.ndarray.Complex; + +import java.util.function.DoubleUnaryOperator; +import java.util.function.UnaryOperator; + +/** + * {@link Complex} transform. + *

+ * Such transforms include {@link FastSineTransform sine transform}, + * {@link FastCosineTransform cosine transform} or {@link + * FastHadamardTransform Hadamard transform}. + */ +public interface ComplexTransform extends UnaryOperator { + /** + * Returns the transform of the specified data set. + * + * @param f the data array to be transformed (signal). + * @return the transformed array (spectrum). + * @throws IllegalArgumentException if the transform cannot be performed. + */ + Array apply(Array f); + + /** + * Returns the transform of the specified data set. + * + * @param f the data array to be transformed (signal). + * @return the transformed array (spectrum). + * @throws IllegalArgumentException if the transform cannot be performed. + */ + Array apply(double[] f); + + /** + * Returns the transform of the specified function. + * + * @param f Function to be sampled and transformed. + * @param min Lower bound (inclusive) of the interval. + * @param max Upper bound (exclusive) of the interval. + * @param n Number of sample points. + * @return the result. + * @throws IllegalArgumentException if the transform cannot be performed. + */ + default Array apply(DoubleUnaryOperator f, + double min, + double max, + int n) { + return apply(TransformUtils.sample(f, min, max, n)); + } +} diff --git a/meteoinfo-math/src/main/java/org/meteoinfo/math/transform/FFT.java b/meteoinfo-math/src/main/java/org/meteoinfo/math/transform/FFT.java new file mode 100644 index 00000000..bc346897 --- /dev/null +++ b/meteoinfo-math/src/main/java/org/meteoinfo/math/transform/FFT.java @@ -0,0 +1,196 @@ +package org.meteoinfo.math.transform; + +import org.meteoinfo.ndarray.Array; +import org.meteoinfo.ndarray.Complex; +import org.meteoinfo.ndarray.DataType; + +public class FFT { + + // compute the fast fourier transform of x[], assuming its length is a power of 2 + public static Complex[] fft(Complex[] x) { + int N = x.length; + + // base case + if (N == 1) return new Complex[] { x[0] }; + + // radix 2 Cooley-Tukey FFT + if (N % 2 != 0) { throw new RuntimeException("N is not a power of 2"); } + + // fft of even terms + Complex[] even = new Complex[N/2]; + for (int k = 0; k < N/2; k++) { + even[k] = x[2*k]; + } + Complex[] q = fft(even); + + // fft of odd terms + Complex[] odd = even; // reuse the array + for (int k = 0; k < N/2; k++) { + odd[k] = x[2*k + 1]; + } + Complex[] r = fft(odd); + + // combine + Complex[] y = new Complex[N]; + for (int k = 0; k < N/2; k++) { + double kth = -2 * k * Math.PI / N; + Complex wk = new Complex(Math.cos(kth), Math.sin(kth)); + y[k] = q[k].add(wk.multiply(r[k])); + y[k + N/2] = q[k].subtract(wk.multiply(r[k])); + } + return y; + } + + // compute the FFT of x[], assuming its length is a power of 2 + public static Array fft(Array x) { + x = x.copyIfView(); + + int N = (int) x.getSize(); + + // base case + if (N == 1) return x.copy(); + + // radix 2 Cooley-Tukey FFT + if (N % 2 != 0) { throw new RuntimeException("N is not a power of 2"); } + + // fft of even terms + Complex[] even = new Complex[N / 2]; + for (int k = 0; k < N / 2; k++) { + even[k] = x.getComplex(2 * k); + } + Complex[] q = fft(even); + + // fft of odd terms + Complex[] odd = even; // reuse the array + for (int k = 0; k < N/2; k++) { + odd[k] = x.getComplex(2 * k + 1); + } + Complex[] r = fft(odd); + + // combine + Array y = Array.factory(DataType.COMPLEX, new int[]{N}); + for (int k = 0; k < N / 2; k++) { + double kth = -2 * k * Math.PI / N; + Complex wk = new Complex(Math.cos(kth), Math.sin(kth)); + y.setComplex(k, q[k].add(wk.multiply(r[k]))); + y.setComplex(k + N / 2, q[k].subtract(wk.multiply(r[k]))); + } + return y; + } + + // compute the inverse FFT of x[], assuming its length is a power of 2 + public static Complex[] ifft(Complex[] x) { + int N = x.length; + Complex[] y = new Complex[N]; + + // take conjugate + for (int i = 0; i < N; i++) { + y[i] = x[i].conj(); + } + + // compute forward FFT + y = fft(y); + + // take conjugate again + for (int i = 0; i < N; i++) { + y[i] = y[i].conj(); + } + + // divide by N + for (int i = 0; i < N; i++) { + y[i] = y[i].multiply(1.0 / N); + } + + return y; + } + + // compute the inverse FFT of x[], assuming its length is a power of 2 + public static Array ifft(Array x) { + x = x.copyIfView(); + + int N = (int) x.getSize(); + Array y = Array.factory(DataType.COMPLEX, new int[]{N}); + + // take conjugate + for (int i = 0; i < N; i++) { + y.setComplex(i, x.getComplex(i).conj()); + } + + // compute forward FFT + y = fft(y); + + // take conjugate again + for (int i = 0; i < N; i++) { + y.setComplex(i, y.getComplex(i).conj()); + } + + // divide by N + for (int i = 0; i < N; i++) { + y.setComplex(i, y.getComplex(i).multiply(1.0 / N)); + } + + return y; + } + + // compute the circular convolution of x and y + public static Complex[] cConvolve(Complex[] x, Complex[] y) { + + // should probably pad x and y with 0s so that they have same length + // and are powers of 2 + if (x.length != y.length) { throw new RuntimeException("Dimensions don't agree"); } + + int N = x.length; + + // compute FFT of each sequence + Complex[] a = fft(x); + Complex[] b = fft(y); + + // point-wise multiply + Complex[] c = new Complex[N]; + for (int i = 0; i < N; i++) { + c[i] = a[i].multiply(b[i]); + } + + // compute inverse FFT + return ifft(c); + } + + + // compute the linear convolution of x and y + public static Complex[] convolve(Complex[] x, Complex[] y) { + Complex ZERO = new Complex(0, 0); + + Complex[] a = new Complex[2*x.length]; + for (int i = 0; i < x.length; i++) a[i] = x[i]; + for (int i = x.length; i < 2*x.length; i++) a[i] = ZERO; + + Complex[] b = new Complex[2*y.length]; + for (int i = 0; i < y.length; i++) b[i] = y[i]; + for (int i = y.length; i < 2*y.length; i++) b[i] = ZERO; + + return cConvolve(a, b); + } + + // display an array of Complex numbers to standard output + public static void show(Complex[] x, String title) { + System.out.println(title); + System.out.println("-------------------"); + for (int i = 0; i < x.length; i++) { + System.out.println(x[i]); + } + System.out.println(); + } + + + public static Complex[] toComplex(double[] SZ) + { + int count = SZ.length; + Complex[] C_SZ = new Complex[count]; + for (int i = 0; i < count; i++) + { + Complex d = new Complex(SZ[i],0); + C_SZ[i] = d; + } + return C_SZ; + } +} diff --git a/meteoinfo-math/src/main/java/org/meteoinfo/math/transform/FastFourierTransform.java b/meteoinfo-math/src/main/java/org/meteoinfo/math/transform/FastFourierTransform.java new file mode 100644 index 00000000..14adf64f --- /dev/null +++ b/meteoinfo-math/src/main/java/org/meteoinfo/math/transform/FastFourierTransform.java @@ -0,0 +1,462 @@ +package org.meteoinfo.math.transform; + +import org.apache.commons.numbers.core.ArithmeticUtils; +import org.meteoinfo.ndarray.Array; +import org.meteoinfo.ndarray.Complex; + +import java.util.Arrays; +import java.util.function.DoubleUnaryOperator; + +/** + * Implements the Fast Fourier Transform for transformation of one-dimensional + * real or complex data sets. For reference, see Applied Numerical Linear + * Algebra, ISBN 0898713897, chapter 6. + *

+ * There are several variants of the discrete Fourier transform, with various + * normalization conventions, which are specified by the parameter + * {@link Norm}. + *

+ * The current implementation of the discrete Fourier transform as a fast + * Fourier transform requires the length of the data set to be a power of 2. + * This greatly simplifies and speeds up the code. Users can pad the data with + * zeros to meet this requirement. There are other flavors of FFT, for + * reference, see S. Winograd, + * On computing the discrete Fourier transform, Mathematics of + * Computation, 32 (1978), 175 - 199. + */ +public class FastFourierTransform implements ComplexTransform { + + /** Number of array slots: 1 for "real" parts 1 for "imaginary" parts. */ + private static final int NUM_PARTS = 2; + /** + * {@code W_SUB_N_R[i]} is the real part of + * {@code exp(- 2 * i * pi / n)}: + * {@code W_SUB_N_R[i] = cos(2 * pi/ n)}, where {@code n = 2^i}. + */ + private static final double[] W_SUB_N_R = { + 0x1.0p0, -0x1.0p0, 0x1.1a62633145c07p-54, 0x1.6a09e667f3bcdp-1, + 0x1.d906bcf328d46p-1, 0x1.f6297cff75cbp-1, 0x1.fd88da3d12526p-1, 0x1.ff621e3796d7ep-1, + 0x1.ffd886084cd0dp-1, 0x1.fff62169b92dbp-1, 0x1.fffd8858e8a92p-1, 0x1.ffff621621d02p-1, + 0x1.ffffd88586ee6p-1, 0x1.fffff62161a34p-1, 0x1.fffffd8858675p-1, 0x1.ffffff621619cp-1, + 0x1.ffffffd885867p-1, 0x1.fffffff62161ap-1, 0x1.fffffffd88586p-1, 0x1.ffffffff62162p-1, + 0x1.ffffffffd8858p-1, 0x1.fffffffff6216p-1, 0x1.fffffffffd886p-1, 0x1.ffffffffff621p-1, + 0x1.ffffffffffd88p-1, 0x1.fffffffffff62p-1, 0x1.fffffffffffd9p-1, 0x1.ffffffffffff6p-1, + 0x1.ffffffffffffep-1, 0x1.fffffffffffffp-1, 0x1.0p0, 0x1.0p0, + 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0, + 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0, + 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0, + 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0, + 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0, + 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0, + 0x1.0p0, 0x1.0p0, 0x1.0p0, 0x1.0p0, + 0x1.0p0, 0x1.0p0, 0x1.0p0 }; + + /** + * {@code W_SUB_N_I[i]} is the imaginary part of + * {@code exp(- 2 * i * pi / n)}: + * {@code W_SUB_N_I[i] = -sin(2 * pi/ n)}, where {@code n = 2^i}. + */ + private static final double[] W_SUB_N_I = { + 0x1.1a62633145c07p-52, -0x1.1a62633145c07p-53, -0x1.0p0, -0x1.6a09e667f3bccp-1, + -0x1.87de2a6aea963p-2, -0x1.8f8b83c69a60ap-3, -0x1.917a6bc29b42cp-4, -0x1.91f65f10dd814p-5, + -0x1.92155f7a3667ep-6, -0x1.921d1fcdec784p-7, -0x1.921f0fe670071p-8, -0x1.921f8becca4bap-9, + -0x1.921faaee6472dp-10, -0x1.921fb2aecb36p-11, -0x1.921fb49ee4ea6p-12, -0x1.921fb51aeb57bp-13, + -0x1.921fb539ecf31p-14, -0x1.921fb541ad59ep-15, -0x1.921fb5439d73ap-16, -0x1.921fb544197ap-17, + -0x1.921fb544387bap-18, -0x1.921fb544403c1p-19, -0x1.921fb544422c2p-20, -0x1.921fb54442a83p-21, + -0x1.921fb54442c73p-22, -0x1.921fb54442cefp-23, -0x1.921fb54442d0ep-24, -0x1.921fb54442d15p-25, + -0x1.921fb54442d17p-26, -0x1.921fb54442d18p-27, -0x1.921fb54442d18p-28, -0x1.921fb54442d18p-29, + -0x1.921fb54442d18p-30, -0x1.921fb54442d18p-31, -0x1.921fb54442d18p-32, -0x1.921fb54442d18p-33, + -0x1.921fb54442d18p-34, -0x1.921fb54442d18p-35, -0x1.921fb54442d18p-36, -0x1.921fb54442d18p-37, + -0x1.921fb54442d18p-38, -0x1.921fb54442d18p-39, -0x1.921fb54442d18p-40, -0x1.921fb54442d18p-41, + -0x1.921fb54442d18p-42, -0x1.921fb54442d18p-43, -0x1.921fb54442d18p-44, -0x1.921fb54442d18p-45, + -0x1.921fb54442d18p-46, -0x1.921fb54442d18p-47, -0x1.921fb54442d18p-48, -0x1.921fb54442d18p-49, + -0x1.921fb54442d18p-50, -0x1.921fb54442d18p-51, -0x1.921fb54442d18p-52, -0x1.921fb54442d18p-53, + -0x1.921fb54442d18p-54, -0x1.921fb54442d18p-55, -0x1.921fb54442d18p-56, -0x1.921fb54442d18p-57, + -0x1.921fb54442d18p-58, -0x1.921fb54442d18p-59, -0x1.921fb54442d18p-60 }; + + /** Type of DFT. */ + protected final Norm normalization; + /** Inverse or forward. */ + protected final boolean inverse; + + /** + * @param normalization Normalization to be applied to the + * transformed data. + * @param inverse Whether to perform the inverse transform. + */ + public FastFourierTransform(final Norm normalization, + final boolean inverse) { + this.normalization = normalization; + this.inverse = inverse; + } + + /** + * @param normalization Normalization to be applied to the + * transformed data. + */ + public FastFourierTransform(final Norm normalization) { + this(normalization, false); + } + + /** + * Constructor + * + * @param inverse Whether to perform the inverse transform. + */ + public FastFourierTransform(final boolean inverse) { + this(Norm.STD, inverse); + } + + /** + * Constructor + */ + public FastFourierTransform() { + this(Norm.STD, false); + } + + /** + * Computes the standard transform of the data. + * Computation is done in place. + * Assumed layout of the input data: + *

    + *
  • {@code dataRI[0][i]}: Real part of the {@code i}-th data point,
    • + *
  • {@code dataRI[1][i]}: Imaginary part of the {@code i}-th data point.
    • + *
+ * + * @param dataRI Two-dimensional array of real and imaginary parts of the data. + * @throws IllegalArgumentException if the number of data points is not + * a power of two, if the number of rows of the specified array is not two, + * or the array is not rectangular. + */ + public void transformInPlace(final double[][] dataRI) { + if (dataRI.length != NUM_PARTS) { + throw new TransformException(TransformException.SIZE_MISMATCH, + dataRI.length, NUM_PARTS); + } + final double[] dataR = dataRI[0]; + final double[] dataI = dataRI[1]; + if (dataR.length != dataI.length) { + throw new TransformException(TransformException.SIZE_MISMATCH, + dataI.length, dataR.length); + } + final int n = dataR.length; + if (!ArithmeticUtils.isPowerOfTwo(n)) { + throw new TransformException(TransformException.NOT_POWER_OF_TWO, + Integer.valueOf(n)); + } + + if (n == 1) { + return; + } else if (n == 2) { + final double srcR0 = dataR[0]; + final double srcI0 = dataI[0]; + final double srcR1 = dataR[1]; + final double srcI1 = dataI[1]; + + // X_0 = x_0 + x_1 + dataR[0] = srcR0 + srcR1; + dataI[0] = srcI0 + srcI1; + // X_1 = x_0 - x_1 + dataR[1] = srcR0 - srcR1; + dataI[1] = srcI0 - srcI1; + + normalizeTransformedData(dataRI); + return; + } + + bitReversalShuffle2(dataR, dataI); + + // Do 4-term DFT. + if (inverse) { + for (int i0 = 0; i0 < n; i0 += 4) { + final int i1 = i0 + 1; + final int i2 = i0 + 2; + final int i3 = i0 + 3; + + final double srcR0 = dataR[i0]; + final double srcI0 = dataI[i0]; + final double srcR1 = dataR[i2]; + final double srcI1 = dataI[i2]; + final double srcR2 = dataR[i1]; + final double srcI2 = dataI[i1]; + final double srcR3 = dataR[i3]; + final double srcI3 = dataI[i3]; + + // 4-term DFT + // X_0 = x_0 + x_1 + x_2 + x_3 + dataR[i0] = srcR0 + srcR1 + srcR2 + srcR3; + dataI[i0] = srcI0 + srcI1 + srcI2 + srcI3; + // X_1 = x_0 - x_2 + j * (x_3 - x_1) + dataR[i1] = srcR0 - srcR2 + (srcI3 - srcI1); + dataI[i1] = srcI0 - srcI2 + (srcR1 - srcR3); + // X_2 = x_0 - x_1 + x_2 - x_3 + dataR[i2] = srcR0 - srcR1 + srcR2 - srcR3; + dataI[i2] = srcI0 - srcI1 + srcI2 - srcI3; + // X_3 = x_0 - x_2 + j * (x_1 - x_3) + dataR[i3] = srcR0 - srcR2 + (srcI1 - srcI3); + dataI[i3] = srcI0 - srcI2 + (srcR3 - srcR1); + } + } else { + for (int i0 = 0; i0 < n; i0 += 4) { + final int i1 = i0 + 1; + final int i2 = i0 + 2; + final int i3 = i0 + 3; + + final double srcR0 = dataR[i0]; + final double srcI0 = dataI[i0]; + final double srcR1 = dataR[i2]; + final double srcI1 = dataI[i2]; + final double srcR2 = dataR[i1]; + final double srcI2 = dataI[i1]; + final double srcR3 = dataR[i3]; + final double srcI3 = dataI[i3]; + + // 4-term DFT + // X_0 = x_0 + x_1 + x_2 + x_3 + dataR[i0] = srcR0 + srcR1 + srcR2 + srcR3; + dataI[i0] = srcI0 + srcI1 + srcI2 + srcI3; + // X_1 = x_0 - x_2 + j * (x_3 - x_1) + dataR[i1] = srcR0 - srcR2 + (srcI1 - srcI3); + dataI[i1] = srcI0 - srcI2 + (srcR3 - srcR1); + // X_2 = x_0 - x_1 + x_2 - x_3 + dataR[i2] = srcR0 - srcR1 + srcR2 - srcR3; + dataI[i2] = srcI0 - srcI1 + srcI2 - srcI3; + // X_3 = x_0 - x_2 + j * (x_1 - x_3) + dataR[i3] = srcR0 - srcR2 + (srcI3 - srcI1); + dataI[i3] = srcI0 - srcI2 + (srcR1 - srcR3); + } + } + + int lastN0 = 4; + int lastLogN0 = 2; + while (lastN0 < n) { + final int n0 = lastN0 << 1; + final int logN0 = lastLogN0 + 1; + final double wSubN0R = W_SUB_N_R[logN0]; + double wSubN0I = W_SUB_N_I[logN0]; + if (inverse) { + wSubN0I = -wSubN0I; + } + + // Combine even/odd transforms of size lastN0 into a transform of size N0 (lastN0 * 2). + for (int destEvenStartIndex = 0; destEvenStartIndex < n; destEvenStartIndex += n0) { + final int destOddStartIndex = destEvenStartIndex + lastN0; + + double wSubN0ToRR = 1; + double wSubN0ToRI = 0; + + for (int r = 0; r < lastN0; r++) { + final int destEvenStartIndexPlusR = destEvenStartIndex + r; + final int destOddStartIndexPlusR = destOddStartIndex + r; + + final double grR = dataR[destEvenStartIndexPlusR]; + final double grI = dataI[destEvenStartIndexPlusR]; + final double hrR = dataR[destOddStartIndexPlusR]; + final double hrI = dataI[destOddStartIndexPlusR]; + + final double a = wSubN0ToRR * hrR - wSubN0ToRI * hrI; + final double b = wSubN0ToRR * hrI + wSubN0ToRI * hrR; + // dest[destEvenStartIndex + r] = Gr + WsubN0ToR * Hr + dataR[destEvenStartIndexPlusR] = grR + a; + dataI[destEvenStartIndexPlusR] = grI + b; + // dest[destOddStartIndex + r] = Gr - WsubN0ToR * Hr + dataR[destOddStartIndexPlusR] = grR - a; + dataI[destOddStartIndexPlusR] = grI - b; + + // WsubN0ToR *= WsubN0R + final double nextWsubN0ToRR = wSubN0ToRR * wSubN0R - wSubN0ToRI * wSubN0I; + final double nextWsubN0ToRI = wSubN0ToRR * wSubN0I + wSubN0ToRI * wSubN0R; + wSubN0ToRR = nextWsubN0ToRR; + wSubN0ToRI = nextWsubN0ToRI; + } + } + + lastN0 = n0; + lastLogN0 = logN0; + } + + normalizeTransformedData(dataRI); + } + + /** + * {@inheritDoc} + * + * @throws IllegalArgumentException if the length of the data array is not a power of two. + */ + @Override + public Array apply(final double[] f) { + final double[][] dataRI = { + Arrays.copyOf(f, f.length), + new double[f.length] + }; + transformInPlace(dataRI); + return TransformUtils.createComplexArray(dataRI); + } + + /** + * {@inheritDoc} + * + * @throws IllegalArgumentException if the number of sample points + * {@code n} is not a power of two, if the lower bound is greater than, + * or equal to the upper bound, if the number of sample points {@code n} + * is negative + */ + @Override + public Array apply(final DoubleUnaryOperator f, + final double min, + final double max, + final int n) { + return apply(TransformUtils.sample(f, min, max, n)); + } + + /** + * {@inheritDoc} + * + * @throws IllegalArgumentException if the length of the data array is + * not a power of two. + */ + @Override + public Array apply(Array f) { + if (!ArithmeticUtils.isPowerOfTwo(f.getSize())) { + f = pad(f); + } + + final double[][] dataRI = TransformUtils.createRealImaginary(f); + transformInPlace(dataRI); + return TransformUtils.createComplexArray(dataRI); + } + + /** + * Applies normalization to the transformed data. + * + * @param dataRI Unscaled transformed data. + */ + private void normalizeTransformedData(final double[][] dataRI) { + if (normalization == Norm.NONE) { + return; + } + + final double[] dataR = dataRI[0]; + final double[] dataI = dataRI[1]; + final int n = dataR.length; + + switch (normalization) { + case STD: + if (inverse) { + final double scaleFactor = 1d / n; + for (int i = 0; i < n; i++) { + dataR[i] *= scaleFactor; + dataI[i] *= scaleFactor; + } + } + + break; + + case UNIT: + final double scaleFactor = 1d / Math.sqrt(n); + for (int i = 0; i < n; i++) { + dataR[i] *= scaleFactor; + dataI[i] *= scaleFactor; + } + + break; + + default: + throw new IllegalStateException(); // Should never happen. + } + } + + /** + * Performs identical index bit reversal shuffles on two arrays of + * identical size. + * Each element in the array is swapped with another element based + * on the bit-reversal of the index. + * For example, in an array with length 16, item at binary index 0011 + * (decimal 3) would be swapped with the item at binary index 1100 + * (decimal 12). + * + * @param a Array to be shuffled. + * @param b Array to be shuffled. + */ + private static void bitReversalShuffle2(double[] a, + double[] b) { + final int n = a.length; + final int halfOfN = n >> 1; + + int j = 0; + for (int i = 0; i < n; i++) { + if (i < j) { + // swap indices i & j + double temp = a[i]; + a[i] = a[j]; + a[j] = temp; + + temp = b[i]; + b[i] = b[j]; + b[j] = temp; + } + + int k = halfOfN; + while (k <= j && k > 0) { + j -= k; + k >>= 1; + } + j += k; + } + } + + /** + * Pads data so that its length is a power of 2 + * + * @param data the 1D array data + * @return the padded 1D array data + */ + private Array pad(Array data) { + int length = (int) data.getSize(); + if (ArithmeticUtils.isPowerOfTwo(length)) { + return data; + } + + int pLength = TransformUtils.nextPowerOfTwo(length); + Array padded = Array.factory(data.getDataType(), new int[]{pLength}); + for (int i = 0; i < length; i++) { + padded.setObject(i, data.getObject(i)); + } + + return padded; + } + + /** + * Normalization types. + */ + public enum Norm { + /** + * Should be passed to the constructor of {@link FastFourierTransform} + * to use the standard normalization convention. This normalization + * convention is defined as follows + *
    + *
  • forward transform: \( y_n = \sum_{k = 0}^{N - 1} x_k e^{-2 \pi i n k / N} \),
    • + *
  • inverse transform: \( x_k = \frac{1}{N} \sum_{n = 0}^{N - 1} y_n e^{2 \pi i n k / N} \),
    • + *
+ * where \( N \) is the size of the data sample. + */ + STD, + + /** + * Should be passed to the constructor of {@link FastFourierTransform} + * to use the unitary normalization convention. This normalization + * convention is defined as follows + *
    + *
  • forward transform: \( y_n = \frac{1}{\sqrt{N}} \sum_{k = 0}^{N - 1} x_k e^{-2 \pi i n k / N} \),
    • + *
  • inverse transform: \( x_k = \frac{1}{\sqrt{N}} \sum_{n = 0}^{N - 1} y_n e^{2 \pi i n k / N} \),
    • + *
+ * where \( N \) is the size of the data sample. + */ + UNIT, + + /** + * Not do normalization + */ + NONE; + } +} diff --git a/meteoinfo-math/src/main/java/org/meteoinfo/math/transform/FastFourierTransform2D.java b/meteoinfo-math/src/main/java/org/meteoinfo/math/transform/FastFourierTransform2D.java new file mode 100644 index 00000000..c1ce43a2 --- /dev/null +++ b/meteoinfo-math/src/main/java/org/meteoinfo/math/transform/FastFourierTransform2D.java @@ -0,0 +1,67 @@ +package org.meteoinfo.math.transform; + +import org.meteoinfo.ndarray.Array; +import org.meteoinfo.ndarray.DataType; +import org.meteoinfo.ndarray.InvalidRangeException; +import org.meteoinfo.ndarray.Range; +import org.meteoinfo.ndarray.math.ArrayMath; +import org.meteoinfo.ndarray.math.ArrayUtil; + +import java.util.Arrays; +import java.util.List; + +public class FastFourierTransform2D extends FastFourierTransform { + + /** + * Constructor + * @param inverse Whether is inverse transform + */ + public FastFourierTransform2D(boolean inverse) { + super(inverse); + } + + @Override + public Array apply(Array f) { + f = f.copyIfView(); + + int[] shape = f.getShape(); + int nRow = shape[0]; + int nCol = shape[1]; + Array r = Array.factory(DataType.COMPLEX, shape); + + try { + FastFourierTransform fastFourierTransform = new FastFourierTransform(this.normalization, this.inverse); + + //Transform rows + Range xRange = new Range(0, nCol - 1, 1); + Range yRange; + List ranges; + for (int i = 0; i < nRow; i++) { + yRange = new Range(i, i); + ranges = Arrays.asList(yRange, xRange); + Array data = ArrayMath.section(f, ranges).copy(); + data = fastFourierTransform.apply(data); + ArrayMath.setSection(r, ranges, data); + } + + //Transform cols + yRange = new Range(0, nRow - 1, 1); + for (int i = 0; i < nCol; i++) { + xRange = new Range(i, i); + ranges = Arrays.asList(yRange, xRange); + Array data = ArrayMath.section(r, ranges).copy(); + data = fastFourierTransform.apply(data); + ArrayMath.setSection(r, ranges, data); + } + } catch (InvalidRangeException e) { + e.printStackTrace(); + } + + return r; + } + + @Override + public Array apply(double[] f) { + return null; + } +} diff --git a/meteoinfo-math/src/main/java/org/meteoinfo/math/transform/TransformException.java b/meteoinfo-math/src/main/java/org/meteoinfo/math/transform/TransformException.java new file mode 100644 index 00000000..fd694ce6 --- /dev/null +++ b/meteoinfo-math/src/main/java/org/meteoinfo/math/transform/TransformException.java @@ -0,0 +1,37 @@ +package org.meteoinfo.math.transform; + +import java.text.MessageFormat; + +/** + * Exception class with constants for frequently used messages. + * Class is package-private (for internal use only). + */ +public class TransformException extends IllegalArgumentException { + + /** Error message for "out of range" condition. */ + public static final String FIRST_ELEMENT_NOT_ZERO = "First element ({0}) must be 0"; + /** Error message for "not strictly positive" condition. */ + public static final String NOT_STRICTLY_POSITIVE = "Number {0} is not strictly positive"; + /** Error message for "too large" condition. */ + public static final String TOO_LARGE = "Number {0} is larger than {1}"; + /** Error message for "size mismatch" condition. */ + public static final String SIZE_MISMATCH = "Size mismatch: {0} != {1}"; + /** Error message for "pow(2, n) + 1". */ + public static final String NOT_POWER_OF_TWO_PLUS_ONE = "{0} is not equal to 1 + pow(2, n), for some n"; + /** Error message for "pow(2, n)". */ + public static final String NOT_POWER_OF_TWO = "{0} is not equal to pow(2, n), for some n"; + + /** Serializable version identifier. */ + private static final long serialVersionUID = 20210522L; + + /** + * Create an exception where the message is constructed by applying + * the {@code format()} method from {@code java.text.MessageFormat}. + * + * @param message Message format (with replaceable parameters). + * @param formatArguments Actual arguments to be displayed in the message. + */ + TransformException(String message, Object... formatArguments) { + super(MessageFormat.format(message, formatArguments)); + } +} diff --git a/meteoinfo-math/src/main/java/org/meteoinfo/math/transform/TransformUtils.java b/meteoinfo-math/src/main/java/org/meteoinfo/math/transform/TransformUtils.java new file mode 100644 index 00000000..905f7224 --- /dev/null +++ b/meteoinfo-math/src/main/java/org/meteoinfo/math/transform/TransformUtils.java @@ -0,0 +1,237 @@ +package org.meteoinfo.math.transform; + +import java.util.function.DoubleUnaryOperator; + +import org.apache.commons.numbers.core.ArithmeticUtils; +import org.meteoinfo.ndarray.Array; +import org.meteoinfo.ndarray.Complex; +import org.meteoinfo.ndarray.DataType; + +/** + * Useful functions for the implementation of various transforms. + * Class is package-private (for internal use only). + */ +final class TransformUtils { + /** Number of array slots: 1 for "real" parts 1 for "imaginary" parts. */ + private static final int NUM_PARTS = 2; + + /** Utility class. */ + private TransformUtils() {} + + /** + * Multiply every component in the given real array by the + * given real number. The change is made in place. + * + * @param f Array to be scaled. + * @param d Scaling coefficient. + * @return a reference to the scaled array. + */ + static double[] scaleInPlace(double[] f, double d) { + for (int i = 0; i < f.length; i++) { + f[i] *= d; + } + return f; + } + + /** + * Multiply every component in the given complex array by the + * given real number. The change is made in place. + * + * @param f Array to be scaled. + * @param d Scaling coefficient. + * @return the scaled array. + */ + static Complex[] scaleInPlace(Complex[] f, double d) { + for (int i = 0; i < f.length; i++) { + f[i] = Complex.ofCartesian(d * f[i].getReal(), d * f[i].getImaginary()); + } + return f; + } + + + /** + * Builds a new two dimensional array of {@code double} filled with the real + * and imaginary parts of the specified {@link Complex} numbers. In the + * returned array {@code dataRI}, the data is laid out as follows + *
    + *
  • {@code dataRI[0][i] = dataC[i].getReal()},
    • + *
  • {@code dataRI[1][i] = dataC[i].getImaginary()}.
    • + *
+ * + * @param dataC Array of {@link Complex} data to be transformed. + * @return a two dimensional array filled with the real and imaginary parts + * of the specified complex input. + */ + static double[][] createRealImaginary(final Complex[] dataC) { + final double[][] dataRI = new double[2][dataC.length]; + final double[] dataR = dataRI[0]; + final double[] dataI = dataRI[1]; + for (int i = 0; i < dataC.length; i++) { + final Complex c = dataC[i]; + dataR[i] = c.getReal(); + dataI[i] = c.getImaginary(); + } + return dataRI; + } + + /** + * Builds a new two dimensional array of {@code double} filled with the real + * and imaginary parts of the specified {@link Complex} numbers. In the + * returned array {@code dataRI}, the data is laid out as follows + *
    + *
  • {@code dataRI[0][i] = dataC[i].getReal()},
    • + *
  • {@code dataRI[1][i] = dataC[i].getImaginary()}.
    • + *
+ * + * @param dataC Array of {@link Complex} data to be transformed. + * @return a two dimensional array filled with the real and imaginary parts + * of the specified complex input. + */ + static double[][] createRealImaginary(final Array dataC) { + dataC.copyIfView(); + + final double[][] dataRI = new double[2][(int) dataC.getSize()]; + final double[] dataR = dataRI[0]; + final double[] dataI = dataRI[1]; + for (int i = 0; i < dataC.getSize(); i++) { + final Complex c = dataC.getComplex(i); + dataR[i] = c.getReal(); + dataI[i] = c.getImaginary(); + } + return dataRI; + } + + /** + * Builds a new array of {@link Complex} from the specified two dimensional + * array of real and imaginary parts. In the returned array {@code dataC}, + * the data is laid out as follows + *
    + *
  • {@code dataC[i].getReal() = dataRI[0][i]},
    • + *
  • {@code dataC[i].getImaginary() = dataRI[1][i]}.
    • + *
+ * + * @param dataRI Array of real and imaginary parts to be transformed. + * @return a {@link Complex} array. + * @throws IllegalArgumentException if the number of rows of the specified + * array is not two, or the array is not rectangular. + */ + static Complex[] createComplex(final double[][] dataRI) { + if (dataRI.length != NUM_PARTS) { + throw new TransformException(TransformException.SIZE_MISMATCH, + dataRI.length, NUM_PARTS); + } + final double[] dataR = dataRI[0]; + final double[] dataI = dataRI[1]; + if (dataR.length != dataI.length) { + throw new TransformException(TransformException.SIZE_MISMATCH, + dataI.length, dataR.length); + } + + final int n = dataR.length; + final Complex[] c = new Complex[n]; + for (int i = 0; i < n; i++) { + c[i] = Complex.ofCartesian(dataR[i], dataI[i]); + } + return c; + } + + /** + * Builds a new array of {@link Complex} from the specified two dimensional + * array of real and imaginary parts. In the returned array {@code dataC}, + * the data is laid out as follows + *
    + *
  • {@code dataC[i].getReal() = dataRI[0][i]},
    • + *
  • {@code dataC[i].getImaginary() = dataRI[1][i]}.
    • + *
+ * + * @param dataRI Array of real and imaginary parts to be transformed. + * @return a {@link Complex} array. + * @throws IllegalArgumentException if the number of rows of the specified + * array is not two, or the array is not rectangular. + */ + static Array createComplexArray(final double[][] dataRI) { + if (dataRI.length != NUM_PARTS) { + throw new TransformException(TransformException.SIZE_MISMATCH, + dataRI.length, NUM_PARTS); + } + final double[] dataR = dataRI[0]; + final double[] dataI = dataRI[1]; + if (dataR.length != dataI.length) { + throw new TransformException(TransformException.SIZE_MISMATCH, + dataI.length, dataR.length); + } + + final int n = dataR.length; + final Array c = Array.factory(DataType.COMPLEX, new int[]{n}); + for (int i = 0; i < n; i++) { + c.setComplex(i, Complex.ofCartesian(dataR[i], dataI[i])); + } + return c; + } + + /** + * Samples the specified univariate real function on the specified interval. + *

+ * The interval is divided equally into {@code n} sections and sample points + * are taken from {@code min} to {@code max - (max - min) / n}; therefore + * {@code f} is not sampled at the upper bound {@code max}.

+ * + * @param f Function to be sampled + * @param min Lower bound of the interval (included). + * @param max Upper bound of the interval (excluded). + * @param n Number of sample points. + * @return the array of samples. + * @throws IllegalArgumentException if the lower bound {@code min} is + * greater than, or equal to the upper bound {@code max}, if the number + * of sample points {@code n} is negative. + */ + static double[] sample(DoubleUnaryOperator f, + double min, + double max, + int n) { + if (n <= 0) { + throw new TransformException(TransformException.NOT_STRICTLY_POSITIVE, + Integer.valueOf(n)); + } + if (min >= max) { + throw new TransformException(TransformException.TOO_LARGE, min, max); + } + + final double[] s = new double[n]; + final double h = (max - min) / n; + for (int i = 0; i < n; i++) { + s[i] = f.applyAsDouble(min + i * h); + } + return s; + } + + /** + * Check whether the number is a power of 2 + * @param number + * @return Whether is a power of 2 + */ + public static boolean isPowerOfTwo(int n) { + return n > 0L && (n & n - 1L) == 0L; + } + + /** + * Find next number which is a power of 2 + * @param number + * @return The next number which is a power of 2 + */ + public static int nextPowerOfTwo(int number) { + if (isPowerOfTwo(number)) { + return number; + } + + int x = number; + x = x - 1; + x |= x >> 1; + x |= x >> 2; + x |= x >> 4; + x |= x >> 8; + x |= x >> 16; + + return x + 1; + } +} diff --git a/meteoinfo-ndarray/src/main/java/org/meteoinfo/ndarray/ArrayComplex.java b/meteoinfo-ndarray/src/main/java/org/meteoinfo/ndarray/ArrayComplex.java index e5a36d0e..2114a05a 100644 --- a/meteoinfo-ndarray/src/main/java/org/meteoinfo/ndarray/ArrayComplex.java +++ b/meteoinfo-ndarray/src/main/java/org/meteoinfo/ndarray/ArrayComplex.java @@ -117,7 +117,11 @@ public class ArrayComplex extends Array { if (data != null) { storage = data; } else { + int n = (int) indexCalc.getSize(); storage = new Complex[(int) indexCalc.getSize()]; + for (int i = 0; i < n; i++) { + storage[i] = Complex.ZERO; + } } } diff --git a/pom.xml b/pom.xml index 369afd86..d635583c 100644 --- a/pom.xml +++ b/pom.xml @@ -35,7 +35,7 @@ UTF-8 1.8 - 3.9.5 + 3.9.6 8 8 8