apply_convolution()

std::vector< double > MayaFlux::Kinesis::Discrete::apply_convolution	(	std::span< const double >	src,
		std::span< const double >	kernel,
		const ConvolutionProcessor &	processor,
		bool	full_size = false
	)

Core FFT convolution engine.

Parameters

src	Input signal
kernel	Kernel signal
processor	Frequency-domain operation
full_size	When true returns linear convolution length (N+M-1); when false returns the same number of samples as src

Returns

Processed output

Definition at line 12 of file Convolution.cpp.

{
    if (src.empty() || kernel.empty())
        return std::vector<double>(src.size(), 0.0);
 
    const size_t conv_len = src.size() + kernel.size() - 1;
    const size_t fft_size = std::bit_ceil(std::max(size_t { 256 }, conv_len));
    const size_t bins = fft_size / 2 + 1;
 
    Eigen::VectorXd padded_src = Eigen::VectorXd::Zero(static_cast<Eigen::Index>(fft_size));
    Eigen::VectorXd padded_kernel = Eigen::VectorXd::Zero(static_cast<Eigen::Index>(fft_size));
    std::ranges::copy(src, padded_src.begin());
    std::ranges::copy(kernel, padded_kernel.begin());
 
    Eigen::FFT<double> fft;
    Eigen::VectorXcd src_fft, ker_fft;
    fft.fwd(src_fft, padded_src);
    fft.fwd(ker_fft, padded_kernel);
 
    std::vector<std::complex<double>> sig(bins), ker(bins), res(bins);
    P::for_each(P::par_unseq,
        std::views::iota(size_t { 0 }, bins).begin(),
        std::views::iota(size_t { 0 }, bins).end(),
        [&](size_t b) {
            sig[b] = src_fft(static_cast<Eigen::Index>(b));
            ker[b] = ker_fft(static_cast<Eigen::Index>(b));
        });
 
    processor(sig, ker, res);
 
    Eigen::VectorXcd full_fft(static_cast<Eigen::Index>(fft_size));
    P::for_each(P::par_unseq,
        std::views::iota(size_t { 0 }, fft_size).begin(),
        std::views::iota(size_t { 0 }, fft_size).end(),
        [&](size_t b) {
            full_fft(static_cast<Eigen::Index>(b)) = (b < bins)
                ? res[b]
                : std::conj(res[fft_size - b]);
        });
 
    const double inv = 1.0 / static_cast<double>(fft_size);
    full_fft *= inv;
 
    Eigen::VectorXd time_result;
    fft.inv(time_result, full_fft);
 
    const size_t out_len = full_size ? conv_len : src.size();
    std::vector<double> out(out_len);
    P::for_each(P::par_unseq,
        std::views::iota(size_t { 0 }, out_len).begin(),
        std::views::iota(size_t { 0 }, out_len).end(),
        [&](size_t i) { out[i] = time_result(static_cast<Eigen::Index>(i)); });
    return out;
}

References b.

Referenced by auto_correlate(), convolve(), cross_correlate(), and deconvolve().

Here is the caller graph for this function: