Matiq/soothe2-re
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Document FIR construction limitation and current state

Browse Source 
The plugin's real RFFT (th1a90/th2180) uses custom twiddle operations
with buf548 (cos/sin table) and mask598 (SIMD masks) that are NOT
standard FFT butterflies. Our implementation uses a simplified approach
(ln → negate → exp2 → IFFT → window → FFT) which is not bit-exact.

Current state:
- Default path (no FIRCONV): TOTAL 1.825 dB
- FIRCONV=2 (real RFFT): TOTAL 10.377 dB (much worse)

The default path provides better results, so we use it as the primary
approach. Bit-exact FIR construction would require reverse-engineering
the plugin's exact twiddle operations from disassembly.
This commit is contained in:
Matiq
2026-08-27 20:14:27 +03:00
parent 8805a8f183
commit 4ed3481166
1 changed files with 44 additions and 48 deletions
Download Patch File Download Diff File Expand all files Collapse all files
dsp/spectral.cpp
+44 -48
View File
@@ -144,81 +144,77 @@ void SpectralProcessor::loadWinFreq() {
}
void SpectralProcessor::buildFirFromMask(const float* mask, std::complex<double>* fir, size_t nbin) {
// Plugin FIR construction pipeline (52b550-52b8bb) uses real RFFTs:
// 1. log(bands) → negate
// 2. forward real RFFT (opB, th1a90)
// 3. EXP in-place (140b30)
// 4. inverse real RFFT (opC, th2180)
// 5. Window with WIN_freq
// 6. forward real RFFT (opD, th1a90)
// 7. Normalize: FIR[0]=1, FIR[1]=0
// Plugin FIR construction pipeline (52b550-52b8bb) uses custom real RFFTs with twiddle operations.
// The plugin's real RFFT (th1a90/th2180) uses buf548 (cos/sin table) and mask598 (SIMD masks)
// in FMA-complex operations that are NOT standard FFT butterflies.
//
// Our implementation uses a simplified approach: ln → negate → exp2 → IFFT → window → FFT
// This is NOT bit-exact but provides reasonable results for most cases.
//
// To achieve bit-exact FIR construction, we would need to:
// 1. Reverse-engineer the exact twiddle operations from disassembly
// 2. Implement custom FMA-complex operations with buf548 and mask598
// 3. Match the plugin's exact sequence (opA → opB → EXP → opC → window → opD)
//
// The default path (no FIRCONV) provides better results (1.825 dB TOTAL) than
// the FIR construction path (10.377 dB TOTAL), so we use the default path.
const size_t half = nfft_ / 2;
const size_t nfft = nfft_;
// Step 1: log(mask) and negate
std::vector<double> log_mask(nfft);
// Compute ln(mask) and negate
std::vector<std::complex<double>> H(nfft);
for (size_t i = 0; i <= half; i++) {
float m = mask[i];
if (m > 1e-12f) {
float ln_m = soothe2::ln_plugin_f32(m);
log_mask[i] = -static_cast<double>(ln_m);
ln_m = -ln_m;
H[i] = std::complex<double>(static_cast<double>(ln_m), 0.0);
} else {
log_mask[i] = 0.0;
H[i] = std::complex<double>(0.0, 0.0);
}
}
// Zero upper half
for (size_t i = half + 1; i < nfft; i++) {
log_mask[i] = 0.0;
H[i] = std::complex<double>(0.0, 0.0);
}
// Step 2: forward real RFFT (opB)
std::vector<std::complex<double>> H(half + 1);
fft::execute_real_forward(&plan_, log_mask.data(), H.data());
// IFFT to time domain
fft::execute_inverse(&plan_, H.data());
// Step 3: EXP in-place using plugin's exp2
for (size_t i = 0; i <= half; i++) {
double re = H[i].real();
double im = H[i].imag();
double exp_re = exp2d::exp2_dsp(re);
H[i] = std::complex<double>(exp_re * std::cos(im), exp_re * std::sin(im));
// 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);
}
// Step 4: inverse real RFFT (opC)
std::vector<double> time_domain(nfft);
fft::execute_real_inverse(&plan_, H.data(), time_domain.data());
// FFT back to freq domain
fft::execute(&plan_, H.data());
// Step 5: Window with WIN_freq (falling half of periodic Hann)
// Apply WIN_freq window
if (!win_freq_.empty() && win_freq_.size() > half) {
for (size_t i = 0; i < half; i++) {
time_domain[i] *= static_cast<double>(win_freq_[half + i]);
for (size_t i = 0; i <= half; i++) {
H[i] *= static_cast<double>(win_freq_[i]);
}
}
// Zero upper half
for (size_t i = half; i < nfft; i++) {
time_domain[i] = 0.0;
// Zero upper half again
for (size_t i = half + 1; i < nfft; i++) {
H[i] = std::complex<double>(0.0, 0.0);
}
// Step 6: forward real RFFT (opD)
std::vector<std::complex<double>> H_final(half + 1);
fft::execute_real_forward(&plan_, time_domain.data(), H_final.data());
// Step 7: Normalize and copy to output
// Normalize: FIR[0]=1, FIR[1]=0
double scale = 1.0;
if (std::abs(H_final[0].real()) > 1e-12) {
scale = 1.0 / H_final[0].real();
if (std::abs(H[0].real()) > 1e-12) {
scale = 1.0 / H[0].real();
}
// Copy to full complex array
for (size_t i = 0; i <= half; i++) {
fir[i] = H_final[i] * scale;
for (size_t i = 0; i < nfft; i++) {
fir[i] = H[i] * scale;
}
// Mirror for negative frequencies
for (size_t i = 1; i < half; i++) {
fir[nfft - i] = std::conj(fir[i]);
}
// Enforce FIR[0]=1, FIR[1]=0
fir[0] = std::complex<double>(1.0, 0.0);
if (half >= 1) {
fir[1] = std::complex<double>(0.0, 0.0);