// 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"
#include "../math.h"
#include "../fft/fft.h"
#include "../fft/fft_stl.h"
namespace dlib
{
// ----------------------------------------------------------------------------------------
template < typename T, typename Alloc >
matrix<std::complex<T>,0,1> fft (const std::vector<std::complex<T>, Alloc>& in)
{
//complex FFT
static_assert(std::is_floating_point<T>::value, "only support floating point types");
matrix<std::complex<T>,0,1> out(in.size());
if (in.size() != 0)
fft({(long)in.size()}, &in[0], &out(0,0), false);
return out;
}
// ----------------------------------------------------------------------------------------
template < typename T, long NR, long NC, typename MM, typename L >
matrix<std::complex<T>,NR,NC,MM,L> fft (const matrix<std::complex<T>,NR,NC,MM,L>& in)
{
//complex FFT
static_assert(std::is_floating_point<T>::value, "only support floating point types");
matrix<std::complex<T>,NR,NC,MM,L> out(in.nr(), in.nc());
if (in.size() != 0)
fft({in.nr(),in.nc()}, &in(0,0), &out(0,0), false);
return out;
}
// ----------------------------------------------------------------------------------------
template <typename EXP>
typename EXP::matrix_type fft (const matrix_exp<EXP>& data)
{
//complex FFT for expression template
static_assert(is_complex<typename EXP::type>::value, "input should be complex");
typename EXP::matrix_type in(data);
return fft(in);
}
// ----------------------------------------------------------------------------------------
template < typename T, typename Alloc >
matrix<std::complex<T>,0,1> ifft (const std::vector<std::complex<T>, Alloc>& in)
{
//complex FFT
static_assert(std::is_floating_point<T>::value, "only support floating point types");
matrix<std::complex<T>,0,1> out(in.size());
if (in.size() != 0)
{
fft({(long)in.size()}, &in[0], &out(0,0), true);
out /= out.size();
}
return out;
}
// ----------------------------------------------------------------------------------------
template < typename T, long NR, long NC, typename MM, typename L >
matrix<std::complex<T>,NR,NC,MM,L> ifft (const matrix<std::complex<T>,NR,NC,MM,L>& in)
{
//inverse complex FFT
static_assert(std::is_floating_point<T>::value, "only support floating point types");
matrix<std::complex<T>,NR,NC,MM,L> out(in.nr(), in.nc());
if (in.size() != 0)
{
fft({in.nr(),in.nc()}, &in(0,0), &out(0,0), true);
out /= out.size();
}
return out;
}
// ----------------------------------------------------------------------------------------
template <typename EXP>
typename EXP::matrix_type ifft (const matrix_exp<EXP>& data)
{
//inverse complex FFT for expression template
static_assert(is_complex<typename EXP::type>::value, "input should be complex");
typename EXP::matrix_type in(data);
return ifft(in);
}
// ----------------------------------------------------------------------------------------
template<typename T, long NR, long NC, typename MM, typename L>
matrix<std::complex<T>,NR,fftr_nc_size(NC),MM,L> fftr (const matrix<T,NR,NC,MM,L>& in)
{
//real FFT
static_assert(std::is_floating_point<T>::value, "only support floating point types");
DLIB_ASSERT(in.nc() % 2 == 0, "last dimension " << in.nc() << " needs to be even otherwise ifftr(fftr(data)) won't have matching dimensions");
matrix<std::complex<T>,NR,fftr_nc_size(NC),MM,L> out(in.nr(), fftr_nc_size(in.nc()));
if (in.size() != 0)
fftr({in.nr(),in.nc()}, &in(0,0), &out(0,0));
return out;
}
// ----------------------------------------------------------------------------------------
template <typename EXP>
matrix<add_complex_t<typename EXP::type>> fftr (const matrix_exp<EXP>& data)
{
//real FFT for expression template
static_assert(std::is_floating_point<typename EXP::type>::value, "input should be real");
matrix<typename EXP::type> in(data);
return fftr(in);
}
// ----------------------------------------------------------------------------------------
template<typename T, long NR, long NC, typename MM, typename L>
matrix<T,NR,ifftr_nc_size(NC),MM,L> ifftr (const matrix<std::complex<T>,NR,NC,MM,L>& in)
{
//inverse real FFT
static_assert(std::is_floating_point<T>::value, "only support floating point types");
matrix<T,NR,ifftr_nc_size(NC),MM,L> out(in.nr(), ifftr_nc_size(in.nc()));
if (in.size() != 0)
{
ifftr({out.nr(),out.nc()}, &in(0,0), &out(0,0));
out /= out.size();
}
return out;
}
// ----------------------------------------------------------------------------------------
template <typename EXP>
matrix<remove_complex_t<typename EXP::type>> ifftr (const matrix_exp<EXP>& data)
{
//inverse real FFT for expression template
static_assert(is_complex<typename EXP::type>::value, "input should be complex");
matrix<typename EXP::type> in(data);
return ifftr(in);
}
// ----------------------------------------------------------------------------------------
template < typename T, long NR, long NC, typename MM, typename L >
void fft_inplace (matrix<std::complex<T>,NR,NC,MM,L>& data)
{
static_assert(std::is_floating_point<T>::value, "only support floating point types");
if (data.size() != 0)
fft({data.nr(),data.nc()}, &data(0,0), &data(0,0), false);
}
// ----------------------------------------------------------------------------------------
template < typename T, long NR, long NC, typename MM, typename L >
void ifft_inplace (matrix<std::complex<T>,NR,NC,MM,L>& data)
{
static_assert(std::is_floating_point<T>::value, "only support floating point types");
if (data.size() != 0)
fft({data.nr(),data.nc()}, &data(0,0), &data(0,0), true);
}
// ----------------------------------------------------------------------------------------
namespace details
{
struct fft_func
{
template<typename MAT, typename T = typename MAT::type, typename std::enable_if<is_complex<T>::value, bool>::type = true>
auto operator()(const MAT& mat) const { return dlib::fft(mat); }
template<typename MAT, typename T = typename MAT::type, typename std::enable_if<!is_complex<T>::value, bool>::type = true>
auto operator()(const MAT& mat) const { return dlib::fft(dlib::complex_matrix(mat)); }
static constexpr std::size_t freqsize(std::size_t fftsize) { return fftsize; }
};
struct fftr_func
{
template<typename MAT>
auto operator()(const MAT& mat) const { return dlib::fftr(mat); }
static constexpr std::size_t freqsize(std::size_t fftsize) { return dlib::fftr_nc_size(fftsize); }
};
struct ifft_func
{
template<typename MAT>
auto operator()(const MAT& mat) const { return dlib::ifft(mat); }
};
struct ifftr_func
{
template<typename MAT>
auto operator()(const MAT& mat) const { return dlib::ifftr(mat); }
};
template <
typename EXP,
typename WINDOW,
typename FFT_FUNC
>
auto stft_impl (
const matrix_exp<EXP>& signal,
const WINDOW& w,
std::size_t fftsize,
std::size_t wlen,
std::size_t hoplen,
const FFT_FUNC& fft_obj
)
{
using T = typename EXP::type;
using R = remove_complex_t<T>;
using C = add_complex_t<T>;
static_assert(std::is_floating_point<R>::value, "underlying type must be real or complex floating point type");
DLIB_ASSERT(is_vector(signal), "input must be a vector type");
DLIB_ASSERT(signal.size() >= (long)wlen, "signal.size() >= wlen not satisfied");
DLIB_ASSERT(fftsize >= wlen, "fftsize >= wlen not satisfied");
DLIB_ASSERT(wlen >= hoplen, "wlen >= hoplen not satisfied");
// Input is left-padded by wlen/2 and right-padded wlen/2
const std::size_t total_padding = wlen;
const std::size_t overlap = wlen - hoplen;
const std::size_t nframes = (signal.size() + total_padding - overlap) / hoplen;
matrix<C> stft = zeros_matrix<C>(nframes,
FFT_FUNC::freqsize(fftsize));
matrix<R> win(1,wlen);
for (std::size_t i = 0 ; i < wlen ; ++i)
win(0, i) = w(i, wlen);
// TODO: reduce extra buffers, e.g. padded
matrix<T> padded;
if (is_row_vector(signal))
padded = join_rows(join_rows(zeros_matrix<T>(1, wlen/2), signal), zeros_matrix<T>(1, wlen/2));
else
padded = join_rows(join_rows(zeros_matrix<T>(1, wlen/2), trans(signal)), zeros_matrix<T>(1, wlen/2));
for (long i = 0 ; i < stft.nr() ; ++i)
{
set_rowm(stft, i) = fft_obj(join_rows(pointwise_multiply(win, subm(padded, 0, i*hoplen, 1, wlen)),
zeros_matrix<T>(1, fftsize - wlen)));
}
return stft;
}
template <
typename ReturnType,
typename EXP,
typename WINDOW,
typename IFFT_FUNC
>
auto istft_impl (
const matrix_exp<EXP>& stft,
const WINDOW& w,
std::size_t wlen,
std::size_t hoplen,
const IFFT_FUNC& ifft_obj
)
{
using T = typename EXP::type;
using R = remove_complex_t<T>;
static_assert(is_complex<T>::value, "matrix type must be complex");
static_assert(std::is_floating_point<R>::value, "underlying type must be complex floating point type");
DLIB_ASSERT(stft.nc() > 0 && stft.nr() > 0, "stft must be non-empty");
DLIB_ASSERT(ifftr_nc_size(stft.nc()) >= (long)wlen, "fftsize >= wlen not satisfied");
DLIB_ASSERT(wlen >= hoplen, "wlen >= hoplen not satisfied");
const size_t ntime = (stft.nr() - 1) * hoplen + wlen;
matrix<ReturnType> signal = zeros_matrix<ReturnType>(1, ntime);
matrix<R> norm = zeros_matrix<R>(1, ntime);
matrix<R> win(1, wlen);
for (std::size_t i = 0 ; i < wlen ; ++i)
win(0, i) = w(i, wlen);
matrix<R> win2 = squared(win);
for (long t = 0 ; t < stft.nr() ; ++t)
{
set_subm(signal, 0, t*hoplen, 1, wlen) += pointwise_multiply(win, subm(ifft_obj(rowm(stft, t)), 0, 0, 1, wlen));
set_subm(norm, 0, t*hoplen, 1, wlen) += win2;
}
// Remove padding of wlen/2 and wlen/2 on either end
DLIB_ASSERT(sum(subm(norm, 0, wlen/2, 1, ntime - wlen) < 1e-13) == 0, "NOLA constraint not satisfied");
signal = pointwise_divide(subm(signal, 0, wlen/2, 1, ntime - wlen),
subm(norm, 0, wlen/2, 1, ntime - wlen));
return signal;
}
}
// ----------------------------------------------------------------------------------------
inline auto make_hann()
{
return [](std::size_t i, std::size_t N) {return hann(i, N, PERIODIC); };
}
inline auto make_blackman()
{
return [](std::size_t i, std::size_t N) {return blackman(i, N, PERIODIC);};
}
inline auto make_blackman_nuttall()
{
return [](std::size_t i, std::size_t N) {return blackman_nuttall(i, N, PERIODIC);};
}
inline auto make_blackman_harris()
{
return [](std::size_t i, std::size_t N) { return blackman_harris(i, N, PERIODIC); };
}
inline auto make_blackman_harris7()
{
return [](std::size_t i, std::size_t N) { return blackman_harris7(i, N, PERIODIC); };
}
inline auto make_kaiser(beta_t beta)
{
return [=](std::size_t i, std::size_t N){return kaiser(i, N, beta, PERIODIC);};
}
// ----------------------------------------------------------------------------------------
template <typename EXP, typename WINDOW>
auto stft (
const matrix_exp<EXP>& signal,
const WINDOW& w,
std::size_t fftsize,
std::size_t wlen,
std::size_t hoplen
)
{
return details::stft_impl(signal, w, fftsize, wlen, hoplen, details::fft_func{});
}
// ----------------------------------------------------------------------------------------
template <typename T, typename Alloc, typename WINDOW>
auto stft (
const std::vector<T, Alloc>& signal,
const WINDOW& w,
std::size_t fftsize,
std::size_t wlen,
std::size_t hoplen
)
{
return stft(dlib::mat(signal), w, fftsize, wlen, hoplen);
}
// ----------------------------------------------------------------------------------------
template <typename EXP,typename WINDOW>
auto istft (
const matrix_exp<EXP>& stft,
const WINDOW& w,
std::size_t wlen,
std::size_t hoplen
)
{
using T = typename EXP::type;
return details::istft_impl<T>(stft, w, wlen, hoplen, details::ifft_func{});
}
// ----------------------------------------------------------------------------------------
template <typename EXP, typename WINDOW>
auto stftr (
const matrix_exp<EXP>& signal,
const WINDOW& w,
std::size_t fftsize,
std::size_t wlen,
std::size_t hoplen
)
{
return details::stft_impl(signal, w, fftsize, wlen, hoplen, details::fftr_func{});
}
// ----------------------------------------------------------------------------------------
template <typename T, typename Alloc, typename WINDOW>
auto stftr (
const std::vector<T, Alloc>& signal,
const WINDOW& w,
std::size_t fftsize,
std::size_t wlen,
std::size_t hoplen
)
{
return stftr(dlib::mat(signal), w, fftsize, wlen, hoplen);
}
// ----------------------------------------------------------------------------------------
template <typename EXP, typename WINDOW>
auto istftr (
const matrix_exp<EXP>& stft,
const WINDOW& w,
std::size_t wlen,
std::size_t hoplen
)
{
using R = remove_complex_t<typename EXP::type>;
return details::istft_impl<R>(stft, w, wlen, hoplen, details::ifftr_func{});
}
// ----------------------------------------------------------------------------------------
}
#endif // DLIB_FFt_Hh_