vtkFFT.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
// SPDX-License-Identifier: BSD-3-Clause
#ifndef vtkFFT_h
#define vtkFFT_h
 
#include "vtkAOSDataArrayTemplate.h" // vtkAOSDataArrayTemplate
#include "vtkCommonMathModule.h"     // export macro
#include "vtkDataArray.h"            // vtkDataArray
#include "vtkDataArrayRange.h"       // vtkDataArrayRange
#include "vtkMath.h"                 // vtkMath::Pi
#include "vtkObject.h"
#include "vtkSMPTools.h" // vtkSMPTools::Transform, vtkSMPTools::For
 
#include "vtk_kissfft.h" // kiss_fft_scalar, kiss_fft_cpx
// clang-format off
#include VTK_KISSFFT_HEADER(kiss_fft.h)
#include VTK_KISSFFT_HEADER(tools/kiss_fftr.h)
// clang-format on
 
#include <array>       // std::array
#include <cmath>       // std::sin, std::cos, std::sqrt
#include <numeric>     // std::accumulate, std::inner_product
#include <type_traits> // std::enable_if, std::iterator_traits
#include <vector>      // std::vector
 
VTK_ABI_NAMESPACE_BEGIN
class VTKCOMMONMATH_EXPORT vtkFFT : public vtkObject
{
public:
 
  using ScalarNumber = kiss_fft_scalar;
  using ComplexNumber = kiss_fft_cpx;
  static_assert(sizeof(ComplexNumber) == 2 * sizeof(ScalarNumber),
    "vtkFFT::ComplexNumber definition is not valid");
 
  using vtkScalarNumberArray = vtkAOSDataArrayTemplate<vtkFFT::ScalarNumber>;
  template <typename T>
  struct isFftType : public std::false_type
  {
  };
 
  enum Octave
  {
    Hz_31_5 = 5,
    Hz_63 = 6,
    Hz_125 = 7,
    Hz_250 = 8,
    Hz_500 = 9,
    kHz_1 = 10,
    kHz_2 = 11,
    kHz_4 = 12,
    kHz_8 = 13,
    kHz_16 = 14
  };
 
  enum OctaveSubdivision
  {
    Full,
    FirstHalf,
    SecondHalf,
    FirstThird,
    SecondThird,
    ThirdThird
  };
 
 
  static std::vector<ComplexNumber> Fft(const std::vector<ScalarNumber>& in);
  static void Fft(ScalarNumber* input, std::size_t size, ComplexNumber* result);
  static std::vector<ComplexNumber> Fft(const std::vector<ComplexNumber>& in);
  static void Fft(ComplexNumber* input, std::size_t size, ComplexNumber* result);
#ifndef __VTK_WRAP__
  static vtkSmartPointer<vtkScalarNumberArray> Fft(vtkScalarNumberArray* input);
#endif
 
 
  static std::vector<ComplexNumber> RFft(const std::vector<ScalarNumber>& in);
  static void RFft(ScalarNumber* input, std::size_t size, ComplexNumber* result);
#ifndef __VTK_WRAP__
  static vtkSmartPointer<vtkScalarNumberArray> RFft(vtkScalarNumberArray* input);
#endif
 
  static std::vector<ComplexNumber> IFft(const std::vector<ComplexNumber>& in);
 
  static std::vector<ScalarNumber> IRFft(const std::vector<ComplexNumber>& in);
 
  static inline ScalarNumber Abs(const ComplexNumber& in);
 
  static inline ScalarNumber SquaredAbs(const ComplexNumber& in);
 
  static inline ComplexNumber Conjugate(const ComplexNumber& in);
 
  static std::vector<ScalarNumber> FftFreq(int windowLength, double sampleSpacing);
 
  static std::vector<ScalarNumber> RFftFreq(int windowLength, double sampleSpacing);
 
  static std::array<double, 2> GetOctaveFrequencyRange(Octave octave,
    OctaveSubdivision octaveSubdivision = OctaveSubdivision::Full, bool baseTwo = true);
 
 
#ifndef __VTK_WRAP__
  template <typename T, typename TW, typename std::enable_if<isFftType<T>::value>::type* = nullptr>
  static std::vector<ComplexNumber> OverlappingFft(const std::vector<T>& signal,
    const std::vector<TW>& window, std::size_t noverlap, bool detrend, bool onesided,
    unsigned int* shape = nullptr);
 
  template <typename TW>
  static vtkFFT::ComplexNumber* OverlappingFft(vtkFFT::vtkScalarNumberArray* signal,
    const std::vector<TW>& window, std::size_t noverlap, bool detrend, bool onesided,
    unsigned int* shape = nullptr);
#endif
 
  enum Scaling : int
  {
    Density = 0, 
    Spectrum     
  };
 
  enum SpectralMode : int
  {
    STFT = 0, 
    PSD       
  };
 
 
#ifndef __VTK_WRAP__
  template <typename T, typename TW, typename std::enable_if<isFftType<T>::value>::type* = nullptr>
  static std::vector<ComplexNumber> Spectrogram(const std::vector<T>& signal,
    const std::vector<TW>& window, double sampleRate, int noverlap, bool detrend, bool onesided,
    vtkFFT::Scaling scaling, vtkFFT::SpectralMode mode, unsigned int* shape = nullptr,
    bool transpose = false);
 
  template <typename TW>
  static vtkSmartPointer<vtkFFT::vtkScalarNumberArray> Spectrogram(
    vtkFFT::vtkScalarNumberArray* signal, const std::vector<TW>& window, double sampleRate,
    int noverlap, bool detrend, bool onesided, vtkFFT::Scaling scaling, vtkFFT::SpectralMode mode,
    unsigned int* shape = nullptr, bool transpose = false);
#endif
 
 
#ifndef __VTK_WRAP__
  template <typename T, typename TW, typename std::enable_if<isFftType<T>::value>::type* = nullptr>
  static std::vector<vtkFFT::ScalarNumber> Csd(const std::vector<T>& signal,
    const std::vector<TW>& window, double sampleRate, int noverlap, bool detrend, bool onesided,
    vtkFFT::Scaling scaling);
 
  template <typename TW>
  static vtkSmartPointer<vtkFFT::vtkScalarNumberArray> Csd(vtkScalarNumberArray* signal,
    const std::vector<TW>& window, double sampleRate, int noverlap, bool detrend, bool onesided,
    vtkFFT::Scaling scaling);
#endif
 
#ifndef __VTK_WRAP__
  template <typename T>
  static void Transpose(T* data, unsigned int* shape);
#endif
 
 
  using WindowGenerator = ScalarNumber (*)(std::size_t, std::size_t);
 
  static inline ScalarNumber HanningGenerator(std::size_t x, std::size_t size);
  static inline ScalarNumber BartlettGenerator(std::size_t x, std::size_t size);
  static inline ScalarNumber SineGenerator(std::size_t x, std::size_t size);
  static inline ScalarNumber BlackmanGenerator(std::size_t x, std::size_t size);
  static inline ScalarNumber RectangularGenerator(std::size_t x, std::size_t size);
 
  template <typename T>
  static void GenerateKernel1D(T* kernel, std::size_t n, WindowGenerator generator);
 
  template <typename T>
  static void GenerateKernel2D(T* kernel, std::size_t n, std::size_t m, WindowGenerator generator);
 
  static vtkFFT* New();
  vtkTypeMacro(vtkFFT, vtkObject);
  void PrintSelf(ostream& os, vtkIndent indent) override;
 
protected:
  vtkFFT() = default;
  ~vtkFFT() override = default;
 
  template <typename T>
  constexpr static T Zero();
 
#ifndef __VTK_WRAP__
  template <typename InputIt>
  static typename std::iterator_traits<InputIt>::value_type ComputeScaling(
    InputIt begin, InputIt end, Scaling scaling, double fs);
 
  template <typename T, typename TW>
  static void PreprocessAndDispatchFft(const T* segment, const std::vector<TW>& window,
    bool detrend, bool onesided, vtkFFT::ComplexNumber* result);
 
  static void RFft(ComplexNumber* input, std::size_t size, ComplexNumber* result);
 
  template <typename TW>
  static void ScaleFft(ComplexNumber* fft, unsigned int shape[2], const std::vector<TW>& window,
    double sampleRate, bool onesided, vtkFFT::Scaling scaling, vtkFFT::SpectralMode mode);
#endif
 
private:
  vtkFFT(const vtkFFT&) = delete;
  void operator=(const vtkFFT&) = delete;
};
 
//------------------------------------------------------------------------------
template <>
struct vtkFFT::isFftType<vtkFFT::ScalarNumber> : public std::true_type
{
};
template <>
struct vtkFFT::isFftType<vtkFFT::ComplexNumber> : public std::true_type
{
};
 
//------------------------------------------------------------------------------
template <typename T>
constexpr T vtkFFT::Zero()
{
  return static_cast<T>(0);
}
template <>
constexpr vtkFFT::ComplexNumber vtkFFT::Zero()
{
  return vtkFFT::ComplexNumber{ 0.0, 0.0 };
}
 
//------------------------------------------------------------------------------
inline vtkFFT::ComplexNumber operator+(
  const vtkFFT::ComplexNumber& lhs, const vtkFFT::ComplexNumber& rhs)
{
  return vtkFFT::ComplexNumber{ lhs.r + rhs.r, lhs.i + rhs.i };
}
inline vtkFFT::ComplexNumber operator-(
  const vtkFFT::ComplexNumber& lhs, const vtkFFT::ComplexNumber& rhs)
{
  return vtkFFT::ComplexNumber{ lhs.r - rhs.r, lhs.i - rhs.i };
}
inline vtkFFT::ComplexNumber operator*(
  const vtkFFT::ComplexNumber& lhs, const vtkFFT::ComplexNumber& rhs)
{
  return vtkFFT::ComplexNumber{ (lhs.r * rhs.r) - (lhs.i * rhs.i),
    (lhs.r * rhs.i) + (lhs.i * rhs.r) };
}
inline vtkFFT::ComplexNumber operator*(
  const vtkFFT::ComplexNumber& lhs, const vtkFFT::ScalarNumber& rhs)
{
  return vtkFFT::ComplexNumber{ lhs.r * rhs, lhs.i * rhs };
}
inline vtkFFT::ComplexNumber operator/(
  const vtkFFT::ComplexNumber& lhs, const vtkFFT::ComplexNumber& rhs)
{
  const double divisor = rhs.r * rhs.r + rhs.i * rhs.i;
  return vtkFFT::ComplexNumber{ ((lhs.r * rhs.r) + (lhs.i * rhs.i)) / divisor,
    ((lhs.i * rhs.r) - (lhs.r * rhs.i)) / divisor };
}
inline vtkFFT::ComplexNumber operator/(
  const vtkFFT::ComplexNumber& lhs, const vtkFFT::ScalarNumber& rhs)
{
  return vtkFFT::ComplexNumber{ lhs.r / rhs, lhs.i / rhs };
}
 
//------------------------------------------------------------------------------
vtkFFT::ScalarNumber vtkFFT::Abs(const ComplexNumber& in)
{
  return std::sqrt(in.r * in.r + in.i * in.i);
}
 
//------------------------------------------------------------------------------
vtkFFT::ScalarNumber vtkFFT::SquaredAbs(const ComplexNumber& in)
{
  return in.r * in.r + in.i * in.i;
}
 
//------------------------------------------------------------------------------
vtkFFT::ComplexNumber vtkFFT::Conjugate(const ComplexNumber& in)
{
  return ComplexNumber{ in.r, -in.i };
}
 
//------------------------------------------------------------------------------
double vtkFFT::HanningGenerator(std::size_t x, std::size_t size)
{
  return 0.5 * (1.0 - std::cos(2.0 * vtkMath::Pi() * x / (size - 1)));
}
 
//------------------------------------------------------------------------------
double vtkFFT::BartlettGenerator(std::size_t x, std::size_t size)
{
  return 2.0 * x / (size - 1);
}
 
//------------------------------------------------------------------------------
double vtkFFT::SineGenerator(std::size_t x, std::size_t size)
{
  return std::sin(vtkMath::Pi() * x / (size - 1));
}
 
//------------------------------------------------------------------------------
double vtkFFT::BlackmanGenerator(std::size_t x, std::size_t size)
{
  const double cosin = std::cos((2.0 * vtkMath::Pi() * x) / (size - 1));
  return 0.42 - 0.5 * cosin + 0.08 * (2.0 * cosin * cosin - 1.0);
}
 
//------------------------------------------------------------------------------
double vtkFFT::RectangularGenerator(std::size_t, std::size_t)
{
  return 1.0;
}
 
//------------------------------------------------------------------------------
template <typename T>
void vtkFFT::GenerateKernel1D(T* kernel, std::size_t n, WindowGenerator generator)
{
  std::size_t half = (n / 2) + (n % 2);
  for (std::size_t i = 0; i < half; ++i)
  {
    kernel[i] = kernel[n - 1 - i] = generator(i, n);
  }
}
 
//------------------------------------------------------------------------------
template <typename T>
void vtkFFT::GenerateKernel2D(T* kernel, std::size_t n, std::size_t m, WindowGenerator generator)
{
  const std::size_t halfX = (n / 2) + (n % 2);
  const std::size_t halfY = (m / 2) + (m % 2);
  for (std::size_t i = 0; i < halfX; ++i)
  {
    for (std::size_t j = 0; j < halfY; ++j)
    {
      // clang-format off
      kernel[i][j]
      = kernel[n - 1 - i][j]
      = kernel[i][m - 1 - j]
      = kernel[n - 1 - i][m - 1 - j]
      = generator(i, n) * generator(j, m);
      // clang-format on
    }
  }
}
 
VTK_ABI_NAMESPACE_END
 
#include "vtkFFT.txx" // complex templated functions not wrapped by python
 
#endif // vtkFFT_h