FFT1D.h

Go to the documentation of this file.

//     ____   ______ __
//    / __ \ / ____// /
//   / /_/ // /    / /
//  / ____// /___ / /___   PixInsight Class Library
// /_/     \____//_____/   PCL 2.7.0
// ----------------------------------------------------------------------------
// pcl/FFT1D.h - Released 2024-06-18T15:48:54Z
// ----------------------------------------------------------------------------
// This file is part of the PixInsight Class Library (PCL).
// PCL is a multiplatform C++ framework for development of PixInsight modules.
//
// Copyright (c) 2003-2024 Pleiades Astrophoto S.L. All Rights Reserved.
//
// Redistribution and use in both source and binary forms, with or without
// modification, is permitted provided that the following conditions are met:
//
// 1. All redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//
// 2. All redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//
// 3. Neither the names "PixInsight" and "Pleiades Astrophoto", nor the names
//    of their contributors, may be used to endorse or promote products derived
//    from this software without specific prior written permission. For written
//    permission, please contact info@pixinsight.com.
//
// 4. All products derived from this software, in any form whatsoever, must
//    reproduce the following acknowledgment in the end-user documentation
//    and/or other materials provided with the product:
//
//    "This product is based on software from the PixInsight project, developed
//    by Pleiades Astrophoto and its contributors (https://pixinsight.com/)."
//
//    Alternatively, if that is where third-party acknowledgments normally
//    appear, this acknowledgment must be reproduced in the product itself.
//
// THIS SOFTWARE IS PROVIDED BY PLEIADES ASTROPHOTO AND ITS CONTRIBUTORS
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
// TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL PLEIADES ASTROPHOTO OR ITS
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, BUSINESS
// INTERRUPTION; PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; AND LOSS OF USE,
// DATA OR PROFITS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
// ----------------------------------------------------------------------------
 
#ifndef __PCL_FFT1D_h
#define __PCL_FFT1D_h
 
 
/*
 * The FFT routines used in this PCL version are based on the KISS FFT library
 * by Mark Borgerding: http://sourceforge.net/projects/kissfft/
 *
 * KISS FFT LICENSE INFORMATION
 * ============================
 *
 * Copyright (c) 2003-2013 Mark Borgerding
 *
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * * Redistributions of source code must retain the above copyright notice,
 *   this list of conditions and the following disclaimer.
 * * Redistributions in binary form must reproduce the above copyright notice,
 *   this list of conditions and the following disclaimer in the documentation
 *   and/or other materials provided with the distribution.
 * * Neither the author nor the names of any contributors may be used to
 *   endorse or promote products derived from this software without specific
 *   prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 */
 
#include <pcl/Defs.h>
 
#include <pcl/Complex.h>
#include <pcl/Vector.h>
 
namespace pcl
{
 
// ----------------------------------------------------------------------------
 
class FFT1DBase
{
protected:
 
   static int   OptimizedLength( int, fcomplex* );
   static int   OptimizedLength( int, dcomplex* );
   static int   OptimizedLength( int, float* );
   static int   OptimizedLength( int, double* );
 
   static void* Create( int, fcomplex* );
   static void* Create( int, dcomplex* );
   static void* Create( int, float* );
   static void* Create( int, double* );
 
   static void* CreateInv( int, fcomplex* );
   static void* CreateInv( int, dcomplex* );
   static void* CreateInv( int, float* );
   static void* CreateInv( int, double* );
 
   static void  Destroy( void* );
 
   static void  Transform( void*, fcomplex*, const fcomplex* );
   static void  Transform( void*, dcomplex*, const dcomplex* );
   static void  Transform( void*, fcomplex*, const float* );
   static void  Transform( void*, dcomplex*, const double* );
   static void  Transform( void*, float*,    const fcomplex* );
   static void  Transform( void*, double*,   const dcomplex* );
};
 
// ----------------------------------------------------------------------------
 
#define PCL_FFT_FORWARD    -1
 
#define PCL_FFT_BACKWARD   +1
 
// ----------------------------------------------------------------------------
 
template <typename T>
class PCL_CLASS AbstractFFT : public FFT1DBase
{
public:
 
   using scalar = T;
 
   using complex = Complex<T>;
 
   using vector = GenericVector<scalar>;
 
   using complex_vector = GenericVector<complex>;
 
   using transform = complex_vector;
 
   AbstractFFT( int length )
      : m_length( length )
   {
   }
 
   virtual ~AbstractFFT()
   {
      Release();
   }
 
   int Length() const
   {
      return m_length;
   }
 
   transform DFT() const
   {
      return m_dft;
   }
 
   transform& DFT()
   {
      return m_dft;
   }
 
   const complex* operator *() const
   {
      return *m_dft;
   }
 
   complex* operator *()
   {
      return *m_dft;
   }
 
   virtual void Release()
   {
      if ( m_handle != nullptr )
         this->Destroy( m_handle ), m_handle = nullptr;
      if ( m_handleInv != nullptr )
         this->Destroy( m_handleInv ), m_handleInv = nullptr;
      m_dft = transform();
   }
 
private:
 
   int m_length;
 
protected:
 
   mutable void*     m_handle = nullptr; // Opaque pointers to internal control structures
   mutable void*     m_handleInv = nullptr;
           transform m_dft; // DFT of complex or real data
};
 
// ----------------------------------------------------------------------------
 
#define m_handle    this->m_handle
#define m_handleInv this->m_handleInv
#define m_dft       this->m_dft
#define m_length    this->Length()
 
// ----------------------------------------------------------------------------
 
template <typename T>
class PCL_CLASS GenericFFT : public AbstractFFT<T>
{
public:
 
   using base = AbstractFFT<T>;
 
   using scalar = typename base::scalar;
 
   using complex = typename base::complex;
 
   using vector = typename base::vector;
 
   using complex_vector = typename base::complex_vector;
 
   using transform = typename base::transform;
 
   GenericFFT( int n )
      : base( n )
   {
   }
 
   ~GenericFFT() override
   {
   }
 
   GenericFFT& operator <<( const complex* x )
   {
      if ( m_dft.IsEmpty() )
         m_dft = transform( m_length );
      if ( m_handle == nullptr )
         m_handle = this->Create( m_length, static_cast<complex*>( nullptr ) );
      this->Transform( m_handle, *m_dft, x );
      return *this;
   }
 
   GenericFFT& operator >>( complex* y ) const
   {
      if ( m_dft.IsEmpty() )
         throw Error( "Invalid out-of-place inverse FFT: No FFT has been performed." );
      if ( m_handleInv == nullptr )
         m_handleInv = this->CreateInv( m_length, static_cast<complex*>( nullptr ) );
      this->Transform( m_handleInv, y, *m_dft );
      return *this;
   }
 
   GenericFFT& operator <<( const complex_vector& x )
   {
      PCL_PRECONDITION( x.Length() >= m_length )
      if ( x.Length() < m_length )
         throw Error( "Invalid FFT input vector length." );
      return operator <<( *x );
   }
 
   GenericFFT& operator >>( complex_vector& y ) const
   {
      if ( y.Length() < m_length )
         throw Error( "Invalid FFT output vector length." );
      return operator >>( *y );
   }
 
   GenericFFT& operator()( complex* y, const complex* x, int dir = PCL_FFT_FORWARD ) const
   {
      if ( dir == PCL_FFT_BACKWARD )
      {
         if ( m_handleInv == nullptr )
            m_handleInv = this->CreateInv( m_length, static_cast<complex*>( nullptr ) );
         this->Transform( m_handleInv, y, x );
      }
      else
      {
         if ( m_handle == nullptr )
            m_handle = this->Create( m_length, static_cast<complex*>( nullptr ) );
         this->Transform( m_handle, y, x );
      }
      return const_cast<GenericFFT&>( *this );
   }
 
   static int OptimizedLength( int n )
   {
      return FFT1DBase::OptimizedLength( n, static_cast<complex*>( nullptr ) );
   }
};
 
// ----------------------------------------------------------------------------
 
template <typename T>
class PCL_CLASS GenericRealFFT : public AbstractFFT<T>
{
public:
 
   using base = AbstractFFT<T>;
 
   using scalar = typename base::scalar;
 
   using complex = typename base::complex;
 
   using vector = typename base::vector;
 
   using complex_vector = typename base::complex_vector;
 
   using transform = typename base::transform;
 
   GenericRealFFT( int length )
      : base( length )
   {
      PCL_PRECONDITION( (length & 1) == 0 )
   }
 
   ~GenericRealFFT() override
   {
   }
 
   GenericRealFFT& operator <<( const scalar* x )
   {
      if ( m_dft.IsEmpty() )
         m_dft = transform( m_length/2 + 1 );
      if ( m_handle == nullptr )
         m_handle = this->Create( m_length, static_cast<scalar*>( nullptr ) );
      this->Transform( m_handle, *m_dft, x );
      return *this;
   }
 
   GenericRealFFT& operator >>( scalar* y ) const
   {
      if ( m_dft.IsEmpty() )
         throw Error( "Invalid out-of-place inverse FFT: No FFT has been performed." );
      if ( m_handleInv == nullptr )
         m_handleInv = this->CreateInv( m_length, static_cast<scalar*>( nullptr ) );
      this->Transform( m_handleInv, y, *m_dft );
      return *this;
   }
 
   GenericRealFFT& operator <<( const vector& x )
   {
      if ( x.Length() < m_length )
         throw Error( "Invalid FFT input vector length." );
      return operator <<( *x );
   }
 
   GenericRealFFT& operator >>( vector& y ) const
   {
      if ( y.Length() < m_length )
         throw Error( "Invalid FFT output vector length." );
      return operator >>( *y );
   }
 
   GenericRealFFT& operator()( complex* y, const scalar* x ) const
   {
      if ( m_handle == nullptr )
         m_handle = this->Create( m_length, static_cast<scalar*>( nullptr ) );
      this->Transform( m_handle, y, x );
      return const_cast<GenericRealFFT&>( *this );
   }
 
   GenericRealFFT& operator()( scalar* y, const complex* x ) const
   {
      if ( m_handleInv == nullptr )
         m_handleInv = this->CreateInv( m_length, static_cast<scalar*>( nullptr ) );
      this->Transform( m_handleInv, y, x );
      return const_cast<GenericRealFFT&>( *this );
   }
 
   static int OptimizedLength( int n )
   {
      return FFT1DBase::OptimizedLength( n, static_cast<scalar*>( nullptr ) );
   }
};
 
// ----------------------------------------------------------------------------
 
#undef m_handle
#undef m_handleInv
#undef m_dft
#undef m_length
 
// ----------------------------------------------------------------------------
 
#ifndef __PCL_NO_FFT1D_INSTANTIATE
 
using FFFT = GenericFFT<float>;
 
using DFFT = GenericFFT<double>;
 
using FRealFFT = GenericRealFFT<float>;
 
using DRealFFT = GenericRealFFT<double>;
 
using FFT = FFFT;
 
using RealFFT = FRealFFT;
 
#endif // __PCL_NO_FFT1D_INSTANTIATE
 
// ----------------------------------------------------------------------------
 
} // pcl
 
#endif   // __PCL_FFT1D_h
 
// ----------------------------------------------------------------------------
// EOF pcl/FFT1D.h - Released 2024-06-18T15:48:54Z