70 lines
2.6 KiB
C++
70 lines
2.6 KiB
C++
#include "fftconv.hpp"
|
|
#include "fft.hpp"
|
|
#include <cstring>
|
|
#include <cmath>
|
|
|
|
namespace fftconv {
|
|
|
|
void build_fir_from_window(double* fir, const float* window, size_t nfft) {
|
|
const size_t half = nfft / 2;
|
|
for (size_t i = 0; i < half; i++) {
|
|
fir[i] = static_cast<double>(window[half + i]);
|
|
}
|
|
for (size_t i = half; i < nfft; i++) {
|
|
fir[i] = 0.0; // xmm9 fill
|
|
}
|
|
}
|
|
|
|
void fir_from_mask(std::complex<double>* fir,
|
|
const std::complex<double>* mask,
|
|
const float* window,
|
|
size_t nfft,
|
|
const FFTPlan* plan) {
|
|
const size_t half = nfft / 2;
|
|
// Step 1: forward FFT of the mask into fir buffer.
|
|
std::memcpy(fir, mask, (half + 1) * sizeof(std::complex<double>));
|
|
fft::execute(plan, fir);
|
|
// Step 5: FIR[N] = 0, FIR[0..N/2-1] = window[N/2..N-1].
|
|
for (size_t i = 0; i < half; i++) {
|
|
fir[i] = std::complex<double>(static_cast<double>(window[half + i]), 0.0);
|
|
}
|
|
for (size_t i = half; i < nfft; i++) {
|
|
fir[i] = std::complex<double>(0.0, 0.0);
|
|
}
|
|
// Step 6b: inverse FFT -> time-domain FIR.
|
|
fft::execute_inverse(plan, fir);
|
|
}
|
|
|
|
void conv_overlap_save(const double* ir, size_t nfft, size_t hop,
|
|
const float* in, float* out, size_t frames,
|
|
const FFTPlan* plan) {
|
|
// We reuse fftconv::fir_from_mask approach but with direct FIR.
|
|
// overlap-save: process block of size nfft, keep tail of hop samples.
|
|
// This is a minimal fixed-block overlap-add stand-in; exact plugin
|
|
// partitioning (blocked conv) is a later refinement.
|
|
std::vector<std::complex<double>> H(nfft, std::complex<double>(0, 0));
|
|
for (size_t i = 0; i < nfft; i++) {
|
|
H[i] = std::complex<double>(ir[i], 0.0);
|
|
}
|
|
fft::execute(plan, H.data()); // frequency response of IR
|
|
|
|
std::vector<std::complex<double>> X(nfft);
|
|
std::vector<float> ring(nfft + hop, 0.0f);
|
|
|
|
for (size_t n = 0; n < frames; n += hop) {
|
|
// shift ring
|
|
std::memmove(ring.data(), ring.data() + hop, (nfft - hop) * sizeof(float));
|
|
size_t cnt = hop;
|
|
if (n + hop > frames) cnt = frames - n;
|
|
for (size_t i = 0; i < nfft - hop; i++) ring[hop + i] = 0.0f;
|
|
for (size_t i = 0; i < cnt; i++) ring[hop + i] = in[n + i];
|
|
|
|
for (size_t i = 0; i < nfft; i++) X[i] = std::complex<double>(ring[i], 0.0);
|
|
fft::execute(plan, X.data());
|
|
for (size_t i = 0; i < nfft; i++) X[i] *= H[i];
|
|
fft::execute_inverse(plan, X.data());
|
|
for (size_t i = 0; i < cnt; i++) out[n + i] = static_cast<float>(X[i].real());
|
|
}
|
|
}
|
|
|
|
} // namespace fftconv
|