add fft package

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<Triangle3D> triangles) {
+        LinkedHashMap<Vector3f, Integer> map = new LinkedHashMap<Vector3f, Integer>();
+        int n = triangles.size();
+        this.vertexIndices = new int[n];
+        Vector3f vector3f;
+        int idx = 0, vertexIdx = 0, triangleIdx = 0, index, ii;
+        List<Integer> 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<Vector3f, Integer> 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

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

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)) {
+                drawTriMeshGraphic(gl, (TriMeshGraphic) 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();
             }
         } 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<PointZ>) 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<PointZ>) 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<PointZ>) 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<PointZ> newPList;
-                gl.glBegin(GL2.GL_LINE_STRIP);
-                for (int h = 0; h < aPG.getHoleLines().size(); h++) {
-                    newPList = (java.util.List<PointZ>) 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()) {

meteoinfo-chart/src/main/java/org/meteoinfo/chart/jogl/Triangle.java

@@ -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<PointZ> getPoints() {
-        List<PointZ> 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 + "]";
-    }
-}

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<Triangle> getTriangles(PolygonZ polygon)
+    public List<Triangle3D> getTriangles(PolygonZ polygon)
             throws TesselationException {
 
         makeTriangles(polygon);
@@ -83,7 +84,7 @@ public class TriangleTessellator {
     }
 
     public interface TessellatorListener {
-        public void onTesselationDone(List<Triangle> triangles);
+        public void onTesselationDone(List<Triangle3D> triangles);
 
         public void onTesselationError(TesselationException err);
     }
@@ -163,9 +164,9 @@ public class TriangleTessellator {
      */
     class TessellationCallback extends GLUtessellatorCallbackAdapter {
 
-        protected List<Triangle> triangles = new ArrayList<>();
+        protected List<Triangle3D> 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;

meteoinfo-lab/milconfig.xml

@@ -1,32 +1,36 @@
 <?xml version="1.0" encoding="UTF-8" standalone="no"?>
 <MeteoInfo File="milconfig.xml" Type="configurefile">
-  <Path OpenPath="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\model">
+  <Path OpenPath="D:\Working\MIScript\Jython\mis\common_math\fft">
-    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\map"/>
+    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\surf"/>
-    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\map\webmap"/>
+    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\text"/>
-    <RecentFolder Folder="D:\Working\MIScript\Jython\mis"/>
+    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\volume"/>
-    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\array"/>
-    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\funny"/>
-    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\animation"/>
-    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types"/>
-    <RecentFolder Folder="D:\Temp\test"/>
-    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\plot"/>
-    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\slice"/>
-    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\geoshow"/>
-    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\mesh"/>
-    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\fill"/>
-    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl"/>
     <RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\model"/>
+    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\isosurface"/>
+    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl"/>
+    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\3d"/>
+    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\pie"/>
+    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types"/>
+    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\map"/>
+    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\map\geoshow"/>
+    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\array"/>
+    <RecentFolder Folder="D:\Working\MIScript\Jython\mis"/>
+    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\common_math"/>
+    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\common_math\fft"/>
   </Path>
   <File>
     <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\plot_types\3d\jogl\model\obj_airplane.py"/>
+      <OpenedFile File="D:\Working\MIScript\Jython\mis\common_math\fft\fft_3.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\plot_types\3d\jogl\model\obj_airplane.py"/>
+      <RecentFile File="D:\Working\MIScript\Jython\mis\common_math\fft\fft_3.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>
@@ -34,5 +38,5 @@
   </Font>
   <LookFeel DockWindowDecorated="true" LafDecorated="true" Name="FlatDarkLaf"/>
   <Figure DoubleBuffering="true"/>
-  <Startup MainFormLocation="-7,0" MainFormSize="1348,825"/>
+  <Startup MainFormLocation="-7,0" MainFormSize="1373,822"/>
 </MeteoInfo>

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__)

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))
 

meteoinfo-lab/pylib/mipylib/numeric/fft/__init__.py

@@ -0,0 +1,2 @@
+from ._fft import *
+from ._helper import *

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)

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

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.
+ * <p>
+ * Such transforms include {@link FastSineTransform sine transform},
+ * {@link FastCosineTransform cosine transform} or {@link
+ * FastHadamardTransform Hadamard transform}.
+ */
+public interface ComplexTransform extends UnaryOperator<Array> {
+    /**
+     * 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));
+    }
+}

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;
+    }
+}

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 <em>Applied Numerical Linear
+ * Algebra</em>, ISBN 0898713897, chapter 6.
+ * <p>
+ * There are several variants of the discrete Fourier transform, with various
+ * normalization conventions, which are specified by the parameter
+ * {@link Norm}.
+ * <p>
+ * 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,
+ * <i>On computing the discrete Fourier transform</i>, 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:
+     * <ul>
+     *   <li>{@code dataRI[0][i]}: Real part of the {@code i}-th data point,</li>
+     *   <li>{@code dataRI[1][i]}: Imaginary part of the {@code i}-th data point.</li>
+     * </ul>
+     *
+     * @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 <em>standard</em> normalization convention. This normalization
+         * convention is defined as follows
+         * <ul>
+         * <li>forward transform: \( y_n = \sum_{k = 0}^{N - 1} x_k e^{-2 \pi i n k / N} \),</li>
+         * <li>inverse transform: \( x_k = \frac{1}{N} \sum_{n = 0}^{N - 1} y_n e^{2 \pi i n k / N} \),</li>
+         * </ul>
+         * where \( N \) is the size of the data sample.
+         */
+        STD,
+
+        /**
+         * Should be passed to the constructor of {@link FastFourierTransform}
+         * to use the <em>unitary</em> normalization convention. This normalization
+         * convention is defined as follows
+         * <ul>
+         * <li>forward transform: \( y_n = \frac{1}{\sqrt{N}} \sum_{k = 0}^{N - 1} x_k e^{-2 \pi i n k / N} \),</li>
+         * <li>inverse transform: \( x_k = \frac{1}{\sqrt{N}} \sum_{n = 0}^{N - 1} y_n e^{2 \pi i n k / N} \),</li>
+         * </ul>
+         * where \( N \) is the size of the data sample.
+         */
+        UNIT,
+
+        /**
+         * Not do normalization
+         */
+        NONE;
+    }
+}

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<Range> 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;
+    }
+}

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));
+    }
+}

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
+     * <ul>
+     * <li>{@code dataRI[0][i] = dataC[i].getReal()},</li>
+     * <li>{@code dataRI[1][i] = dataC[i].getImaginary()}.</li>
+     * </ul>
+     *
+     * @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
+     * <ul>
+     * <li>{@code dataRI[0][i] = dataC[i].getReal()},</li>
+     * <li>{@code dataRI[1][i] = dataC[i].getImaginary()}.</li>
+     * </ul>
+     *
+     * @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
+     * <ul>
+     * <li>{@code dataC[i].getReal() = dataRI[0][i]},</li>
+     * <li>{@code dataC[i].getImaginary() = dataRI[1][i]}.</li>
+     * </ul>
+     *
+     * @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
+     * <ul>
+     * <li>{@code dataC[i].getReal() = dataRI[0][i]},</li>
+     * <li>{@code dataC[i].getImaginary() = dataRI[1][i]}.</li>
+     * </ul>
+     *
+     * @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.
+     * <p>
+     * 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}.</p>
+     *
+     * @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;
+    }
+}

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;
+            }
         }
     }
 

pom.xml

@@ -35,7 +35,7 @@
     <properties>
         <project.build.sourceEncoding>UTF-8</project.build.sourceEncoding>
         <java.version>1.8</java.version>
-        <revision>3.9.5</revision>
+        <revision>3.9.6</revision>
         <maven.compiler.source>8</maven.compiler.source>
         <maven.compiler.target>8</maven.compiler.target>
         <maven.compiler.release>8</maven.compiler.release>