Source

jackpipe-dev / src / FftSineGenerator.h

Full commit
/*
 * Copyright © 2011 Aleksey Kunitskiy <alexey.kv@gmail.com>
 * 
 * This file is part of JackPipe
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
#ifndef _FftSineGenerator_H_
#define _FftSineGenerator_H_

#include <algorithm>
#include <memory>
#include <stdexcept>
#include "Vector.h"
#include "LinearSineGenerator.h"
#include "FftwWrapper.h"

namespace JackPipe {

template<typename T>
class FftSineGenerator : public LinearSineGenerator<T> {

  typedef SignalGenerator<T> sg_t;

  FftwWrapper<T> _FftW;

  typename Vector<T>::vec16c_t _fvDftIn;
  typename Vector<T>::vec16_t _fvDftOut;

  void packPhaseShifts(typename Vector<T>::vec16c_t &p);
  void packPhaseShifts(typename Vector<T>::vec16c_t &p,
      typename Vector<T>::vec16_t const& amp);

public:

  FftSineGenerator(const T f0, const T f1,
      const T fSr, const unsigned uB,
      bool bFreqBase2 = false, const long lRs = 0, bool bRoundUp = true)
    throw(std::invalid_argument);
  virtual ~FftSineGenerator() {}

  virtual void generateSignal(typename Vector<T>::vec16_t &v, bool normalize = false) throw(std::runtime_error);
  virtual void generateSignal(typename Vector<T>::vec16_t &v,
      typename Vector<T>::vec16_t const& amp, bool normalize = false) throw(std::runtime_error);

};
  
template<typename T>
FftSineGenerator<T>::FftSineGenerator(
    const T f0, const T f1,
    const T fSr, const unsigned uB,
    bool bFreqBase2, const long lRs, bool bRoundUp)
  throw(std::invalid_argument)
  : LinearSineGenerator<T>(f0, f1, fSr, uB, bFreqBase2, lRs, bRoundUp)
{
  // TODO: as signal generation is't done too often this
  // can be implemented as in-place transform with phase shifts packing
  // on each generateSignal call

  int fftwOpts = FftwWrapper<T>::destroyInput
                  | FftwWrapper<T>::estimate ; // FIXME: 

  _FftW.init_1d(_fvDftIn, _fvDftOut, sg_t::_uBlockSize, fftwOpts );
}

template<typename T>
void FftSineGenerator<T>::packPhaseShifts(typename Vector<T>::vec16c_t &p)
{
  T normFactor = _FftW.getNormFactor();
  for(unsigned i = 0; i < LinearSineGenerator<T>::_uFreqCount; i++)
  {
    //p[0] is the DC component and must be 0
    p[i+1] = std::complex<T>(
        std::sin(LinearSineGenerator<T>::_fvPhaseShifts[i]) / normFactor,
        -std::cos(LinearSineGenerator<T>::_fvPhaseShifts[i]) / normFactor );
  }
}

//FIXME: reuse code from generateSignal(std::vector, bool)
template<typename T>
void FftSineGenerator<T>::packPhaseShifts(
    typename Vector<T>::vec16c_t &p,
    typename Vector<T>::vec16_t const &amp )
{
  T normFactor = _FftW.getNormFactor();
  for(unsigned i = 0; i < LinearSineGenerator<T>::_uFreqCount; i++)
  {
    //p[0] is the DC component and must be 0
    p[i+1] = std::complex<T>(
        amp[i] * std::sin(LinearSineGenerator<T>::_fvPhaseShifts[i]) / normFactor,
        -amp[i] * std::cos(LinearSineGenerator<T>::_fvPhaseShifts[i]) / normFactor );
  }
}

template<typename T>
void FftSineGenerator<T>::generateSignal(
    typename Vector<T>::vec16_t &v, bool normalize)
throw(std::runtime_error)
{
  packPhaseShifts(_fvDftIn);
  _FftW.execute();

  Vector<T>::resizeIfLess(v, FftSineGenerator<T>::getBufferSize());

  if (normalize)
  {
    T max = static_cast<T>(0.0);
    for(unsigned i = 0; i < LinearSineGenerator<T>::_uBlockSize; i++)
    {
      T a = std::abs(_fvDftOut[i]);
      if(a > max)
        max = a;
    }
    Vector<T>::scaleVecDiv(_fvDftOut, max);
  }

  growVector(_fvDftOut, v, LinearSineGenerator<T>::_uBlockCount);
}

template<typename T>
void FftSineGenerator<T>::generateSignal(
    typename Vector<T>::vec16_t &v,
    typename Vector<T>::vec16_t const &amp,
    bool normalize)
throw(std::runtime_error)
{
  packPhaseShifts(_fvDftIn, amp);
  _FftW.execute();

  Vector<T>::resizeIfLess(v, FftSineGenerator<T>::getBufferSize());

  if (normalize)
  {
    T max = static_cast<T>(0.0);
    for(unsigned i = 0; i < LinearSineGenerator<T>::_uBlockSize; i++)
    {
      T a = std::abs(_fvDftOut[i]);
      if(a > max)
        max = a;
    }
    Vector<T>::scaleVecDiv(_fvDftOut, max);
  }

  growVector(_fvDftOut, v, LinearSineGenerator<T>::_uBlockCount);
}

} // JackPipe

#endif
//  vim:set ts=2 sw=2 tw=78 et: