Merge pull request #77 from sp-nitech/goertzel

Add goertzel

Author: takenori-y

.gitignore

@@ -21,6 +21,7 @@ tools/**/
 
 # misc
 *~
+*.bak
 *.swp
 \#*#
 .\#*

CMakeLists.txt

@@ -60,6 +60,7 @@ set(CC_SOURCES
     ${SOURCE_DIR}/analysis/aperiodicity_extraction_by_world.cc
     ${SOURCE_DIR}/analysis/autocorrelation_analysis.cc
     ${SOURCE_DIR}/analysis/fast_fourier_transform_cepstral_analysis.cc
+    ${SOURCE_DIR}/analysis/goertzel_analysis.cc
     ${SOURCE_DIR}/analysis/mel_cepstral_analysis.cc
     ${SOURCE_DIR}/analysis/mel_filter_bank_analysis.cc
     ${SOURCE_DIR}/analysis/mel_frequency_cepstral_coefficients_analysis.cc
@@ -287,6 +288,7 @@ set(MAIN_SOURCES
     ${SOURCE_DIR}/main/gmm.cc
     ${SOURCE_DIR}/main/gmmp.cc
     ${SOURCE_DIR}/main/gnorm.cc
+    ${SOURCE_DIR}/main/goertzel.cc
     ${SOURCE_DIR}/main/grpdelay.cc
     ${SOURCE_DIR}/main/histogram.cc
     ${SOURCE_DIR}/main/huffman.cc

doc/Doxyfile

@@ -0,0 +1,13 @@
+.. _goertzel:
+
+goertzel
+========
+
+.. doxygenfile:: goertzel.cc
+
+.. seealso::
+
+    :ref:`fftr`
+
+.. doxygenclass:: sptk::GoertzelAnalysis
+    :members:

doc/main/goertzel.rst

@@ -0,0 +1,86 @@
+// ------------------------------------------------------------------------ //
+// Copyright 2021 SPTK Working Group                                        //
+//                                                                          //
+// Licensed under the Apache License, Version 2.0 (the "License");          //
+// you may not use this file except in compliance with the License.         //
+// You may obtain a copy of the License at                                  //
+//                                                                          //
+//     http://www.apache.org/licenses/LICENSE-2.0                           //
+//                                                                          //
+// Unless required by applicable law or agreed to in writing, software      //
+// distributed under the License is distributed on an "AS IS" BASIS,        //
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. //
+// See the License for the specific language governing permissions and      //
+// limitations under the License.                                           //
+// ------------------------------------------------------------------------ //
+
+#ifndef SPTK_ANALYSIS_GOERTZEL_ANALYSIS_H_
+#define SPTK_ANALYSIS_GOERTZEL_ANALYSIS_H_
+
+#include <vector>  // std::vector
+
+#include "SPTK/utils/sptk_utils.h"
+
+namespace sptk {
+
+/**
+ * Calculate DFT values at specified frequencies using the Goertzel algorithm.
+ *
+ * The input is the @f$L@f$-length waveform signals:
+ * @f[
+ *   \begin{array}{cccc}
+ *     x(0), & x(1), & \ldots, & x(L-1).
+ *   \end{array}
+ * @f]
+ * The outputs are the real and imaginary parts of the DFT values:
+ * @f[
+ *   \begin{array}{cccc}
+ *   \mathrm{Re}(X(0)), & \mathrm{Re}(X(1)), & \ldots, & \mathrm{Re}(X(K-1)), \\
+ *   \mathrm{Im}(X(0)), & \mathrm{Im}(X(1)), & \ldots, & \mathrm{Im}(X(K-1)),
+ *   \end{array}
+ * @f]
+ * where @f$K@f$ is the number of frequencies to be analyzed.
+ */
+class GoertzelAnalysis {
+ public:
+  /**
+   * @param[in] sampling_rate Sampling rate in Hz.
+   * @param[in] frequencies @f$K@f$ frequencies in Hz to be analyzed.
+   * @param[in] fft_length Number of points assumed in DFT. This is used to
+   *            determine frequency bin resolution.
+   */
+  GoertzelAnalysis(double sampling_rate, const std::vector<double>& frequencies,
+                   int fft_length);
+
+  virtual ~GoertzelAnalysis() {
+  }
+
+  /**
+   * @return True if this object is valid.
+   */
+  bool IsValid() const {
+    return is_valid_;
+  }
+
+  /**
+   * @param[in] signals @f$L@f$-length waveform signals.
+   * @param[out] real_part_output @f$K@f$ real parts.
+   * @param[out] imag_part_output @f$K@f$ imaginary parts.
+   * @return True on success, false on failure.
+   */
+  bool Run(const std::vector<double>& signals,
+           std::vector<double>* real_part_output,
+           std::vector<double>* imag_part_output) const;
+
+ private:
+  bool is_valid_;
+
+  std::vector<double> sine_table_;
+  std::vector<double> cosine_table_;
+
+  DISALLOW_COPY_AND_ASSIGN(GoertzelAnalysis);
+};
+
+}  // namespace sptk
+
+#endif  // SPTK_ANALYSIS_GOERTZEL_ANALYSIS_H_

include/SPTK/analysis/goertzel_analysis.h

@@ -0,0 +1,91 @@
+// ------------------------------------------------------------------------ //
+// Copyright 2021 SPTK Working Group                                        //
+//                                                                          //
+// Licensed under the Apache License, Version 2.0 (the "License");          //
+// you may not use this file except in compliance with the License.         //
+// You may obtain a copy of the License at                                  //
+//                                                                          //
+//     http://www.apache.org/licenses/LICENSE-2.0                           //
+//                                                                          //
+// Unless required by applicable law or agreed to in writing, software      //
+// distributed under the License is distributed on an "AS IS" BASIS,        //
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. //
+// See the License for the specific language governing permissions and      //
+// limitations under the License.                                           //
+// ------------------------------------------------------------------------ //
+
+#include "SPTK/analysis/goertzel_analysis.h"
+
+#include <cmath>    // std::round
+#include <cstddef>  // std::size_t
+#include <vector>   // std::vector
+
+namespace sptk {
+
+GoertzelAnalysis::GoertzelAnalysis(double sampling_rate,
+                                   const std::vector<double>& frequencies,
+                                   int fft_length)
+    : is_valid_(true) {
+  if (sampling_rate <= 0.0 || frequencies.empty() || fft_length <= 0) {
+    is_valid_ = false;
+    return;
+  }
+
+  const int table_size(frequencies.size());
+  sine_table_.resize(table_size);
+  cosine_table_.resize(table_size);
+  for (int k(0); k < table_size; ++k) {
+    if (frequencies[k] < 0.0 || 0.5 * sampling_rate <= frequencies[k]) {
+      is_valid_ = false;
+      break;
+    }
+    const double n(std::round(fft_length * frequencies[k] / sampling_rate));
+    const double omega(sptk::kTwoPi * n / fft_length);
+    sine_table_[k] = std::sin(omega);
+    cosine_table_[k] = std::cos(omega);
+  }
+}
+
+bool GoertzelAnalysis::Run(const std::vector<double>& signals,
+                           std::vector<double>* real_part,
+                           std::vector<double>* imag_part) const {
+  // Check inputs.
+  if (!is_valid_ || NULL == real_part || NULL == imag_part) {
+    return false;
+  }
+
+  // Prepare memories.
+  const int table_size(sine_table_.size());
+  if (real_part->size() != static_cast<std::size_t>(table_size)) {
+    real_part->resize(table_size);
+  }
+  if (imag_part->size() != static_cast<std::size_t>(table_size)) {
+    imag_part->resize(table_size);
+  }
+
+  const int signal_size(signals.size());
+  const double* x(&(signals[0]));
+  double* real(&((*real_part)[0]));
+  double* imag(&((*imag_part)[0]));
+
+  for (int k(0); k < table_size; ++k) {
+    // Perform second-order IIR filtering.
+    const double c(2.0 * cosine_table_[k]);
+    double s1(0.0), s2(0.0);
+    {
+      double s0;
+      for (int t(0); t < signal_size; ++t) {
+        s0 = x[t] + c * s1 - s2;
+        s2 = s1;
+        s1 = s0;
+      }
+    }
+    // Perform first-order FIR filtering.
+    real[k] = s1 * cosine_table_[k] - s2;
+    imag[k] = s1 * sine_table_[k];
+  }
+
+  return true;
+}
+
+}  // namespace sptk