54 lines
2.0 KiB
C++
54 lines
2.0 KiB
C++
#include <cstdio>
|
|
#include <cmath>
|
|
#include <cstring>
|
|
#include <vector>
|
|
#include <complex>
|
|
#include "fft.hpp"
|
|
#include "fftconv.hpp"
|
|
#include "tables_data.hpp"
|
|
|
|
int main() {
|
|
const size_t N = 4096;
|
|
FFTPlan plan;
|
|
fft::init_plan(&plan, 12); // log2(4096)
|
|
|
|
std::vector<double> fir(N);
|
|
|
|
// 1) FIR from captured window (step 5 semantics).
|
|
fftconv::build_fir_from_window(fir.data(), WIN_WINDOW, N);
|
|
double esum = 0.0;
|
|
for (size_t i = 0; i < N; i++) esum += fir[i] * fir[i];
|
|
std::printf("step5 FIR: half-sum=%.3f energy=%.3f fir[0]=%.4f fir[2047]=%.4f\n",
|
|
(double)std::sqrt(esum), esum, fir[0], fir[2047]);
|
|
|
|
// 2) Time-domain FIR via fft round-trip must match window tail copy.
|
|
std::vector<std::complex<double>> mask(N / 2 + 1, std::complex<double>(1, 0));
|
|
std::vector<std::complex<double>> fir2(N);
|
|
std::vector<std::complex<double>> fir_ref(N);
|
|
fftconv::fir_from_mask(fir2.data(), mask.data(), WIN_WINDOW, N, &plan);
|
|
// inverse FFT then normalize by N (radix-2 inv has 1/N?) — check factor.
|
|
double peak = 0.0;
|
|
for (size_t i = 0; i < N; i++) {
|
|
double r = std::fabs(fir2[i].real());
|
|
if (r > peak) peak = r;
|
|
}
|
|
std::printf("fir_from_mask peak=%.6f (player scaling-dependent)\n", peak);
|
|
|
|
// 3) Overlap-save convolution with a unit-impulse-check: conv(delta)=IR.
|
|
{
|
|
std::vector<float> in(N, 0.0f), out(N, 0.0f);
|
|
in[0] = 1.0f;
|
|
fftconv::conv_overlap_save(fir.data(), N, N / 2,
|
|
in.data(), out.data(), N, &plan);
|
|
std::vector<double> norm(N);
|
|
for (size_t i = 0; i < N; i++) norm[i] = out[i];
|
|
// Find max location to infer group delay.
|
|
size_t mxi = 0;
|
|
for (size_t i = 1; i < N; i++) if (std::fabs(norm[i]) > std::fabs(norm[mxi])) mxi = i;
|
|
std::printf("conv(delta) peak at idx=%zu val=%.4f (was %.4f) — group delay check\n",
|
|
mxi, norm[mxi], fir[mxi]);
|
|
}
|
|
|
|
std::printf("fftconv integration check done\n");
|
|
return 0;
|
|
} |