spectrum.h

Go to the documentation of this file.

// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.
#ifndef SP_SPECTRUM_H
#define SP_SPECTRUM_H
namespace sp
{

template <class T1>
arma::cx_vec spectrum(const arma::Col<T1> &x, const arma::vec &W)
{
  arma::cx_vec Pxx(x.size());
  double wc = sum(W);          // Window correction factor
  Pxx       = fft(x % W) / wc; // FFT calc
  return Pxx;
}

template <class T1> arma::vec psd(const arma::Col<T1> &x, const arma::vec &W)
{
  arma::cx_vec X(x.size());
  arma::vec Pxx(x.size());
  X   = spectrum(x, W);    // FFT calc
  Pxx = real(X % conj(X)); // Calc power spectra
  return Pxx;
}

template <class T1> arma::vec psd(const arma::Col<T1> &x)
{
  arma::vec W;
  W = hamming(x.size());
  return psd(x, W);
}

template <class T1>
arma::cx_mat specgram_cx(const arma::Col<T1> &x, const arma::uword Nfft = 512,
                         const arma::uword Noverl = 256)
{
  arma::cx_mat Pw;

  // Def params
  arma::uword N = x.size();
  arma::uword D = Nfft - Noverl;
  arma::uword m = 0;
  if (N > Nfft)
  {
    arma::Col<T1> xk(Nfft);
    arma::vec W(Nfft);

    W             = hamming(Nfft);
    arma::uword U = static_cast<arma::uword>(floor((N - Noverl) / double(D)));
    Pw.set_size(Nfft, U);
    Pw.zeros();

    // Avg loop
    for (arma::uword k = 0; k < N - Nfft; k += D)
    {
      xk          = x.rows(k, k + Nfft - 1); // Pick out chunk
      Pw.col(m++) = spectrum(xk, W);         // Calculate spectrum
    }
  }
  else
  {
    arma::vec W(N);
    W = hamming(N);
    Pw.set_size(N, 1);
    Pw = spectrum(x, W);
  }
  return Pw;
}

template <class T1>
arma::mat specgram(const arma::Col<T1> &x, const arma::uword Nfft = 512,
                   const arma::uword Noverl = 256)
{
  arma::cx_mat Pw;
  arma::mat Sg;
  Pw = specgram_cx(x, Nfft, Noverl);
  Sg = real(Pw % conj(Pw)); // Calculate power spectrum
  return Sg;
}

template <class T1>
arma::mat specgram_ph(const arma::Col<T1> &x, const arma::uword Nfft = 512,
                      const arma::uword Noverl = 256)
{
  arma::cx_mat Pw;
  arma::mat Sg;
  Pw = specgram_cx(x, Nfft, Noverl);
  Sg = angle(Pw); // Calculate phase spectrum
  return Sg;
}

template <class T1>
arma::vec pwelch_ph(const arma::Col<T1> &x, const arma::uword Nfft = 512,
                    const arma::uword Noverl = 256)
{
  arma::mat Ph;
  Ph = specgram_ph(x, Nfft, Noverl);
  return arma::mean(Ph, 1);
}

template <class T1>
arma::vec pwelch(const arma::Col<T1> &x, const arma::uword Nfft = 512,
                 const arma::uword Noverl = 256)
{
  arma::mat Pxx;
  Pxx = specgram(x, Nfft, Noverl);
  return arma::mean(Pxx, 1);
}

template <class T1>
std::complex<double> goertzel(const arma::Col<T1> &x, const double f)
{
  // Constants
  arma::uword N = x.size();
  double Q      = f / N;
  double A      = PI_2 * Q;
  double B      = 2 * cos(A);
  std::complex<double> C(cos(A), -sin(A));
  // States
  T1 s0 = 0;
  T1 s1 = 0;
  T1 s2 = 0;

  // Accumulate data
  for (arma::uword n = 0; n < N; n++)
  {
    // Update filter
    s0 = x(n) + B * s1 - s2;

    // Shift buffer
    s2 = s1;
    s1 = s0;
  }
  // Update output state
  s0 = B * s1 - s2;

  // Return the complex DFT output
  return s0 - s1 * C;
}

template <class T1>
arma::cx_vec goertzel(const arma::Col<T1> &x, const arma::vec f)
{
  arma::uword N = f.size();
  arma::cx_vec P(N);
  for (arma::uword n = 0; n < N; n++)
  {
    P(n) = goertzel(x, f(n));
  }

  // Return the complex DFT output vector
  return P;
}

} // namespace sp
#endif