mirror of
https://gitcode.com/gh_mirrors/ope/OpenFace.git
synced 2026-07-17 18:47:48 +00:00
279 lines
9.1 KiB
C++
279 lines
9.1 KiB
C++
// Copyright (C) 2013 Davis E. King (davis@dlib.net)
|
|
// License: Boost Software License See LICENSE.txt for the full license.
|
|
#ifndef DLIB_FFt_Hh_
|
|
#define DLIB_FFt_Hh_
|
|
|
|
#include "matrix_fft_abstract.h"
|
|
#include "matrix_utilities.h"
|
|
#include "../hash.h"
|
|
#include "../algs.h"
|
|
|
|
#ifdef DLIB_USE_FFTW
|
|
#include <fftw3.h>
|
|
#endif // DLIB_USE_FFTW
|
|
|
|
namespace dlib
|
|
{
|
|
|
|
// ----------------------------------------------------------------------------------------
|
|
|
|
namespace impl
|
|
{
|
|
inline unsigned long reverse_bits (
|
|
unsigned long val,
|
|
unsigned long num
|
|
)
|
|
{
|
|
unsigned long temp = 0;
|
|
for (unsigned long i = 0; i < num; ++i)
|
|
{
|
|
temp <<= 1;
|
|
temp |= val&0x1;
|
|
val >>= 1;
|
|
}
|
|
return temp;
|
|
}
|
|
|
|
template <typename EXP>
|
|
void permute (
|
|
const matrix_exp<EXP>& data,
|
|
typename EXP::matrix_type& outdata
|
|
)
|
|
{
|
|
outdata.set_size(data.size());
|
|
if (data.size() == 0)
|
|
return;
|
|
|
|
const unsigned long num = static_cast<unsigned long>(std::log((double)data.size())/std::log(2.0) + 0.5);
|
|
for (unsigned long i = 0; i < (unsigned long)data.size(); ++i)
|
|
{
|
|
outdata(impl::reverse_bits(i,num)) = data(i);
|
|
}
|
|
}
|
|
|
|
}
|
|
|
|
// ----------------------------------------------------------------------------------------
|
|
|
|
inline bool is_power_of_two (
|
|
const unsigned long& value
|
|
)
|
|
{
|
|
if (value == 0)
|
|
return true;
|
|
else
|
|
return count_bits(value) == 1;
|
|
}
|
|
|
|
// ----------------------------------------------------------------------------------------
|
|
|
|
template <typename EXP>
|
|
typename EXP::matrix_type fft (
|
|
const matrix_exp<EXP>& data
|
|
)
|
|
{
|
|
if (data.size() == 0)
|
|
return data;
|
|
|
|
// You have to give a complex matrix
|
|
COMPILE_TIME_ASSERT(is_complex<typename EXP::type>::value);
|
|
// make sure requires clause is not broken
|
|
DLIB_CASSERT(is_vector(data) && is_power_of_two(data.size()),
|
|
"\t void ifft(data)"
|
|
<< "\n\t data must be a vector with a size that is a power of two."
|
|
<< "\n\t is_vector(data): " << is_vector(data)
|
|
<< "\n\t data.size(): " << data.size()
|
|
);
|
|
|
|
typedef typename EXP::type::value_type T;
|
|
|
|
typename EXP::matrix_type outdata(data);
|
|
|
|
const long half = outdata.size()/2;
|
|
|
|
typedef std::complex<T> ct;
|
|
matrix<ct,0,1,typename EXP::mem_manager_type> twiddle_factors(half);
|
|
|
|
// compute the complex root of unity w
|
|
const T temp = -2.0*pi/outdata.size();
|
|
ct w = ct(std::cos(temp),std::sin(temp));
|
|
|
|
ct w_pow = 1;
|
|
|
|
// compute the twiddle factors
|
|
for (long j = 0; j < twiddle_factors.size(); ++j)
|
|
{
|
|
twiddle_factors(j) = w_pow;
|
|
w_pow *= w;
|
|
}
|
|
|
|
|
|
// now compute the decimation in frequency. This first
|
|
// outer loop loops log2(outdata.size()) number of times
|
|
long skip = 1;
|
|
for (long step = half; step != 0; step >>= 1)
|
|
{
|
|
// do blocks of butterflies in this loop
|
|
for (long j = 0; j < outdata.size(); j += step*2)
|
|
{
|
|
// do step butterflies
|
|
for (long k = 0; k < step; ++k)
|
|
{
|
|
const long a_idx = j+k;
|
|
const long b_idx = j+k+step;
|
|
const ct a = outdata(a_idx) + outdata(b_idx);
|
|
const ct b = (outdata(a_idx) - outdata(b_idx))*twiddle_factors(k*skip);
|
|
outdata(a_idx) = a;
|
|
outdata(b_idx) = b;
|
|
}
|
|
}
|
|
skip *= 2;
|
|
}
|
|
|
|
typename EXP::matrix_type outperm;
|
|
impl::permute(outdata, outperm);
|
|
return outperm;
|
|
}
|
|
|
|
// ----------------------------------------------------------------------------------------
|
|
|
|
template <typename EXP>
|
|
typename EXP::matrix_type ifft (
|
|
const matrix_exp<EXP>& data
|
|
)
|
|
{
|
|
if (data.size() == 0)
|
|
return data;
|
|
|
|
// You have to give a complex matrix
|
|
COMPILE_TIME_ASSERT(is_complex<typename EXP::type>::value);
|
|
// make sure requires clause is not broken
|
|
DLIB_CASSERT(is_vector(data) && is_power_of_two(data.size()),
|
|
"\t void ifft(data)"
|
|
<< "\n\t data must be a vector with a size that is a power of two."
|
|
<< "\n\t is_vector(data): " << is_vector(data)
|
|
<< "\n\t data.size(): " << data.size()
|
|
);
|
|
|
|
|
|
typedef typename EXP::type::value_type T;
|
|
|
|
typename EXP::matrix_type outdata;
|
|
impl::permute(data,outdata);
|
|
|
|
const long half = outdata.size()/2;
|
|
|
|
typedef std::complex<T> ct;
|
|
matrix<ct,0,1,typename EXP::mem_manager_type> twiddle_factors(half);
|
|
|
|
// compute the complex root of unity w
|
|
const T temp = 2.0*pi/outdata.size();
|
|
ct w = ct(std::cos(temp),std::sin(temp));
|
|
|
|
ct w_pow = 1;
|
|
|
|
// compute the twiddle factors
|
|
for (long j = 0; j < twiddle_factors.size(); ++j)
|
|
{
|
|
twiddle_factors(j) = w_pow;
|
|
w_pow *= w;
|
|
}
|
|
|
|
// now compute the inverse decimation in frequency. This first
|
|
// outer loop loops log2(outdata.size()) number of times
|
|
long skip = half;
|
|
for (long step = 1; step <= half; step <<= 1)
|
|
{
|
|
// do blocks of butterflies in this loop
|
|
for (long j = 0; j < outdata.size(); j += step*2)
|
|
{
|
|
// do step butterflies
|
|
for (long k = 0; k < step; ++k)
|
|
{
|
|
const long a_idx = j+k;
|
|
const long b_idx = j+k+step;
|
|
outdata(b_idx) *= twiddle_factors(k*skip);
|
|
const ct a = outdata(a_idx) + outdata(b_idx);
|
|
const ct b = outdata(a_idx) - outdata(b_idx);
|
|
outdata(a_idx) = a;
|
|
outdata(b_idx) = b;
|
|
}
|
|
}
|
|
skip /= 2;
|
|
}
|
|
|
|
outdata /= outdata.size();
|
|
return outdata;
|
|
}
|
|
|
|
// ----------------------------------------------------------------------------------------
|
|
|
|
#ifdef DLIB_USE_FFTW
|
|
|
|
template <long NR, long NC, typename MM, typename L>
|
|
matrix<std::complex<double>,NR,NC,MM,L> call_fftw_fft(
|
|
const matrix<std::complex<double>,NR,NC,MM,L>& data
|
|
)
|
|
{
|
|
// make sure requires clause is not broken
|
|
DLIB_CASSERT(is_vector(data) && is_power_of_two(data.size()),
|
|
"\t void fft(data)"
|
|
<< "\n\t data must be a vector with a size that is a power of two."
|
|
<< "\n\t is_vector(data): " << is_vector(data)
|
|
<< "\n\t data.size(): " << data.size()
|
|
);
|
|
|
|
matrix<std::complex<double>,NR,NC,MM,L> m2(data.nr(),data.nc());
|
|
fftw_complex *in, *out;
|
|
fftw_plan p;
|
|
in = (fftw_complex*)&data(0);
|
|
out = (fftw_complex*)&m2(0);
|
|
p = fftw_plan_dft_1d(data.size(), in, out, FFTW_FORWARD, FFTW_ESTIMATE);
|
|
fftw_execute(p);
|
|
fftw_destroy_plan(p);
|
|
return m2;
|
|
}
|
|
|
|
template <long NR, long NC, typename MM, typename L>
|
|
matrix<std::complex<double>,NR,NC,MM,L> call_fftw_ifft(
|
|
const matrix<std::complex<double>,NR,NC,MM,L>& data
|
|
)
|
|
{
|
|
// make sure requires clause is not broken
|
|
DLIB_CASSERT(is_vector(data) && is_power_of_two(data.size()),
|
|
"\t void ifft(data)"
|
|
<< "\n\t data must be a vector with a size that is a power of two."
|
|
<< "\n\t is_vector(data): " << is_vector(data)
|
|
<< "\n\t data.size(): " << data.size()
|
|
);
|
|
|
|
matrix<std::complex<double>,NR,NC,MM,L> m2(data.nr(),data.nc());
|
|
fftw_complex *in, *out;
|
|
fftw_plan p;
|
|
in = (fftw_complex*)&data(0);
|
|
out = (fftw_complex*)&m2(0);
|
|
p = fftw_plan_dft_1d(data.size(), in, out, FFTW_BACKWARD, FFTW_ESTIMATE);
|
|
fftw_execute(p);
|
|
fftw_destroy_plan(p);
|
|
return m2/data.size();
|
|
}
|
|
|
|
// ----------------------------------------------------------------------------------------
|
|
|
|
// call FFTW for these cases:
|
|
inline matrix<std::complex<double>,0,1> fft (const matrix<std::complex<double>,0,1>& data) {return call_fftw_fft(data);}
|
|
inline matrix<std::complex<double>,0,1> ifft(const matrix<std::complex<double>,0,1>& data) {return call_fftw_ifft(data);}
|
|
inline matrix<std::complex<double>,1,0> fft (const matrix<std::complex<double>,1,0>& data) {return call_fftw_fft(data);}
|
|
inline matrix<std::complex<double>,1,0> ifft(const matrix<std::complex<double>,1,0>& data) {return call_fftw_ifft(data);}
|
|
inline matrix<std::complex<double> > fft (const matrix<std::complex<double> >& data) {return call_fftw_fft(data);}
|
|
inline matrix<std::complex<double> > ifft(const matrix<std::complex<double> >& data) {return call_fftw_ifft(data);}
|
|
|
|
#endif // DLIB_USE_FFTW
|
|
|
|
// ----------------------------------------------------------------------------------------
|
|
|
|
}
|
|
|
|
#endif // DLIB_FFt_Hh_
|
|
|