conndetect_features.cpp

Go to the documentation of this file.

#include "conndetect_features.h"
#include <iostream>
#include <libraries/math_neon/math_neon.h>
 
using namespace std;
const float number_pi = M_PI;
const float sq2 = sqrtf_neon(2.0f);
 
#define max(a, b) (((a) > (b) ? a : b))
 
static void createHannWindow(Channel &ch) {
  for (us i = 0; i < BLOCKSIZE; i++) {
    float tmp = sinf_neon(number_pi * i / (BLOCKSIZE - 1));
    ch[i] = tmp * tmp;
  }
}
 
Features::Features(const float fs) {
  createHannWindow(_window);
  fft.setup(BLOCKSIZE);
  // Compute frequency array, single sided spectrum
  for (us i = 0; i < BLOCKSIZE / 2 + 1; i++) {
 
    // Dynamically allocated, as this is the input required for the Fft class.
    _wbd.resize(BLOCKSIZE);
    _freq[i] = ((float)i) * fs / ((float)BLOCKSIZE);
  }
}
 
void Features::setChannel(const Channel &channel) {
  _rawdata = &channel;
  static_assert(sizeof(Channel) / sizeof(float) == BLOCKSIZE, "Oops");
  for (us i = 0; i < BLOCKSIZE; i++) {
    _wbd[i] = channel[i] * _window[i];
    _signal_squared[i] = _wbd[i] * _wbd[i];
  }
  _ms_computed = false;
  _fft_computed = false;
  _delta_computed = false;
}
 
float Features::peak_dB() {
  float peak = 0.0f;
 
  // Operate on signal squared, as then we do not have to compute the abs of
  // the data anymore.
  for (us i = 0; i < _wbd.size(); i++) {
    peak = (_signal_squared[i] > peak) ? _signal_squared[i] : peak;
  }
  return 10 * log10f_neon(peak);
}
float Features::rms_dB() {
  compute_ms();
  // 10*log10(_ms) = 20*log10(sqrt(_ms))
  return 10 * log10f_neon(_ms);
}
float Features::kurtosis() {
  compute_ms();
 
  float kurtosis = 0;
  for(us i = 0; i < BLOCKSIZE; i++) {
    kurtosis += _signal_squared[i] * _signal_squared[i];
  }
  kurtosis /= (BLOCKSIZE * _ms);
 
  return kurtosis;
}
 
void Features::compute_ms() {
  if (_ms_computed)
    return;
  float ms = 0.0f;
  for (us i = 0; i < BLOCKSIZE; i++) {
    ms += _signal_squared[i];
  }
  ms /= BLOCKSIZE;
  _ms = ms;
  _ms_computed = true;
}
void Features::compute_fft() {
 
  if (_fft_computed)
    return;
  fft.fft(_wbd);
 
  _Ssq_tot = 0.0f;
  _fc = 0.0f;
 
  float curpeak = 0.0f;
  us idx_peak = 0;
 
  for (us i = 0; i < BLOCKSIZE / 2 + 1; i++) {
    _Ssq[i] = fft.fdr(i) * fft.fdr(i) + fft.fdi(i) * fft.fdi(i);
    _Ssq_tot += _Ssq[i];
    _fc += _Ssq[i] * _freq[i];
    if (curpeak < _Ssq[i]) {
      idx_peak = i;
      curpeak = _Ssq[i];
    }
  }
  _fc /= _Ssq_tot;
  _idx_peak = idx_peak;
 
  _fft_computed = true;
}
float Features::peak_freq() {
  compute_fft();
  return _freq[_idx_peak];
}
 
void Features::compute_delta() {
  if (_delta_computed)
    return;
  for (us i = 0; i < BLOCKSIZE - 1; i++) {
    _absdelta[i] = fabsf_neon(_wbd[i + 1] - _wbd[i]);
  }
  _absdelta[BLOCKSIZE - 1] = 0;
  _delta_computed = true;
}
float Features::wamp(float wampl) {
  us wamp_ctr = 0;
  compute_delta();
 
  for (us i = 0; i < BLOCKSIZE; i++) {
    if (_absdelta[i] > wampl)
      wamp_ctr++;
  }
  return (float)wamp_ctr;
}
float Features::zc(float zcl) {
  us zc_ctr = 0;
  compute_delta();
 
  for (us i = 0; i < BLOCKSIZE - 1; i++) {
    if (_wbd[i] * _wbd[i + 1] < 0 && _absdelta[i] > zcl) {
      zc_ctr++;
    }
  }
  return (float)zc_ctr;
}
float Features::TM3() {
  float tmp = 0.0f;
 
  for (us i = 0; i < BLOCKSIZE; i++) {
    tmp += _signal_squared[i] * (*_rawdata)[i];
  }
  return tmp / BLOCKSIZE;
}
float Features::TM5() {
  float tmp = 0.0f;
  for (us i = 0; i < BLOCKSIZE; i++) {
    tmp += _signal_squared[i] * _signal_squared[i] * (*_rawdata)[i];
  }
  return tmp / BLOCKSIZE;
}
float Features::inverse_variance() {
  compute_fft();
 
  float denom = 0.0f;
  float num = _Ssq_tot * _fc * _fc;
 
  for (us i = 0; i < BLOCKSIZE / 2 + 1; i++) {
    denom += _Ssq[i] * (_freq[i] - _fc) * (_freq[i] - _fc);
  }
 
  return num / denom;
}