GCC Code Coverage Report

Directory: ./
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 94.0% 267 / 0 / 284
Functions: 90.9% 10 / 0 / 11
Branches: 58.5% 131 / 0 / 224
src/cpu/stft.cpp
Line Branch Exec Source
1 // ─── CPU STFT / iSTFT (brosoundml CHUNK 2) ─────────────────────────────────
2 //
3 // Short-time Fourier transform and its inverse, plus their adjoints. CPU
4 // backend, FP32-only. No external libraries — the per-frame DFT reuses the
5 // hand-rolled mixed-radix + Bluestein engine from detail/cpu/fft_core.h
6 // (so n_fft = 400 and prime sizes both work, exactly as in fft.cpp).
7 //
8 // Ops implemented here:
9 // stft / stft_backward real signal <-> complex spectrogram
10 // istft / istft_backward complex spectrogram <-> real signal (COLA OLA)
11 //
12 // ── Layout (see the doc comments in ops.h for the full contract) ────────────
13 // signal: REAL (N, signal_len) — N batched signals, one / row.
14 // window: REAL (1, win_length) — caller-supplied.
15 // spec: interleaved-complex (N*frames, 2*bins), bins = n_fft/2+1. Each
16 // frame is a row; the N signals' frame blocks are stacked in order.
17 //
18 // ── Frame model ─────────────────────────────────────────────────────────────
19 // Frame f of signal b takes n_fft samples starting at padded position
20 // f*hop_length, multiplies the central win_length of them by `window`, and
21 // rfft's the n_fft buffer. The window sits centred in the n_fft buffer
22 // (pad = (n_fft-win_length)/2 zeros each side). When center == true the
23 // signal is reflect-padded by n_fft/2 each side first; otherwise the raw
24 // signal is used. `padded_index` below maps a padded position back to a raw
25 // signal index (reflecting at the borders when center == true) so the forward
26 // op and its adjoint share one indexing rule and stay exact transposes.
27 //
28 // ── Normalisation ───────────────────────────────────────────────────────────
29 // rfft uses the "backward" convention (forward unscaled). normalized == true
30 // multiplies the forward spectrum by 1/sqrt(n_fft) (istft divides by it).
31 //
32 // ── Gradient design ─────────────────────────────────────────────────────────
33 // stft and istft are linear but NOT mutual adjoints (window + COLA). Each
34 // backward op is the exact transpose of its own forward linear map — see the
35 // ops.h header note. They are the minimal correct set for the
36 // multi-resolution STFT loss.
37
38 #include <brotensor/detail/cpu/fft_core.h>
39 #include <brotensor/detail/cpu/thread_pool.h>
40 #include <brotensor/tensor.h>
41
42 #include <algorithm>
43 #include <cmath>
44 #include <cstddef>
45 #include <stdexcept>
46 #include <string>
47 #include <vector>
48
49 namespace brotensor::detail::cpu {
50
51 using fftcore::Cd;
52 using fftcore::dft_1d;
53
54 namespace {
55
56 ✗ [[noreturn]] void fail(const char* op, const std::string& reason) {
57 ✗ throw std::runtime_error(std::string("brotensor: ") + op + ": " + reason);
58 ✗ }
59
60 12838 void require_fp32_host(const char* op, const ::brotensor::Tensor& t,
61 const char* name) {
62
1/2
✓ Branch 0 taken 12838 times.
✗ Branch 1 not taken.
 12838 if (t.device != ::brotensor::Device::CPU) {
63 ✗ fail(op, std::string(name) + " must be a CPU tensor");
64 }
65
1/2
✓ Branch 0 taken 12838 times.
✗ Branch 1 not taken.
 12838 if (t.dtype != ::brotensor::Dtype::FP32) {
66 ✗ fail(op, std::string(name) + " must be FP32 (CPU is FP32-only)");
67 }
68 12838 }
69
70 // Common parameter validation + derived sizes for all four ops.
71 6419 struct StftGeom {
72 6419 int bins = 0; // n_fft/2 + 1
73 6419 int frames = 0; // frames per signal
74 6419 int padded_len = 0; // signal length the frame loop indexes into
75 6419 int pad_lo = 0; // (n_fft - win_length) / 2 — window offset in buffer
76 };
77
78 6419 StftGeom check_geom(const char* op, int N, int signal_len, int n_fft,
79 int hop_length, int win_length, bool center) {
80
1/4
✓ Branch 0 taken 6419 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 6419 if (N < 0) fail(op, "N must be >= 0");
81
1/4
✓ Branch 0 taken 6419 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 6419 if (n_fft < 1) fail(op, "n_fft must be >= 1");
82
1/4
✓ Branch 0 taken 6419 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 6419 if (hop_length < 1) fail(op, "hop_length must be >= 1");
83
1/2
✓ Branch 0 taken 6419 times.
✗ Branch 1 not taken.
 6419 if (win_length < 1 || win_length > n_fft) {
84 ✗ fail(op, "win_length must satisfy 1 <= win_length <= n_fft");
85 }
86
1/4
✓ Branch 0 taken 6419 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 6419 if (signal_len < 1) fail(op, "signal_len must be >= 1");
87
88 6419 StftGeom g;
89 6419 g.bins = n_fft / 2 + 1;
90 6419 g.pad_lo = (n_fft - win_length) / 2;
91
92
2/2
✓ Branch 0 taken 3661 times.
✓ Branch 1 taken 2758 times.
 6419 if (center) {
93 // Reflect padding by n_fft/2 each side. numpy/torch 'reflect' mode
94 // needs at least 2 samples (the reflected index must stay in range);
95 // require enough signal to fill the n_fft/2 pad.
96
1/2
✓ Branch 0 taken 3661 times.
✗ Branch 1 not taken.
 3661 if (signal_len < n_fft / 2 + 1) {
97 ✗ fail(op, "center=true needs signal_len >= n_fft/2 + 1");
98 }
99 3661 g.padded_len = signal_len + n_fft;
100 3661 g.frames = 1 + signal_len / hop_length;
101 3661 } else {
102
1/2
✓ Branch 0 taken 2758 times.
✗ Branch 1 not taken.
 2758 if (signal_len < n_fft) {
103 ✗ fail(op, "center=false needs signal_len >= n_fft");
104 }
105 2758 g.padded_len = signal_len;
106 2758 g.frames = 1 + (signal_len - n_fft) / hop_length;
107 }
108 6419 return g;
109 ✗ }
110
111 // Map a padded position p in [0, padded_len) to a raw signal index in
112 // [0, signal_len). center == false is the identity; center == true reflects
113 // at the borders (numpy 'reflect': edge sample not repeated).
114 //
115 // Reflection over [0, L-1] with period 2*(L-1): fold q into that range.
116 184512 inline int reflect_index(int q, int L) {
117
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 184512 times.
 184512 if (L == 1) return 0;
118 184512 const int period = 2 * (L - 1);
119 184512 int m = q % period;
120
2/2
✓ Branch 0 taken 173636 times.
✓ Branch 1 taken 10876 times.
 184512 if (m < 0) m += period;
121
2/2
✓ Branch 0 taken 166890 times.
✓ Branch 1 taken 17622 times.
 184512 return (m < L) ? m : period - m;
122 184512 }
123
124 318027 inline int padded_index(int p, int signal_len, int n_fft, bool center) {
125
2/2
✓ Branch 0 taken 134155 times.
✓ Branch 1 taken 183872 times.
 318027 if (!center) return p;
126 183872 return reflect_index(p - n_fft / 2, signal_len);
127 318027 }
128
129 } // namespace
130
131 // ════════════════════════════════════════════════════════════════════════════
132 // stft — real signal -> complex spectrogram
133 // ════════════════════════════════════════════════════════════════════════════
134 1286 void stft(const ::brotensor::Tensor& signal, const ::brotensor::Tensor& window,
135 int N, int n_fft, int hop_length, int win_length,
136 bool center, bool normalized, ::brotensor::Tensor& spec) {
137 1286 require_fp32_host("stft", signal, "signal");
138 1286 require_fp32_host("stft", window, "window");
139
1/2
✓ Branch 0 taken 1286 times.
✗ Branch 1 not taken.
 1286 if (signal.rows != N) {
140 ✗ fail("stft", "signal.rows must equal N");
141 }
142 1286 const int signal_len = signal.cols;
143
1/2
✓ Branch 0 taken 1286 times.
✗ Branch 1 not taken.
 1286 if (window.rows != 1 || window.cols != win_length) {
144 ✗ fail("stft", "window must be a (1, win_length) tensor");
145 }
146 2572 const StftGeom g = check_geom("stft", N, signal_len, n_fft, hop_length,
147 1286 win_length, center);
148 1286 const int out_rows = N * g.frames;
149 1286 const int out_cols = 2 * g.bins;
150
3/4
✓ Branch 0 taken 1256 times.
✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1256 times.
 1286 if (spec.rows != out_rows || spec.cols != out_cols) {
151 30 spec.resize(out_rows, out_cols);
152 30 }
153
1/2
✓ Branch 0 taken 1286 times.
✗ Branch 1 not taken.
 1286 if (out_rows == 0) return;
154
155 1286 const float* sig = signal.host_f32();
156 1286 const float* win = window.host_f32();
157 1286 float* sp = spec.host_f32_mut();
158
2/2
✓ Branch 0 taken 637 times.
✓ Branch 1 taken 649 times.
 1286 const double norm = normalized
159 637 ? 1.0 / std::sqrt(static_cast<double>(n_fft))
160 : 1.0;
161
162 // One frame per work item. Frames are independent — each reads its own slice
163 // of the signal and writes its own row of `spec` — so this is the pool's
164 // ordinary case, and it is worth taking: a mel front end is thousands of
165 // frames of a double-precision DFT and it was the *host* half of an ASR
166 // encoder's cost. Measured on a FastConformer over 18 s of 16 kHz audio
167 // (1 801 frames, n_fft 512): 212 ms on one core.
168 //
169 // The scratch buffers are per item rather than hoisted, because two workers
170 // sharing them is a data race. fft_core's twiddle table is thread_local and
171 // keyed on (N, sign), so each worker builds it once and every frame after
172 // the first reuses it — which is why the buffers being fresh costs nothing.
173
2/4
✓ Branch 0 taken 1286 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1286 times.
✗ Branch 3 not taken.
 2572 parallel_for(static_cast<std::size_t>(N) * static_cast<std::size_t>(g.frames),
174 17623 [&](std::size_t item) {
175 16337 const int b = static_cast<int>(item / static_cast<std::size_t>(g.frames));
176 16337 const int f = static_cast<int>(item % static_cast<std::size_t>(g.frames));
177 16337 const float* srow = sig + static_cast<std::size_t>(b) * signal_len;
178
179 16337 std::vector<Cd> buf(static_cast<std::size_t>(n_fft)), out;
180 16337 const int base = f * hop_length; // padded-position start
181
2/2
✓ Branch 0 taken 16324 times.
✓ Branch 1 taken 315335 times.
 331659 for (int j = 0; j < win_length; ++j) {
182 315335 const int i = g.pad_lo + j;
183 315335 const int p = base + i;
184
2/2
✓ Branch 0 taken 315322 times.
✓ Branch 1 taken 13 times.
 315335 const int s = padded_index(p, signal_len, n_fft, center);
185 315322 buf[static_cast<std::size_t>(i)] =
186 315322 {static_cast<double>(srow[s]) * static_cast<double>(win[j]), 0.0};
187 315322 }
188
2/2
✓ Branch 0 taken 16311 times.
✓ Branch 1 taken 13 times.
 16324 dft_1d(buf, out, -1); // unscaled forward DFT
189 16311 float* dst = sp + item * static_cast<std::size_t>(out_cols);
190
2/2
✓ Branch 0 taken 177862 times.
✓ Branch 1 taken 16311 times.
 194173 for (int k = 0; k < g.bins; ++k) {
191 177862 dst[2 * k] = static_cast<float>(
192 177862 out[static_cast<std::size_t>(k)].re * norm);
193 177862 dst[2 * k + 1] = static_cast<float>(
194 177862 out[static_cast<std::size_t>(k)].im * norm);
195 177862 }
196 16337 });
197 1286 }
198
199 // ════════════════════════════════════════════════════════════════════════════
200 // stft_backward — adjoint of stft
201 // ════════════════════════════════════════════════════════════════════════════
202 //
203 // stft is the linear map spec = R * W * P * signal where P scatters the
204 // signal into frame buffers (with reflect padding folded in), W multiplies by
205 // the window, and R is the truncated forward DFT. Its adjoint applied to
206 // dSpec is P^T * W^T * R^T * dSpec :
207 // * R^T per frame is exactly rfft_backward's adjoint (the +1-sign unscaled
208 // DFT of the zero-padded n_fft spectrum, real part);
209 // * W^T is the same window multiply (diagonal — self-transpose);
210 // * P^T accumulates each frame sample back into the signal (the same index
211 // map, summed), so overlapping frames add — NO COLA division here (that
212 // belongs to istft, a different map).
213 // dSignal is *overwritten* (zeroed then accumulated).
214 11 void stft_backward(const ::brotensor::Tensor& dSpec,
215 const ::brotensor::Tensor& window,
216 int N, int signal_len, int n_fft, int hop_length,
217 int win_length, bool center, bool normalized,
218 ::brotensor::Tensor& dSignal) {
219 11 require_fp32_host("stft_backward", dSpec, "dSpec");
220 11 require_fp32_host("stft_backward", window, "window");
221
1/2
✓ Branch 0 taken 11 times.
✗ Branch 1 not taken.
 11 if (window.rows != 1 || window.cols != win_length) {
222 ✗ fail("stft_backward", "window must be a (1, win_length) tensor");
223 }
224 22 const StftGeom g = check_geom("stft_backward", N, signal_len, n_fft,
225 11 hop_length, win_length, center);
226 11 const int exp_rows = N * g.frames;
227 11 const int exp_cols = 2 * g.bins;
228
1/2
✓ Branch 0 taken 11 times.
✗ Branch 1 not taken.
 11 if (dSpec.rows != exp_rows || dSpec.cols != exp_cols) {
229 ✗ fail("stft_backward", "dSpec shape must match the stft output shape");
230 }
231
3/4
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 8 times.
 11 if (dSignal.rows != N || dSignal.cols != signal_len) {
232 3 dSignal.resize(N, signal_len);
233 3 }
234
1/2
✓ Branch 0 taken 11 times.
✗ Branch 1 not taken.
 11 if (dSignal.size() != 0) {
235 11 float* z = dSignal.host_f32_mut();
236
2/2
✓ Branch 0 taken 1120 times.
✓ Branch 1 taken 11 times.
 1131 for (int i = 0; i < dSignal.size(); ++i) z[i] = 0.0f;
237 11 }
238
1/2
✓ Branch 0 taken 11 times.
✗ Branch 1 not taken.
 11 if (exp_rows == 0) return;
239
240 11 const float* gp = dSpec.host_f32();
241 11 const float* win = window.host_f32();
242 11 float* dsig = dSignal.host_f32_mut();
243
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 6 times.
 11 const double norm = normalized
244 5 ? 1.0 / std::sqrt(static_cast<double>(n_fft))
245 : 1.0;
246
247 // Per frame: spec[k] = norm * (truncated DFT)[k]. The adjoint of the
248 // truncated forward DFT is: zero-pad dSpec to length n_fft, run an
249 // unscaled +1-sign DFT, take the real part (== rfft_backward's core).
250 11 std::vector<Cd> spec(static_cast<std::size_t>(n_fft)), tbuf;
251
2/2
✓ Branch 0 taken 23 times.
✓ Branch 1 taken 11 times.
 34 for (int b = 0; b < N; ++b) {
252 23 float* drow = dsig + static_cast<std::size_t>(b) * signal_len;
253
2/2
✓ Branch 0 taken 23 times.
✓ Branch 1 taken 251 times.
 274 for (int f = 0; f < g.frames; ++f) {
254 502 const float* grow = gp + static_cast<std::size_t>(
255 502 static_cast<std::size_t>(b) *
256 502 g.frames +
257 502 f) * exp_cols;
258
2/2
✓ Branch 0 taken 3168 times.
✓ Branch 1 taken 251 times.
 3419 for (int k = 0; k < n_fft; ++k) spec[static_cast<std::size_t>(k)] = Cd{};
259
2/2
✓ Branch 0 taken 1835 times.
✓ Branch 1 taken 251 times.
 2086 for (int k = 0; k < g.bins; ++k) {
260 1835 spec[static_cast<std::size_t>(k)] =
261 3670 {static_cast<double>(grow[2 * k]) * norm,
262 1835 static_cast<double>(grow[2 * k + 1]) * norm};
263 1835 }
264
1/2
✓ Branch 0 taken 251 times.
✗ Branch 1 not taken.
 251 dft_1d(spec, tbuf, +1); // adjoint of truncated forward DFT
265 // W^T (window) then P^T (scatter-add into the signal).
266 251 const int base = f * hop_length;
267
2/2
✓ Branch 0 taken 251 times.
✓ Branch 1 taken 2916 times.
 3167 for (int j = 0; j < win_length; ++j) {
268 2916 const int i = g.pad_lo + j;
269 2916 const int p = base + i;
270
1/2
✓ Branch 0 taken 2916 times.
✗ Branch 1 not taken.
 2916 const int s = padded_index(p, signal_len, n_fft, center);
271 2916 drow[s] += static_cast<float>(
272 5832 tbuf[static_cast<std::size_t>(i)].re *
273 2916 static_cast<double>(win[j]));
274 2916 }
275 251 }
276 23 }
277 11 }
278
279 // ════════════════════════════════════════════════════════════════════════════
280 // istft — complex spectrogram -> real signal (windowed overlap-add + COLA)
281 // ════════════════════════════════════════════════════════════════════════════
282 //
283 // Per frame: irfft the n_fft spectrum, multiply by the window, scatter-add
284 // into the output. Then divide each output sample by the overlap-added
285 // squared window (the COLA envelope) so a COLA-satisfying window+hop makes
286 // istft(stft(x)) == x. Samples with a ~0 envelope (edges with no frame
287 // coverage) stay 0.
288 5111 void istft(const ::brotensor::Tensor& spec, const ::brotensor::Tensor& window,
289 int N, int signal_len, int n_fft, int hop_length, int win_length,
290 bool center, bool normalized, ::brotensor::Tensor& signal) {
291 5111 require_fp32_host("istft", spec, "spec");
292 5111 require_fp32_host("istft", window, "window");
293
1/2
✓ Branch 0 taken 5111 times.
✗ Branch 1 not taken.
 5111 if (window.rows != 1 || window.cols != win_length) {
294 ✗ fail("istft", "window must be a (1, win_length) tensor");
295 }
296 10222 const StftGeom g = check_geom("istft", N, signal_len, n_fft, hop_length,
297 5111 win_length, center);
298 5111 const int exp_rows = N * g.frames;
299 5111 const int exp_cols = 2 * g.bins;
300
1/2
✓ Branch 0 taken 5111 times.
✗ Branch 1 not taken.
 5111 if (spec.rows != exp_rows || spec.cols != exp_cols) {
301 ✗ fail("istft", "spec shape must match the stft output shape");
302 }
303
3/4
✓ Branch 0 taken 5096 times.
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 5096 times.
 5111 if (signal.rows != N || signal.cols != signal_len) {
304 15 signal.resize(N, signal_len);
305 15 }
306
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5111 times.
 5111 if (signal.size() != 0) {
307 5111 float* z = signal.host_f32_mut();
308
2/2
✓ Branch 0 taken 451096 times.
✓ Branch 1 taken 5111 times.
 456207 for (int i = 0; i < signal.size(); ++i) z[i] = 0.0f;
309 5111 }
310
1/2
✓ Branch 0 taken 5111 times.
✗ Branch 1 not taken.
 5111 if (exp_rows == 0) return;
311
312 5111 const float* sp = spec.host_f32();
313 5111 const float* win = window.host_f32();
314 5111 float* sig = signal.host_f32_mut();
315 // istft inverts stft's optional 1/sqrt(n_fft): multiply the spectrum by
316 // sqrt(n_fft) so the irfft 1/n_fft scaling lands at the right amplitude.
317
2/2
✓ Branch 0 taken 2555 times.
✓ Branch 1 taken 2556 times.
 5111 const double norm = normalized ? std::sqrt(static_cast<double>(n_fft))
318 : 1.0;
319 5111 const double invN = 1.0 / static_cast<double>(n_fft);
320
321 // COLA envelope: overlap-added squared window, in padded coordinates.
322 5111 std::vector<double> env(static_cast<std::size_t>(g.padded_len), 0.0);
323
2/2
✓ Branch 0 taken 60119 times.
✓ Branch 1 taken 5111 times.
 65230 for (int f = 0; f < g.frames; ++f) {
324 60119 const int base = f * hop_length;
325
2/2
✓ Branch 0 taken 811508 times.
✓ Branch 1 taken 60119 times.
 871627 for (int j = 0; j < win_length; ++j) {
326 811508 const int p = base + g.pad_lo + j;
327 811508 const double w = static_cast<double>(win[j]);
328 811508 env[static_cast<std::size_t>(p)] += w * w;
329 811508 }
330 60119 }
331
332 // Per signal: overlap-add the windowed irfft frames, then COLA-divide.
333
1/2
✓ Branch 0 taken 5111 times.
✗ Branch 1 not taken.
 5111 std::vector<Cd> full(static_cast<std::size_t>(n_fft)), out;
334
1/2
✓ Branch 0 taken 5111 times.
✗ Branch 1 not taken.
 5111 std::vector<double> acc(static_cast<std::size_t>(g.padded_len));
335
2/2
✓ Branch 0 taken 10223 times.
✓ Branch 1 taken 5111 times.
 15334 for (int b = 0; b < N; ++b) {
336
1/2
✓ Branch 0 taken 10223 times.
✗ Branch 1 not taken.
 10223 std::fill(acc.begin(), acc.end(), 0.0);
337
2/2
✓ Branch 0 taken 120259 times.
✓ Branch 1 taken 10223 times.
 130482 for (int f = 0; f < g.frames; ++f) {
338 240518 const float* srow = sp + static_cast<std::size_t>(
339 240518 static_cast<std::size_t>(b) *
340 240518 g.frames +
341 240518 f) * exp_cols;
342 // Rebuild the Hermitian-symmetric n_fft spectrum, irfft it.
343
2/2
✓ Branch 0 taken 932019 times.
✓ Branch 1 taken 120259 times.
 1052278 for (int k = 0; k < g.bins; ++k) {
344 932019 full[static_cast<std::size_t>(k)] =
345 1864038 {static_cast<double>(srow[2 * k]) * norm,
346 932019 static_cast<double>(srow[2 * k + 1]) * norm};
347 932019 }
348
2/2
✓ Branch 0 taken 691501 times.
✓ Branch 1 taken 120259 times.
 811760 for (int k = 1; k < n_fft - g.bins + 1; ++k) {
349 691501 const Cd c = full[static_cast<std::size_t>(k)];
350 691501 full[static_cast<std::size_t>(n_fft - k)] = {c.re, -c.im};
351 691501 }
352
1/2
✓ Branch 0 taken 120259 times.
✗ Branch 1 not taken.
 120259 dft_1d(full, out, +1); // inverse DFT, still needs *1/n_fft
353 120259 const int base = f * hop_length;
354
2/2
✓ Branch 0 taken 1623268 times.
✓ Branch 1 taken 120259 times.
 1743527 for (int j = 0; j < win_length; ++j) {
355 1623268 const int i = g.pad_lo + j;
356 1623268 const int p = base + i;
357 1623268 const double t = out[static_cast<std::size_t>(i)].re * invN;
358 1623268 acc[static_cast<std::size_t>(p)] +=
359 1623268 t * static_cast<double>(win[j]);
360 1623268 }
361 120259 }
362 // COLA-divide and strip centre padding back to the raw signal.
363 10223 float* drow = sig + static_cast<std::size_t>(b) * signal_len;
364
2/2
✓ Branch 0 taken 6009 times.
✓ Branch 1 taken 4214 times.
 10223 const int shift = center ? n_fft / 2 : 0;
365
2/2
✓ Branch 0 taken 451096 times.
✓ Branch 1 taken 10223 times.
 461319 for (int n = 0; n < signal_len; ++n) {
366 451096 const int p = n + shift;
367 451096 const double e = env[static_cast<std::size_t>(p)];
368
2/2
✓ Branch 0 taken 437444 times.
✓ Branch 1 taken 13652 times.
 451096 drow[n] = (e > 1e-10)
369 874888 ? static_cast<float>(acc[static_cast<std::size_t>(p)]
370 437444 / e)
371 : 0.0f;
372 451096 }
373 10223 }
374 5111 }
375
376 // ════════════════════════════════════════════════════════════════════════════
377 // istft_backward — adjoint of istft
378 // ════════════════════════════════════════════════════════════════════════════
379 //
380 // istft is the linear map signal = D * E^{-1} * P * W * I * spec where I is
381 // the per-frame inverse DFT (1/n_fft scaled), W the window multiply, P the
382 // overlap-add scatter, E^{-1} the per-sample COLA division, and D the
383 // centre-padding strip. The COLA envelope E depends only on the window/hop,
384 // not on the spectrum, so E^{-1} is a (data-independent) diagonal — and its
385 // transpose is itself. The adjoint applied to dSignal is therefore
386 // I^T * W^T * P^T * E^{-1} * D^T * dSignal :
387 // * D^T scatters dSignal back into padded coordinates;
388 // * E^{-1} divides by the same COLA envelope (diagonal, self-transpose);
389 // * P^T gathers each frame's window_length samples;
390 // * W^T is the window multiply again;
391 // * I^T is the adjoint of the inverse DFT — which is irfft_backward's core
392 // (forward-sign DFT of the gathered frame, 1/n_fft scaling, and the
393 // interior-bin doubling from the Hermitian fold).
394 // dSpec is *overwritten*.
395 11 void istft_backward(const ::brotensor::Tensor& dSignal,
396 const ::brotensor::Tensor& window,
397 int N, int signal_len, int n_fft, int hop_length,
398 int win_length, bool center, bool normalized,
399 ::brotensor::Tensor& dSpec) {
400 11 require_fp32_host("istft_backward", dSignal, "dSignal");
401 11 require_fp32_host("istft_backward", window, "window");
402
1/2
✓ Branch 0 taken 11 times.
✗ Branch 1 not taken.
 11 if (dSignal.rows != N || dSignal.cols != signal_len) {
403 ✗ fail("istft_backward", "dSignal must be a (N, signal_len) tensor");
404 }
405
1/2
✓ Branch 0 taken 11 times.
✗ Branch 1 not taken.
 11 if (window.rows != 1 || window.cols != win_length) {
406 ✗ fail("istft_backward", "window must be a (1, win_length) tensor");
407 }
408 22 const StftGeom g = check_geom("istft_backward", N, signal_len, n_fft,
409 11 hop_length, win_length, center);
410 11 const int out_rows = N * g.frames;
411 11 const int out_cols = 2 * g.bins;
412
3/4
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 8 times.
 11 if (dSpec.rows != out_rows || dSpec.cols != out_cols) {
413 3 dSpec.resize(out_rows, out_cols);
414 3 }
415
1/2
✓ Branch 0 taken 11 times.
✗ Branch 1 not taken.
 11 if (out_rows == 0) return;
416
417 11 const float* dsig = dSignal.host_f32();
418 11 const float* win = window.host_f32();
419 11 float* gp = dSpec.host_f32_mut();
420
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 6 times.
 11 const double norm = normalized ? std::sqrt(static_cast<double>(n_fft))
421 : 1.0;
422 11 const double invN = 1.0 / static_cast<double>(n_fft);
423
424 // Same COLA envelope as istft (window/hop only).
425 11 std::vector<double> env(static_cast<std::size_t>(g.padded_len), 0.0);
426
2/2
✓ Branch 0 taken 159 times.
✓ Branch 1 taken 11 times.
 170 for (int f = 0; f < g.frames; ++f) {
427 159 const int base = f * hop_length;
428
2/2
✓ Branch 0 taken 1844 times.
✓ Branch 1 taken 159 times.
 2003 for (int j = 0; j < win_length; ++j) {
429 1844 const int p = base + g.pad_lo + j;
430 1844 const double w = static_cast<double>(win[j]);
431 1844 env[static_cast<std::size_t>(p)] += w * w;
432 1844 }
433 159 }
434
435 11 const bool even = (n_fft % 2 == 0);
436
1/2
✓ Branch 0 taken 11 times.
✗ Branch 1 not taken.
 11 std::vector<double> gacc(static_cast<std::size_t>(g.padded_len));
437
1/2
✓ Branch 0 taken 11 times.
✗ Branch 1 not taken.
 11 std::vector<Cd> frame(static_cast<std::size_t>(n_fft)), spec;
438
2/2
✓ Branch 0 taken 23 times.
✓ Branch 1 taken 11 times.
 34 for (int b = 0; b < N; ++b) {
439 // D^T then E^{-1}: scatter dSignal into padded coords, COLA-divide.
440
1/2
✓ Branch 0 taken 23 times.
✗ Branch 1 not taken.
 23 std::fill(gacc.begin(), gacc.end(), 0.0);
441 23 const float* drow = dsig + static_cast<std::size_t>(b) * signal_len;
442
2/2
✓ Branch 0 taken 13 times.
✓ Branch 1 taken 10 times.
 23 const int shift = center ? n_fft / 2 : 0;
443
2/2
✓ Branch 0 taken 1120 times.
✓ Branch 1 taken 23 times.
 1143 for (int n = 0; n < signal_len; ++n) {
444 1120 const int p = n + shift;
445 1120 const double e = env[static_cast<std::size_t>(p)];
446 1120 gacc[static_cast<std::size_t>(p)] =
447
2/2
✓ Branch 0 taken 1092 times.
✓ Branch 1 taken 28 times.
 1120 (e > 1e-10) ? static_cast<double>(drow[n]) / e : 0.0;
448 1120 }
449
2/2
✓ Branch 0 taken 23 times.
✓ Branch 1 taken 339 times.
 362 for (int f = 0; f < g.frames; ++f) {
450 // P^T (gather) then W^T (window): the frame's n_fft time buffer.
451
2/2
✓ Branch 0 taken 4192 times.
✓ Branch 1 taken 339 times.
 4531 for (int k = 0; k < n_fft; ++k) frame[static_cast<std::size_t>(k)] = Cd{};
452 339 const int base = f * hop_length;
453
2/2
✓ Branch 0 taken 3940 times.
✓ Branch 1 taken 339 times.
 4279 for (int j = 0; j < win_length; ++j) {
454 3940 const int i = g.pad_lo + j;
455 3940 const int p = base + i;
456 3940 frame[static_cast<std::size_t>(i)] =
457 7880 {gacc[static_cast<std::size_t>(p)] *
458 3940 static_cast<double>(win[j]),
459 0.0};
460 3940 }
461 // I^T: adjoint of the 1/n_fft inverse DFT — forward-sign DFT,
462 // 1/n_fft scaling, interior-bin doubling for the Hermitian fold.
463
1/2
✓ Branch 0 taken 339 times.
✗ Branch 1 not taken.
 339 dft_1d(frame, spec, -1);
464 678 float* grow = gp + static_cast<std::size_t>(
465 678 static_cast<std::size_t>(b) * g.frames +
466 678 f) * out_cols;
467
2/2
✓ Branch 0 taken 2435 times.
✓ Branch 1 taken 339 times.
 2774 for (int k = 0; k < g.bins; ++k) {
468 2435 double s = 2.0;
469
2/2
✓ Branch 0 taken 2096 times.
✓ Branch 1 taken 339 times.
 2435 if (k == 0) s = 1.0;
470
3/4
✓ Branch 0 taken 2435 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2096 times.
✓ Branch 3 taken 339 times.
 2435 if (even && k == n_fft / 2) s = 1.0;
471 2435 const double scale = s * invN * norm;
472 2435 grow[2 * k] = static_cast<float>(
473 2435 scale * spec[static_cast<std::size_t>(k)].re);
474 2435 grow[2 * k + 1] = static_cast<float>(
475 2435 scale * spec[static_cast<std::size_t>(k)].im);
476 2435 }
477 339 }
478 23 }
479 11 }
480
481 } // namespace brotensor::detail::cpu
482