Package net.algart.matrices.spectra

Class SeparableFastHartleyTransform

java.lang.Object
net.algart.matrices.spectra.AbstractSpectralTransform
net.algart.matrices.spectra.SeparableFastHartleyTransform
All Implemented Interfaces:
SpectralTransform
public class SeparableFastHartleyTransform extends AbstractSpectralTransform implements SpectralTransform

Fast Hartley transform (FHT) (in multidimensional case — the separable fast Hartley transform). This class implements standard one-dimensional FHT algorithm over an abstract SampleArray. It is generalized to multidimensional case by the simplest way, implemented in AbstractSpectralTransform class (applying the transform separably to each dimension); the resulting transformation for 2- or multidimensional AlgART numeric matrices is usually called separable fast Hartley transform (SFHT). The samples, processed by this class, can be both real or complex (areComplexSamplesRequired() method returns false). This class is especially useful in a case of real samples. In this case it is performed faster than the classic FFT, FastFourierTransform class, because there are no needs to allocate and process arrays of imaginary parts. The simple relation between Hartley and Fourier transform (see below) allows to use this transform almost in all areas where Fourier transform is applicable.

More precisely, this class implements the classic fast "butterfly" algorithm (FHT) for calculating discrete Hartley transform (DHT), described at http://en.wikipedia.org/wiki/Discrete_Hartley_transform.

Namely, let x0,x1,...,xN−1 are some real or complex samples (represented by abstract SampleArray), and H0,H1,...,HN−1 are their Hartley spectrum: the result of DHT. Let's designate cas θ = cos θ + sin θ. This class implements two possible definitions of DHT:

  1. direct transform is Hk = ∑(0≤n<N) xn cas(2knπ/N), inverse transform is xn = N −1 ∑(0≤k<N) Hk cas(2knπ/N).
  2. direct transform is Hk = N −1 ∑(0≤n<N) xn cas(2knπ/N), inverse transform is xn = ∑(0≤k<N) Hk cas(2knπ/N).

The only difference is when to normalize the result: while inverse transform (case 1) or direct transform (case 2). The Wikipedia offers formulas of the 1st case. This class allows to calculate both variants: the 1st case is chosen if the normalizeDirectTransform argument of the constructors is false or if this class is created by a constructor without this argument (it is the default behaviour), the 2nd case is chosen if the normalizeDirectTransform argument of the constructors is true.

The very useful feature of DHT is that for real samples xk the Hartley spectrum Hk is also real — unlike DFT, when even real samples lead to complex spectrum. As a result, the transformation algorithms in this class can process real arrays and matrices, without imaginary parts. In this case, they work in two and even more times faster than FFT algorithms, implemented in FastFourierTransform, and do not require allocating additional memory for storing imaginary parts of the complex numbers.

The formulas above correspond to one-dimensional transforms and specify the results of directTransform / inverseTransform methods. They are generalized to multidimensional case by default algorithms, implemented in AbstractSpectralTransform class, i.e. by applying the transform separably to each dimension. It leads to so-called multidimensional separable discrete Hartley transformations (SDHT). Below are the formulas for 2-dimensional separable discrete Hartley transformation of the matrix xij (0≤i<M, 0≤j<N) for the case 1 (normalizing the inverse transform):

direct: Hij = ∑(0≤m<M) ∑(0≤n<N) xmn cas(2imπ/M) cas(2jnπ/N),
inverse: xmn = (MN) −1 ∑(0≤i<M) ∑(0≤j<N) xij cas(2imπ/M) cas(2jnπ/N).

There is the simple relation between classic DFT (discrete Fourier transform) and SDHT (separable discrete Hartley transform).

Let's consider one-dimensional case (usual DHT). Let x0,x1,...,xN−1 are some real or complex samples, F0,F1,...,FN−1 are their Fourier spectrum and H0,H1,...,HN−1 are their Hartley spectrum. Let i is the usual imaginary unit. For simplicity, let's consider that F−k=FN−k, H−k=HN−k, k=1,2,... Then:

Fk = (Hk+H−k)/2 − i (Hk−H−k)/2,
Hk = (Fk+F−k)/2 + i (Fk−F−k)/2,
in a case of real samples: Hk = Re Fk − Im Fk

(of course, we consider the same definition, 1 or 2, for both DFT and SDHT spectra).

In 2-dimensional case, the relation between DFT and SDHT is the following (we similarly suppose that F−i, j=FM−i, j, Fi,−j=Fi, N−j, H−i, j=HM−i, j, Hi,−j=Hi, N−j):

Fi, j = (Hi,−j+H−i, j)/2 − i (Hi, j−H−i,−j)/2,
Hi, j = (Fi,−j+F−i, j)/2 + i (Fi, j−F−i,−j)/2,
in a case of real samples: Hi, j = Re Fi,−j − Im Fi, j.

In the common n-dimensional case, there are similar formulas, which express Fi, j,...,k through a linear combination of 2n numbers H± i,± j,...,± k and, vice versa, express Hi, j,...,k through a linear combination of 2n numbers F± i,± j,...,± k.

This class contains the ready for use methods, allowing to convert n-dimensional separable Hartley spectrum to Fourier one and vice versa, n=1,2,3,...:

If it is necessary to get the Fourier spectrum of some real matrix, probably process it and transform the Fourier spectrum back to the real matrix, you can use a combination of SHFT, provided by this class, and the conversion methods listed above (cases of real matrices). But if all that you need is to calculate a convolution of two real matrices, there is a better way: see below.

One-dimensional Hartley transform, defined by the formulas 1 and 2 above, complies with the convolution theorem. Namely, let p0,p1,...,pN−1 is the first complex or real numeric function and q0,q1,...,qN−1 is the second complex or real function, and c0,c1,...,cN−1 is their (complex or real) convolution, defined as:

ck = ∑(0≤n<N) pnqk−n

(here and below we consider that Z−k=ZN−k for all samples and spectra). Also, let P0,P1,...,PN−1, Q0,Q1,...,QN−1 and C0,C1,...,CN−1 are Hartley spectra of these functions. Then:

  1. Ck = (PkQ−k+P−kQk)/2 + (PkQk−P−kQ−k)/2, if the spectra were calculated according formula 1 above (default method);
  2. Ck = N ((PkQ−k+P−kQk)/2 + (PkQk−P−kQ−k)/2), if the spectra were calculated according formula 2 above.

There are similar formulas in the common n-dimensional case, allowing to express the separable Hartley spectrum of the convolution of two n-dimensional matrices via the spectra of the source matrices. In particular, in the 2-dimensional case:

  1. Ci, j = ((Pi, j+P−i,−j) (Qi, j+Q−i,−j) − (Pi,−j−P−i, j) (Qi,−j−Q−i, j) + (Pi, j−P−i,−j) (Qi,−j+Q−i, j) + (Pi,−j+P−i, j) (Qi, j−Q−i,−j))/4, if the spectra were calculated according formula 1 above (default method);
  2. Ci, j = the same expression multiplied by MN (M and N are the dimensions of the matrices), if the spectra were calculated according formula 2 above.

This class contains the ready for use methods, allowing to calculate a spectrum of convolution C on the base of the given spectra P and Q of two source numeric matrices x and y according the formulas A, C and their generalization for any number of dimensions:

So, if you need to calculate a convolution of some real matrices, for example, for goals of linear filtering, you can use the SFHT transform and the spectrumOfConvolution method, provided by this class: it is much better idea than using FastFourierTransform class.

Please note: in the one-dimensional case, the spectral transofmation algorithms, implemented by directTransformMatrix / inverseTransformMatrix methods of this class, work with normal (i.e. high) performance only if the passed one-dimensional AlgART matrices are stored in SimpleMemoryModel (more precisely, if they are directly accessible). In other case, each access to every sample leads to calling accessing methods getDouble and setDouble, which can work slowly in non-simple memory models like LargeMemoryModel. There is the same problem for directTransform / inverseTransform methods, if the passed sample arrays are created via RealScalarSampleArray.asSampleArray or ComplexScalarSampleArray.asSampleArray methods on the base of updatable AlgART arrays, created by memory model other than SimpleMemoryModel.

For n-dimensional matrices (n≥2), this problem usually does not occur at all, even for non-simple memory models, if you use standard implementations of directTransformMatrix / inverseTransformMatrix from AbstractSpectralTransform class: these implementations automatically download necessary parts of the matrix into SimpleMemoryModel. This problem also does not occur while using conversion methods separableHartleyToFourier(ArrayContext, Matrix, Matrix, Matrix), separableHartleyToFourier(ArrayContext, Matrix, Matrix, Matrix, Matrix), fourierToSeparableHartley(ArrayContext, Matrix, Matrix, Matrix), fourierToSeparableHartley(ArrayContext, Matrix, Matrix, Matrix, Matrix) and methods of calculation of the spectrum of convolution spectrumOfConvolution(ArrayContext, Matrix, Matrix, Matrix) and spectrumOfConvolution(ArrayContext, Matrix, Matrix, Matrix, Matrix, Matrix, Matrix), if all processed matrices have the same float or double element types.

  • Constructor Details

  • Method Details

    • separableHartleyToFourier

      public void separableHartleyToFourier(ArrayContext context, Matrix<? extends UpdatablePNumberArray> fRe, Matrix<? extends UpdatablePNumberArray> fIm, Matrix<? extends PNumberArray> h)
      Converts the separable Hartley spectrum H of some real n-dimensional matrix into the (complex) Fourier spectrum F of the same matrix. See the comments to this class about the relation formulas between separable Hartley and Fourier spectra.

      The complex matrix F is represented as a pair of AlgART matrices (fRe,fIm): the corresponding elements of these 2 matrices contain the real and imaginary parts of the corresponding elements of the complex matrix F. The real matrix H is passed as an AlgART matrix h.

      All matrices, passed to this method, must have equal dimensions. The element type of the passed matrices can be different, but we recommend using the same float or double element type for all matrices. There are no restrictions for the dimensions of the passed matrices: isLengthAllowed(long) method is not used here.

      This method works correctly, if you pass the same matrix as fRe / fIm and h.

      If you need to convert spectrum in a case of one-dimensional numeric AlgART arrays, you just need to convert them into one-dimensional AlgART matrices by Matrices.matrix(Array, long...) call, for example: Matrices.matrix(array, array.length()).

      Parameters:
      context - the context that will be used by this algorithm; may be null (see comments to SpectralTransform).
      fRe - the real parts of the elements of the resulting matrix (Fourier spectrum).
      fIm - the imaginary parts of the elements of the resulting matrix (Fourier spectrum).
      h - the source real matrix (separable Hartley spectrum).
      Throws:
      NullPointerException - if one of fRe, fIm, h arguments is null.
      SizeMismatchException - if some of the passed matrices have different dimensions.
      See Also:

    • separableHartleyToFourier

      public void separableHartleyToFourier(ArrayContext context, Matrix<? extends UpdatablePNumberArray> fRe, Matrix<? extends UpdatablePNumberArray> fIm, Matrix<? extends PNumberArray> hRe, Matrix<? extends PNumberArray> hIm)
      Converts the separable Hartley spectrum H of some complex n-dimensional matrix into the (complex) Fourier spectrum F of the same matrix. See the comments to this class about the relation formulas between separable Hartley and Fourier spectra.

      The complex matrix F is represented as a pair of AlgART matrices (fRe,fIm): the corresponding elements of these 2 matrices contain the real and imaginary parts of the corresponding elements of the complex matrix F. Similarly, the complex matrix H is represented as a pair of AlgART matrices (hRe,hIm).

      All matrices, passed to this method, must have equal dimensions. The element type of the passed matrices can be different, but we recommend using the same float or double element type for all matrices. There are no restrictions for the dimensions of the passed matrices: isLengthAllowed(long) method is not used here.

      This method works correctly, if you pass the same complex matrix as F and H. So, you can calculate and return the result in the source matrices.

      If you need to convert spectrum in a case of one-dimensional numeric AlgART arrays, you just need to convert them into one-dimensional AlgART matrices by Matrices.matrix(Array, long...) call, for example: Matrices.matrix(array, array.length()).

      Parameters:
      context - the context that will be used by this algorithm; may be null (see comments to SpectralTransform).
      fRe - the real parts of the elements of the resulting matrix (Fourier spectrum).
      fIm - the imaginary parts of the elements of the resulting matrix (Fourier spectrum).
      hRe - the real parts of the elements of the source matrix (separable Hartley spectrum).
      hIm - the imaginary parts of the elements of the source matrix (separable Hartley spectrum).
      Throws:
      NullPointerException - if one of fRe, fIm, hRe, hIm arguments is null.
      SizeMismatchException - if some of the passed matrices have different dimensions.
      See Also:

    • fourierToSeparableHartley

      public void fourierToSeparableHartley(ArrayContext context, Matrix<? extends UpdatablePNumberArray> h, Matrix<? extends PNumberArray> fRe, Matrix<? extends PNumberArray> fIm)
      Converts the Fourier spectrum F of some real n-dimensional matrix into the (real) separable Hartley spectrum H of the same matrix. See the comments to this class about the relation formulas between separable Hartley and Fourier spectra. If the passed Fourier spectrum is not a spectrum of a real matrix (in other words, if the inverse Fourier transform of F matrix contains nonzero imaginary parts), then this method still correctly calculates the real parts of the separable Hartley spectrum H.

      The complex matrix F is represented as a pair of AlgART matrices (fRe,fIm): the corresponding elements of these 2 matrices contain the real and imaginary parts of the corresponding elements of the complex matrix F. The real matrix H (or the real parts of H, if the passed F matrix is not a spectrum of a real matrix) is passed as an AlgART matrix h.

      All matrices, passed to this method, must have equal dimensions. The element type of the passed matrices can be different, but we recommend using the same float or double element type for all matrices. There are no restrictions for the dimensions of the passed matrices: isLengthAllowed(long) method is not used here.

      This method works correctly, if you pass the same matrix as fRe / fIm and h.

      If you need to convert spectrum in a case of one-dimensional numeric AlgART arrays, you just need to convert them into one-dimensional AlgART matrices by Matrices.matrix(Array, long...) call, for example: Matrices.matrix(array, array.length()).

      Parameters:
      context - the context that will be used by this algorithm; may be null (see comments to SpectralTransform).
      h - the resulting real matrix (separable Hartley spectrum).
      fRe - the real parts of the elements of the source matrix (Fourier spectrum).
      fIm - the imaginary parts of the elements of the source matrix (Fourier spectrum).
      Throws:
      NullPointerException - if one of h, fRe, fIm arguments is null.
      SizeMismatchException - if some of the passed matrices have different dimensions.
      See Also:

    • fourierToSeparableHartley

      public void fourierToSeparableHartley(ArrayContext context, Matrix<? extends UpdatablePNumberArray> hRe, Matrix<? extends UpdatablePNumberArray> hIm, Matrix<? extends PNumberArray> fRe, Matrix<? extends PNumberArray> fIm)
      Converts the Fourier spectrum F of some complex n-dimensional matrix into the (complex) separable Hartley spectrum H of the same matrix. See the comments to this class about the relation formulas between separable Hartley and Fourier spectra.

      The complex matrix F is represented as a pair of AlgART matrices (fRe,fIm): the corresponding elements of these 2 matrices contain the real and imaginary parts of the corresponding elements of the complex matrix F. Similarly, the complex matrix H is represented as a pair of AlgART matrices (hRe,hIm).

      All matrices, passed to this method, must have equal dimensions. The element type of the passed matrices can be different, but we recommend using the same float or double element type for all matrices. There are no restrictions for the dimensions of the passed matrices: isLengthAllowed(long) method is not used here.

      This method works correctly, if you pass the same complex matrix as F and H. So, you can calculate and return the result in the source matrices.

      If you need to convert spectrum in a case of one-dimensional numeric AlgART arrays, you just need to convert them into one-dimensional AlgART matrices by Matrices.matrix(Array, long...) call, for example: Matrices.matrix(array, array.length()).

      Parameters:
      context - the context that will be used by this algorithm; may be null (see comments to SpectralTransform).
      hRe - the real parts of the elements of the resulting matrix (separable Hartley spectrum).
      hIm - the imaginary parts of the elements of the resulting matrix (separable Hartley spectrum).
      fRe - the real parts of the elements of the source matrix (Fourier spectrum).
      fIm - the imaginary parts of the elements of the source matrix (Fourier spectrum).
      Throws:
      NullPointerException - if one of hRe, hIm, fRe, fIm arguments is null.
      SizeMismatchException - if some of the passed matrices have different dimensions.
      See Also:

    • spectrumOfConvolution

      public void spectrumOfConvolution(ArrayContext context, Matrix<? extends UpdatablePNumberArray> c, Matrix<? extends PNumberArray> p, Matrix<? extends PNumberArray> q)
      Calculates C, the separable Hartley spectrum of the convolution of some two real matrices, on the base of P and Q — the separable Hartley spectra of these two real matrices.

      The real matrices P, Q, C are passed as AlgART matrices p, q, c.

      All matrices, passed to this method, must have equal dimensions. The element type of the passed matrices can be different, but we recommend using the same float or double element type for all matrices. There are no restrictions for the dimensions of the passed matrices: isLengthAllowed(long) method is not used here.

      This method works correctly, if you pass the same complex matrix as P and Q, or as P and C, or as Q and C, or even as all three matrices. So, you can calculate and return the result in one of the source matrices.

      If you need to calculate the Hartley spectrum of convolution for a case of one-dimensional numeric AlgART arrays, you just need to convert them into one-dimensional AlgART matrices by Matrices.matrix(Array, long...) call, for example: Matrices.matrix(array, array.length()).

      Parameters:
      context - the context that will be used by this algorithm; may be null (see comments to SpectralTransform).
      c - the resulting matrix (spectrum of the convolution).
      p - the spectrum of the 1st matrix.
      q - the spectrum of the 2nd matrix.
      Throws:
      NullPointerException - if one of c, p, q arguments is null.
      SizeMismatchException - if some of the passed matrices have different dimensions.
      See Also:

    • spectrumOfConvolution

      public void spectrumOfConvolution(ArrayContext context, Matrix<? extends UpdatablePNumberArray> cRe, Matrix<? extends UpdatablePNumberArray> cIm, Matrix<? extends PNumberArray> pRe, Matrix<? extends PNumberArray> pIm, Matrix<? extends PNumberArray> qRe, Matrix<? extends PNumberArray> qIm)
      Calculates C, the separable Hartley spectrum of the convolution of some two complex matrices, on the base of P and Q — the separable Hartley spectra of these two complex matrices.

      The complex matrix P is represented as a pair of AlgART matrices (pRe,pIm): the corresponding elements of these 2 matrices contain the real and imaginary parts of the corresponding elements of the complex matrix P. Similarly, the complex matrix Q is represented as a pair of AlgART matrices (qRe,qIm), and the complex matrix C is represented as a pair of AlgART matrices (cRe,cIm).

      All matrices, passed to this method, must have equal dimensions. The element type of the passed matrices can be different, but we recommend using the same float or double element type for all matrices. There are no restrictions for the dimensions of the passed matrices: isLengthAllowed(long) method is not used here.

      This method works correctly, if you pass the same complex matrix as P and Q, or as P and C, or as Q and C, or even as all three matrices. So, you can calculate and return the result in one of the source matrices.

      If you need to calculate the Hartley spectrum of convolution for a case of one-dimensional numeric AlgART arrays, you just need to convert them into one-dimensional AlgART matrices by Matrices.matrix(Array, long...) call, for example: Matrices.matrix(array, array.length()).

      Parameters:
      context - the context that will be used by this algorithm; may be null (see comments to SpectralTransform).
      cRe - the real parts of the elements of the resulting matrix (spectrum of the convolution).
      cIm - the imaginary parts of the elements of the resulting matrix (spectrum of the convolution).
      pRe - the real parts of the elements of the spectrum of the 1st matrix.
      pIm - the imaginary parts of the elements of the spectrum of the 1st matrix.
      qRe - the real parts of the elements of the spectrum of the 2nd matrix.
      qIm - the imaginary parts of the elements of the spectrum of the 2nd matrix.
      Throws:
      NullPointerException - if one of cRe, cIm, pRe, pIm, qRe, qIm arguments is null.
      SizeMismatchException - if some of the passed matrices have different dimensions.
      See Also:

    • isLengthAllowed

      public final boolean isLengthAllowed(long length)
      Description copied from interface: SpectralTransform
      Returns true if the specified argument is an allowed dimension for arrays or matrices, transformed by directTransform, inverseTransform, directTransformMatrix or inverseTransformMatrix method.

      More precisely, if this method returns false for the length of a sample array, passed to 1st or 2nd methods, or for some dimension of some matrix, passed to 3rd or 4th method, then those methods throw IllegalArgumentException. In other case, those methods will process that passed data.

      In both implementations of this interface, offered by this package, this method returns true if the passed length is a power of two (2k).

      If the length argument is negative, the result of this method is unspecified. It is not a problem, because lengths of sample arrays and dimensions of AlgART matrices cannot be negative.

      Specified by:
      isLengthAllowed in interface SpectralTransform
      Specified by:
      isLengthAllowed in class AbstractSpectralTransform
      Parameters:
      length - the checked length or matrix dimension.
      Returns:
      whether the specified argument is an allowed dimension for arrays or matrices, trasformed by this transformation.

    • areComplexSamplesRequired

      public boolean areComplexSamplesRequired()
      Description copied from interface: SpectralTransform
      Returns true if the transformation methods of this class (directTransform, inverseTransform, directTransformMatrix, inverseTransformMatrix) can process only complex samples, false if the real samples are also allowed.

      More precisely, if this method returns true, then the methods directTransform / inverseTransform checks, whether SampleArray.isComplex() method returns true for the samples argument, and the methods directTransformMatrix / inverseTransformMatrix checks, whether the matrixIm argument is not null. If this condition is not fulfilled, these methods throw UnsupportedOperationException. In other case, these methods work normally.

      In implementations, offered by this package, this method returns true in FastFourierTransform class and false in SeparableFastHartleyTransform class.

      Specified by:
      areComplexSamplesRequired in interface SpectralTransform
      Specified by:
      areComplexSamplesRequired in class AbstractSpectralTransform
      Returns:
      true if this class can transform complex samples only, false if real samples can be transformed too.

    • unallowedLengthMessage

      protected String unallowedLengthMessage()
      Description copied from class: AbstractSpectralTransform
      Retrurns a message used while throwing IllegalArgumentException by methods of this class in a case, when the length of the samples array or some of the matrix dimensions is not allowed according to AbstractSpectralTransform.isLengthAllowed(long) method. Typical examples of this message (implemented in FastFourierTransform and SeparableFastHartleyTransform classes): "FFT algorithm can process only 2^k elements" or "FHT algorithm can process only 2^k elements".
      Specified by:
      unallowedLengthMessage in class AbstractSpectralTransform
      Returns:
      a message used while thrown exception if AbstractSpectralTransform.isLengthAllowed(long) method returns false.

    • transform

      protected final void transform(ArrayContext context, SampleArray samples, boolean inverse)
      Description copied from class: AbstractSpectralTransform
      Actually performs the 1-dimensional transform of the sample array, direct or inverse.

      It is called from directTransform / inverseTransform methods. In this case, there is a guarantee that: 1) samples!=null; 2) if AbstractSpectralTransform.areComplexSamplesRequired(), then samples.isComplex() returns true; 3) AbstractSpectralTransform.isLengthAllowed(long) returns true for samples.length().

      Specified by:
      transform in class AbstractSpectralTransform
      Parameters:
      context - the context that will be used by this algorithm; may be null (see comments to SpectralTransform).
      samples - the transformed samples.
      inverse - true if this method implements the inverse transform, false if this method implements the direct transform.