Matiq
f689023089
- log2_ln.hpp/cpp: Plugin's exact ln(float) polynomial from 535a70 (0x1802a24c0). IEEE 754 bit extraction + Horner evaluation. Coefficients extracted from binary at 0x181f81f80..0x181f821c0. Max error ~3e-6 for typical inputs. - spectral.cpp: Updated buildFirFromMask to use plugin's ln→negate→exp2 pipeline instead of naive 1/mask reciprocal. - exp2_tables.hpp/cpp: Already contains plugin's exp2 tables (0x26b820). Remaining: twiddle stages (ops B/C/D with cos/sin tables from buf548) are the missing piece for bit-exact FIR construction. These are FFT butterflies already implemented in fft.hpp but need integration into the FIR pipeline.
301 lines
11 KiB
C++
301 lines
11 KiB
C++
#include "spectral.hpp"
|
||
#include "fftconv.hpp"
|
||
#include "log2_ln.hpp"
|
||
#include "exp2_tables.hpp"
|
||
#include "exp2.hpp"
|
||
#include <cmath>
|
||
#include <cstring>
|
||
#include <vector>
|
||
#include <cstdlib>
|
||
#include <cstdio>
|
||
|
||
SpectralProcessor::SpectralProcessor(size_t nfft, size_t hop, float sample_rate)
|
||
: nfft_(nfft), hop_(hop), frame_count_(0), output_pos_(0),
|
||
detector_(nfft, sample_rate) {
|
||
window_ = new double[nfft_];
|
||
computeWindow();
|
||
fft::init_plan(&plan_, static_cast<uint32_t>(std::log2(nfft_)));
|
||
buf_ = new std::complex<double>[nfft_];
|
||
tmp_buf_ = new std::complex<double>[nfft_];
|
||
fir_buf_ = new std::complex<double>[nfft_];
|
||
fir_freq_ = new std::complex<double>[nfft_];
|
||
overlap_.resize(nfft_, 0.0f);
|
||
mask_.resize(nfft_, 1.0f);
|
||
|
||
// Build FIR window: falling half of periodic Hann(4096).
|
||
// Plugin reads window[N/2..N-1] of periodic Hann (rising 0→1).
|
||
fir_window_.resize(nfft_);
|
||
for (size_t i = 0; i < nfft_; i++) {
|
||
fir_window_[i] = 0.5 * (1.0 - std::cos(2.0 * M_PI * i / nfft_));
|
||
}
|
||
}
|
||
|
||
SpectralProcessor::~SpectralProcessor() {
|
||
delete[] window_;
|
||
delete[] buf_;
|
||
delete[] tmp_buf_;
|
||
delete[] fir_buf_;
|
||
delete[] fir_freq_;
|
||
}
|
||
|
||
void SpectralProcessor::setDetectorParams(const std::vector<DetectorBand>& bands) {
|
||
detector_.setParams(bands);
|
||
loadWinFreq();
|
||
}
|
||
|
||
void SpectralProcessor::computeWindow() {
|
||
// RT_WIN: 0=symmetric hann (legacy), 1=periodic hann, 2=rectangular
|
||
static const int winmode = getenv("RT_WIN") ? atoi(getenv("RT_WIN")) : 0;
|
||
for (size_t i = 0; i < nfft_; i++) {
|
||
double v;
|
||
if (winmode == 1) v = 0.5 * (1.0 - std::cos(2.0 * M_PI * i / nfft_));
|
||
else if (winmode == 2) v = 1.0;
|
||
else v = 0.5 * (1.0 - std::cos(2.0 * M_PI * i / (nfft_ - 1)));
|
||
window_[i] = v;
|
||
}
|
||
}
|
||
|
||
void SpectralProcessor::stftFrame(const float* in, std::complex<double>* out) {
|
||
for (size_t i = 0; i < nfft_; i++) {
|
||
out[i] = std::complex<double>(static_cast<double>(in[i]) * window_[i], 0.0);
|
||
}
|
||
fft::execute(&plan_, out);
|
||
}
|
||
|
||
void SpectralProcessor::istftFrame(std::complex<double>* in, float* out, float* overlap) {
|
||
memcpy(tmp_buf_, in, nfft_ * sizeof(std::complex<double>));
|
||
fft::execute_inverse(&plan_, tmp_buf_);
|
||
static bool wola_computed = false;
|
||
static float wola_norm = 1.0f;
|
||
// RT_SYN: 0=synthesis window = analysis window (WOLA), 1=none
|
||
static const int synmode = getenv("RT_SYN") ? atoi(getenv("RT_SYN")) : 0;
|
||
if (!wola_computed) {
|
||
double wola_sum = 0.0;
|
||
for (size_t i = 0; i < nfft_; i++) {
|
||
double w = (synmode == 1) ? 1.0 : window_[i];
|
||
wola_sum += window_[i] * w;
|
||
}
|
||
wola_norm = static_cast<float>(wola_sum / hop_);
|
||
wola_computed = true;
|
||
}
|
||
for (size_t i = 0; i < nfft_; i++) {
|
||
double w = (synmode == 1) ? 1.0f : window_[i];
|
||
overlap[i] += static_cast<float>(tmp_buf_[i].real() * w);
|
||
}
|
||
for (size_t i = 0; i < hop_; i++) {
|
||
out[i] = overlap[i] / wola_norm;
|
||
}
|
||
for (size_t i = 0; i < nfft_ - hop_; i++) {
|
||
overlap[i] = overlap[i + hop_];
|
||
}
|
||
for (size_t i = nfft_ - hop_; i < nfft_; i++) {
|
||
overlap[i] = 0.0f;
|
||
}
|
||
}
|
||
|
||
void SpectralProcessor::loadWinFreq() {
|
||
if (win_freq_loaded_) return;
|
||
win_freq_loaded_ = true;
|
||
// Try to load WIN_freq from live capture (handoff/rtwin_freq_44100.npy)
|
||
FILE* f = fopen("handoff/rtwin_freq_44100.npy", "rb");
|
||
if (!f) {
|
||
// Fallback: compute periodic Hann, second half (0.5→1.0 rising)
|
||
win_freq_.resize(nfft_ / 2 + 1);
|
||
for (size_t i = 0; i <= nfft_ / 2; i++) {
|
||
win_freq_[i] = static_cast<float>(0.5 * (1.0 - std::cos(2.0 * M_PI * i / nfft_)));
|
||
}
|
||
return;
|
||
}
|
||
// Read numpy header
|
||
char header[128];
|
||
if (fread(header, 1, 6, f) != 6) { fclose(f); return; }
|
||
// Skip to data (numpy format: magic + header_len + desc)
|
||
fseek(f, 0, SEEK_END);
|
||
long fsize = ftell(f);
|
||
fseek(f, 0, SEEK_SET);
|
||
// Simple approach: skip header until '\n' appears, then read raw float32
|
||
fseek(f, 0, SEEK_SET);
|
||
int c;
|
||
while ((c = fgetc(f)) != '\n' && c != EOF) {}
|
||
// Read count (should be 8193 for 44100)
|
||
int32_t count = 0;
|
||
fread(&count, 4, 1, f);
|
||
// Actually numpy header is more complex; just read all remaining as float32
|
||
fseek(f, 0, SEEK_SET);
|
||
// Skip to data: find first 'N' (for 'astype') then skip past it
|
||
fseek(f, 6, SEEK_SET);
|
||
while ((c = fgetc(f)) != '\n' && c != EOF) {}
|
||
// Now at data start. Read until we have enough floats
|
||
std::vector<float> raw;
|
||
float val;
|
||
while (fread(&val, 4, 1, f) == 1) {
|
||
raw.push_back(val);
|
||
}
|
||
fclose(f);
|
||
if (raw.size() > 0) {
|
||
win_freq_ = raw;
|
||
} else {
|
||
// Fallback
|
||
win_freq_.resize(nfft_ / 2 + 1);
|
||
for (size_t i = 0; i <= nfft_ / 2; i++) {
|
||
win_freq_[i] = static_cast<float>(0.5 * (1.0 - std::cos(2.0 * M_PI * i / nfft_)));
|
||
}
|
||
}
|
||
}
|
||
|
||
void SpectralProcessor::buildFirFromMask(const float* mask, std::complex<double>* fir, size_t nbin) {
|
||
// Exact plugin FIR construction pipeline (52b550-52b8bb):
|
||
// 1. bands *= s888 (wet scale) - already applied to mask
|
||
// 2. th2270 - scalar transform - already in detector
|
||
// 3. 535a70: scratch = log(bands) - NATURAL LOG via plugin polynomial
|
||
// 4. Sign inversion: FIR[1..n/2] /= -1 (negate log = 1/bands after exp)
|
||
// 5. Zero upper half
|
||
// 6. opB: FMA twiddle (FFT butterfly with cos/sin)
|
||
// 7. BIGKERNEL 140b30: EXP in-place (exp2 via plugin tables)
|
||
// 8. opC: FMA twiddle
|
||
// 9. Window with WIN_freq
|
||
// 10. Zero upper half
|
||
// 11. opD: FMA twiddle
|
||
// 12. FIR[0]=1, FIR[1]=0
|
||
// 13. Scale by wet (already in mask)
|
||
// 14. df0: complex multiply FIR × audio spectrum
|
||
|
||
// Implementation matching plugin's log→negate→exp2 pipeline:
|
||
// mask → ln → negate → exp2 → IFFT → causal window → FFT → normalize
|
||
|
||
const size_t half = nfft_ / 2;
|
||
const size_t nfft = nfft_;
|
||
|
||
// Step 1-3: Compute ln(mask) using plugin's exact ln polynomial
|
||
// Then negate (sign inversion) → ln(1/mask)
|
||
// Then exp2 → 1/mask (reciprocal)
|
||
std::vector<float> log_mask(half + 1);
|
||
std::vector<float> recip_mask(half + 1);
|
||
|
||
for (size_t i = 0; i <= half; i++) {
|
||
float m = mask[i];
|
||
if (m > 1e-12f) {
|
||
// Plugin's ln polynomial
|
||
float ln_m = soothe2::ln_plugin_f32(m);
|
||
// Negate (sign inversion = divide by -1)
|
||
ln_m = -ln_m;
|
||
// Plugin's exp2 (exact from 0x26b820)
|
||
recip_mask[i] = static_cast<float>(exp2d::exp2_dsp(ln_m));
|
||
} else {
|
||
recip_mask[i] = 1.0f;
|
||
}
|
||
}
|
||
|
||
// Step 4-5: Zero upper half (Hermitian symmetry)
|
||
std::vector<std::complex<double>> H(nfft);
|
||
for (size_t i = 0; i <= half; i++) {
|
||
H[i] = std::complex<double>(static_cast<double>(recip_mask[i]), 0.0);
|
||
}
|
||
for (size_t i = half + 1; i < nfft; i++) {
|
||
H[i] = std::complex<double>(0.0, 0.0);
|
||
}
|
||
|
||
// Step 6-8: The twiddle ops (B/C/D) + EXP are effectively
|
||
// minimum-phase FIR design: IFFT → causal window → FFT
|
||
// Our fft::execute already matches plugin's FFT butterflies
|
||
|
||
// IFFT to time domain
|
||
fft::execute_inverse(&plan_, H.data());
|
||
|
||
// Causal window: keep first half, apply rising Hann (0.5→1.0)
|
||
for (size_t i = 0; i < half; i++) {
|
||
double win = 0.5 * (1.0 - std::cos(2.0 * M_PI * i / nfft));
|
||
H[i] *= win;
|
||
}
|
||
for (size_t i = half; i < nfft; i++) {
|
||
H[i] = std::complex<double>(0.0, 0.0);
|
||
}
|
||
|
||
// FFT back to freq domain
|
||
fft::execute(&plan_, H.data());
|
||
|
||
// Apply WIN_freq window (falling half of periodic Hann)
|
||
// But WIN_freq[n/2..n-1] is all 1.0, so this is no-op for lower half
|
||
if (!win_freq_.empty() && win_freq_.size() > half) {
|
||
for (size_t i = 0; i <= half; i++) {
|
||
H[i] *= static_cast<double>(win_freq_[i]);
|
||
}
|
||
}
|
||
|
||
// Zero upper half again
|
||
for (size_t i = half + 1; i < nfft; i++) {
|
||
H[i] = std::complex<double>(0.0, 0.0);
|
||
}
|
||
|
||
// Normalize: FIR[0]=1, FIR[1]=0
|
||
double scale = 1.0;
|
||
if (std::abs(H[0].real()) > 1e-12) {
|
||
scale = 1.0 / H[0].real();
|
||
}
|
||
for (size_t i = 0; i < nfft; i++) {
|
||
fir[i] = H[i] * scale;
|
||
}
|
||
fir[0] = std::complex<double>(1.0, 0.0);
|
||
if (half >= 1) {
|
||
fir[1] = std::complex<double>(0.0, 0.0);
|
||
}
|
||
}
|
||
|
||
void SpectralProcessor::processBlock(float* in, float* out, size_t num_samples, size_t num_channels) {
|
||
memset(out, 0, num_samples * sizeof(float));
|
||
if (num_samples == 0 || num_samples < nfft_) {
|
||
return;
|
||
}
|
||
|
||
static const int firconv = []() {
|
||
const char* e = getenv("RT_FIRCONV");
|
||
return e ? atoi(e) : 0;
|
||
}();
|
||
|
||
size_t nframes = (num_samples - nfft_) / hop_ + 1;
|
||
|
||
for (size_t f = 0; f < nframes; f++) {
|
||
size_t offset = f * hop_;
|
||
if (offset + nfft_ > num_samples) break;
|
||
stftFrame(in + offset, buf_);
|
||
|
||
detector_.processFrame(buf_, mask_.data());
|
||
|
||
if (firconv == 3) {
|
||
// RT_FIRCONV=3 (NOTES 24k): plugin application law decoded live:
|
||
// applied_gain = 1.019 * V^1.8345 per bin (rms 0.0025 dB over
|
||
// 8 drive levels). V = band curve (detector output); here M.
|
||
for (size_t i = 0; i < nfft_; i++) {
|
||
double m = std::max(static_cast<double>(mask_[i]), 1e-12);
|
||
double a = 1.019 * std::pow(m, 1.8345);
|
||
buf_[i] *= a;
|
||
}
|
||
} else if (firconv) {
|
||
// RT_FIRCONV=2: Full FIR construction pipeline (52b550-52b8bb).
|
||
// mask → reciprocal (1/mask) → window → normalize → complex multiply.
|
||
// This replicates the plugin's FFT-conv FIR design path.
|
||
buildFirFromMask(mask_.data(), fir_freq_, nfft_);
|
||
// Complex multiply FIR × audio spectrum
|
||
for (size_t i = 0; i < nfft_; i++) {
|
||
buf_[i] *= fir_freq_[i];
|
||
}
|
||
} else if (firconv == 1) {
|
||
// RT_FIRCONV=1: Simple frequency-domain mask multiply (legacy).
|
||
for (size_t i = 0; i < nfft_; i++) {
|
||
fir_freq_[i] = std::complex<double>(
|
||
static_cast<double>(mask_[i % (nfft_/2+1)]), 0.0);
|
||
}
|
||
for (size_t i = 0; i < nfft_; i++) {
|
||
buf_[i] *= fir_freq_[i];
|
||
}
|
||
} else {
|
||
// Default path: simple frequency-domain mask multiply.
|
||
for (size_t i = 0; i < nfft_; i++) {
|
||
buf_[i] *= mask_[i];
|
||
}
|
||
}
|
||
|
||
istftFrame(buf_, out + offset, overlap_.data());
|
||
}
|
||
}
|