VTK  9.3.20230924
vtkFFT.h
Go to the documentation of this file.
1 // SPDX-FileCopyrightText: Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
2 // SPDX-License-Identifier: BSD-3-Clause
23 #ifndef vtkFFT_h
24 #define vtkFFT_h
25 
26 #include "vtkAOSDataArrayTemplate.h" // vtkAOSDataArrayTemplate
27 #include "vtkCommonMathModule.h" // export macro
28 #include "vtkDataArray.h" // vtkDataArray
29 #include "vtkDataArrayRange.h" // vtkDataArrayRange
30 #include "vtkMath.h" // vtkMath::Pi
31 #include "vtkObject.h"
32 #include "vtkSMPTools.h" // vtkSMPTools::Transform, vtkSMPTools::For
33 
34 #include "vtk_kissfft.h" // kiss_fft_scalar, kiss_fft_cpx
35 // clang-format off
36 #include VTK_KISSFFT_HEADER(kiss_fft.h)
37 #include VTK_KISSFFT_HEADER(tools/kiss_fftr.h)
38 // clang-format on
39 
40 #include <array> // std::array
41 #include <cmath> // std::sin, std::cos, std::sqrt
42 #include <numeric> // std::accumulate, std::inner_product
43 #include <type_traits> // std::enable_if, std::iterator_traits
44 #include <vector> // std::vector
45 
46 VTK_ABI_NAMESPACE_BEGIN
47 class VTKCOMMONMATH_EXPORT vtkFFT : public vtkObject
48 {
49 public:
51 
73  using ScalarNumber = kiss_fft_scalar;
74  using ComplexNumber = kiss_fft_cpx;
75  static_assert(sizeof(ComplexNumber) == 2 * sizeof(ScalarNumber),
76  "vtkFFT::ComplexNumber definition is not valid");
77 
79  template <typename T>
80  struct isFftType : public std::false_type
81  {
82  };
84 
89  enum Octave
90  {
91  Hz_31_5 = 5,
92  Hz_63 = 6,
93  Hz_125 = 7,
94  Hz_250 = 8,
95  Hz_500 = 9,
96  kHz_1 = 10,
97  kHz_2 = 11,
98  kHz_4 = 12,
99  kHz_8 = 13,
100  kHz_16 = 14
101  };
102 
108  {
114  ThirdThird
115  };
116 
118 
125  static std::vector<ComplexNumber> Fft(const std::vector<ScalarNumber>& in);
126  static void Fft(ScalarNumber* input, std::size_t size, ComplexNumber* result);
127  static std::vector<ComplexNumber> Fft(const std::vector<ComplexNumber>& in);
128  static void Fft(ComplexNumber* input, std::size_t size, ComplexNumber* result);
129 #ifndef __VTK_WRAP__
131 #endif
133 
135 
141  static std::vector<ComplexNumber> RFft(const std::vector<ScalarNumber>& in);
142  static void RFft(ScalarNumber* input, std::size_t size, ComplexNumber* result);
143 #ifndef __VTK_WRAP__
145 #endif
147 
158  static std::vector<ComplexNumber> IFft(const std::vector<ComplexNumber>& in);
159 
168  static std::vector<ScalarNumber> IRFft(const std::vector<ComplexNumber>& in);
169 
173  static inline ScalarNumber Abs(const ComplexNumber& in);
174 
178  static inline ScalarNumber SquaredAbs(const ComplexNumber& in);
179 
183  static inline ComplexNumber Conjugate(const ComplexNumber& in);
184 
188  static std::vector<ScalarNumber> FftFreq(int windowLength, double sampleSpacing);
189 
194  static std::vector<ScalarNumber> RFftFreq(int windowLength, double sampleSpacing);
195 
203  static std::array<double, 2> GetOctaveFrequencyRange(Octave octave,
204  OctaveSubdivision octaveSubdivision = OctaveSubdivision::Full, bool baseTwo = true);
205 
207 
226 #ifndef __VTK_WRAP__
228  static std::vector<ComplexNumber> OverlappingFft(const std::vector<T>& signal,
229  const std::vector<TW>& window, std::size_t noverlap, bool detrend, bool onesided,
230  unsigned int* shape = nullptr);
231 
232  template <typename TW>
234  const std::vector<TW>& window, std::size_t noverlap, bool detrend, bool onesided,
235  unsigned int* shape = nullptr);
236 #endif
238 
242  enum Scaling : int
243  {
244  Density = 0,
245  Spectrum
246  };
247 
251  enum SpectralMode : int
252  {
253  STFT = 0,
254  PSD
255  };
256 
258 
283 #ifndef __VTK_WRAP__
285  static std::vector<ComplexNumber> Spectrogram(const std::vector<T>& signal,
286  const std::vector<TW>& window, double sampleRate, int noverlap, bool detrend, bool onesided,
287  vtkFFT::Scaling scaling, vtkFFT::SpectralMode mode, unsigned int* shape = nullptr,
288  bool transpose = false);
289 
290  template <typename TW>
292  vtkFFT::vtkScalarNumberArray* signal, const std::vector<TW>& window, double sampleRate,
293  int noverlap, bool detrend, bool onesided, vtkFFT::Scaling scaling, vtkFFT::SpectralMode mode,
294  unsigned int* shape = nullptr, bool transpose = false);
295 #endif
297 
299 
320 #ifndef __VTK_WRAP__
322  static std::vector<vtkFFT::ScalarNumber> Csd(const std::vector<T>& signal,
323  const std::vector<TW>& window, double sampleRate, int noverlap, bool detrend, bool onesided,
324  vtkFFT::Scaling scaling);
325 
326  template <typename TW>
328  const std::vector<TW>& window, double sampleRate, int noverlap, bool detrend, bool onesided,
329  vtkFFT::Scaling scaling);
330 #endif
332 
343 #ifndef __VTK_WRAP__
344  template <typename T>
345  static void Transpose(T* data, unsigned int* shape);
346 #endif
347 
349 
358  using WindowGenerator = ScalarNumber (*)(std::size_t, std::size_t);
359 
360  static inline ScalarNumber HanningGenerator(std::size_t x, std::size_t size);
361  static inline ScalarNumber BartlettGenerator(std::size_t x, std::size_t size);
362  static inline ScalarNumber SineGenerator(std::size_t x, std::size_t size);
363  static inline ScalarNumber BlackmanGenerator(std::size_t x, std::size_t size);
364  static inline ScalarNumber RectangularGenerator(std::size_t x, std::size_t size);
366 
371  template <typename T>
372  static void GenerateKernel1D(T* kernel, std::size_t n, WindowGenerator generator);
373 
378  template <typename T>
379  static void GenerateKernel2D(T* kernel, std::size_t n, std::size_t m, WindowGenerator generator);
380 
381  static vtkFFT* New();
382  vtkTypeMacro(vtkFFT, vtkObject);
383  void PrintSelf(ostream& os, vtkIndent indent) override;
384 
385 protected:
386  vtkFFT() = default;
387  ~vtkFFT() override = default;
388 
392  template <typename T>
393  constexpr static T Zero();
394 
395 #ifndef __VTK_WRAP__
400  template <typename InputIt>
401  static typename std::iterator_traits<InputIt>::value_type ComputeScaling(
402  InputIt begin, InputIt end, Scaling scaling, double fs);
403 
408  template <typename T, typename TW>
409  static void PreprocessAndDispatchFft(const T* segment, const std::vector<TW>& window,
410  bool detrend, bool onesided, vtkFFT::ComplexNumber* result);
411 
420  static void RFft(ComplexNumber* input, std::size_t size, ComplexNumber* result);
421 
425  template <typename TW>
426  static void ScaleFft(ComplexNumber* fft, unsigned int shape[2], const std::vector<TW>& window,
427  double sampleRate, bool onesided, vtkFFT::Scaling scaling, vtkFFT::SpectralMode mode);
428 #endif
429 
430 private:
431  vtkFFT(const vtkFFT&) = delete;
432  void operator=(const vtkFFT&) = delete;
433 };
434 
435 //------------------------------------------------------------------------------
436 template <>
437 struct vtkFFT::isFftType<vtkFFT::ScalarNumber> : public std::true_type
438 {
439 };
440 template <>
441 struct vtkFFT::isFftType<vtkFFT::ComplexNumber> : public std::true_type
442 {
443 };
444 
445 //------------------------------------------------------------------------------
446 template <typename T>
447 constexpr T vtkFFT::Zero()
448 {
449  return static_cast<T>(0);
450 }
451 template <>
453 {
454  return vtkFFT::ComplexNumber{ 0.0, 0.0 };
455 }
456 
457 //------------------------------------------------------------------------------
459  const vtkFFT::ComplexNumber& lhs, const vtkFFT::ComplexNumber& rhs)
460 {
461  return vtkFFT::ComplexNumber{ lhs.r + rhs.r, lhs.i + rhs.i };
462 }
464  const vtkFFT::ComplexNumber& lhs, const vtkFFT::ComplexNumber& rhs)
465 {
466  return vtkFFT::ComplexNumber{ lhs.r - rhs.r, lhs.i - rhs.i };
467 }
469  const vtkFFT::ComplexNumber& lhs, const vtkFFT::ComplexNumber& rhs)
470 {
471  return vtkFFT::ComplexNumber{ (lhs.r * rhs.r) - (lhs.i * rhs.i),
472  (lhs.r * rhs.i) + (lhs.i * rhs.r) };
473 }
475  const vtkFFT::ComplexNumber& lhs, const vtkFFT::ScalarNumber& rhs)
476 {
477  return vtkFFT::ComplexNumber{ lhs.r * rhs, lhs.i * rhs };
478 }
480  const vtkFFT::ComplexNumber& lhs, const vtkFFT::ComplexNumber& rhs)
481 {
482  const double divisor = rhs.r * rhs.r + rhs.i * rhs.i;
483  return vtkFFT::ComplexNumber{ ((lhs.r * rhs.r) + (lhs.i * rhs.i)) / divisor,
484  ((lhs.i * rhs.r) - (lhs.r * rhs.i)) / divisor };
485 }
487  const vtkFFT::ComplexNumber& lhs, const vtkFFT::ScalarNumber& rhs)
488 {
489  return vtkFFT::ComplexNumber{ lhs.r / rhs, lhs.i / rhs };
490 }
491 
492 //------------------------------------------------------------------------------
494 {
495  return std::sqrt(in.r * in.r + in.i * in.i);
496 }
497 
498 //------------------------------------------------------------------------------
500 {
501  return in.r * in.r + in.i * in.i;
502 }
503 
504 //------------------------------------------------------------------------------
506 {
507  return ComplexNumber{ in.r, -in.i };
508 }
509 
510 //------------------------------------------------------------------------------
511 double vtkFFT::HanningGenerator(std::size_t x, std::size_t size)
512 {
513  return 0.5 * (1.0 - std::cos(2.0 * vtkMath::Pi() * x / (size - 1)));
514 }
515 
516 //------------------------------------------------------------------------------
517 double vtkFFT::BartlettGenerator(std::size_t x, std::size_t size)
518 {
519  return 2.0 * x / (size - 1);
520 }
521 
522 //------------------------------------------------------------------------------
523 double vtkFFT::SineGenerator(std::size_t x, std::size_t size)
524 {
525  return std::sin(vtkMath::Pi() * x / (size - 1));
526 }
527 
528 //------------------------------------------------------------------------------
529 double vtkFFT::BlackmanGenerator(std::size_t x, std::size_t size)
530 {
531  const double cosin = std::cos((2.0 * vtkMath::Pi() * x) / (size - 1));
532  return 0.42 - 0.5 * cosin + 0.08 * (2.0 * cosin * cosin - 1.0);
533 }
534 
535 //------------------------------------------------------------------------------
536 double vtkFFT::RectangularGenerator(std::size_t, std::size_t)
537 {
538  return 1.0;
539 }
540 
541 //------------------------------------------------------------------------------
542 template <typename T>
543 void vtkFFT::GenerateKernel1D(T* kernel, std::size_t n, WindowGenerator generator)
544 {
545  std::size_t half = (n / 2) + (n % 2);
546  for (std::size_t i = 0; i < half; ++i)
547  {
548  kernel[i] = kernel[n - 1 - i] = generator(i, n);
549  }
550 }
551 
552 //------------------------------------------------------------------------------
553 template <typename T>
554 void vtkFFT::GenerateKernel2D(T* kernel, std::size_t n, std::size_t m, WindowGenerator generator)
555 {
556  const std::size_t halfX = (n / 2) + (n % 2);
557  const std::size_t halfY = (m / 2) + (m % 2);
558  for (std::size_t i = 0; i < halfX; ++i)
559  {
560  for (std::size_t j = 0; j < halfY; ++j)
561  {
562  // clang-format off
563  kernel[i][j]
564  = kernel[n - 1 - i][j]
565  = kernel[i][m - 1 - j]
566  = kernel[n - 1 - i][m - 1 - j]
567  = generator(i, n) * generator(j, m);
568  // clang-format on
569  }
570  }
571 }
572 
573 VTK_ABI_NAMESPACE_END
574 
575 #include "vtkFFT.txx" // complex templated functions not wrapped by python
576 
577 #endif // vtkFFT_h
Array-Of-Structs implementation of vtkGenericDataArray.
perform Discrete Fourier Transforms
Definition: vtkFFT.h:48
static ScalarNumber BlackmanGenerator(std::size_t x, std::size_t size)
Window generator functions.
Definition: vtkFFT.h:529
static ScalarNumber SineGenerator(std::size_t x, std::size_t size)
Window generator functions.
Definition: vtkFFT.h:523
static void Fft(ComplexNumber *input, std::size_t size, ComplexNumber *result)
Compute the one-dimensional DFT for complex input.
ScalarNumber(*)(std::size_t, std::size_t) WindowGenerator
Window generator functions.
Definition: vtkFFT.h:358
SpectralMode
Spectral modes for Spectrogram and Csd functions.
Definition: vtkFFT.h:252
kiss_fft_scalar ScalarNumber
Useful type definitions and utilities.
Definition: vtkFFT.h:73
OctaveSubdivision
Enum specifying which octave band we want to compute.
Definition: vtkFFT.h:108
@ Full
Definition: vtkFFT.h:109
@ FirstThird
Definition: vtkFFT.h:112
@ FirstHalf
Definition: vtkFFT.h:110
@ SecondHalf
Definition: vtkFFT.h:111
@ SecondThird
Definition: vtkFFT.h:113
static std::vector< ComplexNumber > RFft(const std::vector< ScalarNumber > &in)
Compute the one-dimensional DFT for real input.
static ScalarNumber RectangularGenerator(std::size_t x, std::size_t size)
Window generator functions.
Definition: vtkFFT.h:536
static void ScaleFft(ComplexNumber *fft, unsigned int shape[2], const std::vector< TW > &window, double sampleRate, bool onesided, vtkFFT::Scaling scaling, vtkFFT::SpectralMode mode)
Scale a fft according to its window and some mode.
vtkFFT()=default
static std::vector< ScalarNumber > IRFft(const std::vector< ComplexNumber > &in)
Compute the inverse of RFft.
static std::vector< ComplexNumber > Fft(const std::vector< ComplexNumber > &in)
Compute the one-dimensional DFT for complex input.
static ScalarNumber HanningGenerator(std::size_t x, std::size_t size)
Window generator functions.
Definition: vtkFFT.h:511
static void RFft(ComplexNumber *input, std::size_t size, ComplexNumber *result)
XXX(c++17): This function should NOT exist and is here just for the sake template unfolding purposes.
static std::iterator_traits< InputIt >::value_type ComputeScaling(InputIt begin, InputIt end, Scaling scaling, double fs)
For a given window defined by begin and end, compute the scaling needed to apply to the resulting FFT...
static ScalarNumber Abs(const ComplexNumber &in)
Return the absolute value (also known as norm, modulus, or magnitude) of complex number.
Definition: vtkFFT.h:493
static void Transpose(T *data, unsigned int *shape)
Transpose in place an inlined 2D matrix.
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)
Compute consecutive Fourier transforms Welch method without averaging nor scaling the result.
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)
Compute the Cross Spectral Density of a given signal.
static std::vector< ComplexNumber > IFft(const std::vector< ComplexNumber > &in)
Compute the inverse of Fft.
static std::vector< ScalarNumber > RFftFreq(int windowLength, double sampleSpacing)
Return the DFT sample frequencies for the real version of the dft (see Rfft).
static void GenerateKernel2D(T *kernel, std::size_t n, std::size_t m, WindowGenerator generator)
Given a window generator function, create a symmetric 2D kernel.
Definition: vtkFFT.h:554
static ComplexNumber Conjugate(const ComplexNumber &in)
Return the conjugate of the given complex number.
Definition: vtkFFT.h:505
constexpr static T Zero()
Templated zero value, specialized for vtkFFT::ComplexNumber.
static void GenerateKernel1D(T *kernel, std::size_t n, WindowGenerator generator)
Given a window generator function, create a symmetric 1D kernel.
Definition: vtkFFT.h:543
static void RFft(ScalarNumber *input, std::size_t size, ComplexNumber *result)
Compute the one-dimensional DFT for real input.
static void PreprocessAndDispatchFft(const T *segment, const std::vector< TW > &window, bool detrend, bool onesided, vtkFFT::ComplexNumber *result)
Dispatch the signal to the right FFT function according to the given parameters.
static vtkSmartPointer< vtkScalarNumberArray > Fft(vtkScalarNumberArray *input)
Compute the one-dimensional DFT for complex input.
static std::array< double, 2 > GetOctaveFrequencyRange(Octave octave, OctaveSubdivision octaveSubdivision=OctaveSubdivision::Full, bool baseTwo=true)
Return lower and upper frequency from a octave band number / nominal midband frequency.
static ScalarNumber BartlettGenerator(std::size_t x, std::size_t size)
Window generator functions.
Definition: vtkFFT.h:517
~vtkFFT() override=default
Octave
Enum containing octave band numbers, named upon their nominal midband frequency.
Definition: vtkFFT.h:90
static vtkFFT * New()
static ScalarNumber SquaredAbs(const ComplexNumber &in)
Return the squared absolute value of the complex number.
Definition: vtkFFT.h:499
kiss_fft_cpx ComplexNumber
Useful type definitions and utilities.
Definition: vtkFFT.h:74
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)
Compute a spectrogram with consecutive Fourier transforms using Welch method.
static vtkFFT::ComplexNumber * OverlappingFft(vtkFFT::vtkScalarNumberArray *signal, const std::vector< TW > &window, std::size_t noverlap, bool detrend, bool onesided, unsigned int *shape=nullptr)
Compute consecutive Fourier transforms Welch method without averaging nor scaling the result.
static std::vector< ComplexNumber > Fft(const std::vector< ScalarNumber > &in)
Compute the one-dimensional DFT for complex input.
Scaling
Scaling modes for Spectrogram and Csd functions.
Definition: vtkFFT.h:243
static vtkSmartPointer< vtkScalarNumberArray > RFft(vtkScalarNumberArray *input)
Compute the one-dimensional DFT for real input.
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)
Compute a spectrogram with consecutive Fourier transforms using Welch method.
static std::vector< ScalarNumber > FftFreq(int windowLength, double sampleSpacing)
Return the DFT sample frequencies.
void PrintSelf(ostream &os, vtkIndent indent) override
Methods invoked by print to print information about the object including superclasses.
static vtkSmartPointer< vtkFFT::vtkScalarNumberArray > Csd(vtkScalarNumberArray *signal, const std::vector< TW > &window, double sampleRate, int noverlap, bool detrend, bool onesided, vtkFFT::Scaling scaling)
Compute the Cross Spectral Density of a given signal.
static void Fft(ScalarNumber *input, std::size_t size, ComplexNumber *result)
Compute the one-dimensional DFT for complex input.
a simple class to control print indentation
Definition: vtkIndent.h:108
static constexpr double Pi()
A mathematical constant.
Definition: vtkMath.h:197
abstract base class for most VTK objects
Definition: vtkObject.h:161
Hold a reference to a vtkObjectBase instance.
@ mode
Definition: vtkX3D.h:247
@ value
Definition: vtkX3D.h:220
@ sampleRate
Definition: vtkX3D.h:472
@ type
Definition: vtkX3D.h:516
@ size
Definition: vtkX3D.h:253
@ data
Definition: vtkX3D.h:315
STL-compatible iterable ranges that provide access to vtkDataArray elements.
vtkFFT::ComplexNumber operator*(const vtkFFT::ComplexNumber &lhs, const vtkFFT::ComplexNumber &rhs)
Definition: vtkFFT.h:468
vtkFFT::ComplexNumber operator+(const vtkFFT::ComplexNumber &lhs, const vtkFFT::ComplexNumber &rhs)
Definition: vtkFFT.h:458
vtkFFT::ComplexNumber operator-(const vtkFFT::ComplexNumber &lhs, const vtkFFT::ComplexNumber &rhs)
Definition: vtkFFT.h:463
vtkFFT::ComplexNumber operator/(const vtkFFT::ComplexNumber &lhs, const vtkFFT::ComplexNumber &rhs)
Definition: vtkFFT.h:479