tyler/rspektrum
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

refactor: extract shared ComputeSegment for STFT passes

Browse Source 
ComputeSTFTIncremental (overview, skipFactor-strided) and
ComputeNextHighResChunk (high-res fill) had ~50 lines of identical
per-segment code (windowing, FFT, bin fill for both the normal and
derivative spectra). Extract it into ComputeSegment() with a reusable
SegScratch buffer set (allocated once per pass, no per-segment malloc).
Each caller keeps only its own skip logic. Behavior-preserving — verified
the rendered spectrogram is unchanged.

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
This commit is contained in:
Tyler
2026-05-25 02:02:46 -07:00
parent 7240cf8ecf
commit 487e3ad85b
1 changed files with 90 additions and 125 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/stft.c
+90 -125
View File
@@ -135,74 +135,97 @@ void SaveToCache(void)
app.fftSize, app.stft.numSegments); app.fftSize, app.stft.numSegments);
} }
// ===== Background high-res computation ===== // ===== Per-segment STFT (shared by the overview and high-res passes) =====
int ComputeNextHighResChunk(AudioSignal* signal, StftResult* result,
int fftSize, int startSeg, int endSeg) // Scratch buffers reused across every segment in one pass, so we don't malloc
// per segment. Allocate once, hand to ComputeSegment, free when the pass ends.
typedef struct {
float* windowed;
float* derivWindowed;
float complex* fftIn;
float complex* fftOut;
} SegScratch;
static SegScratch AllocSegScratch(int fftSize)
{
SegScratch sc;
sc.windowed = (float*)malloc(fftSize * sizeof(float));
sc.derivWindowed = (float*)malloc(fftSize * sizeof(float));
sc.fftIn = (float complex*)malloc(fftSize * sizeof(float complex));
sc.fftOut = (float complex*)malloc(fftSize * sizeof(float complex));
return sc;
}
static void FreeSegScratch(SegScratch* sc)
{
free(sc->windowed);
free(sc->derivWindowed);
free(sc->fftIn);
free(sc->fftOut);
}
// Compute one STFT segment (normal V_f + derivative-window V_fd spectra) into
// result->segments[seg]. Caller ensures the segment isn't already computed.
static void ComputeSegment(AudioSignal* signal, StftResult* result, int fftSize, int seg, SegScratch* sc)
{ {
int hopSize = fftSize / HOP_RATIO; int hopSize = fftSize / HOP_RATIO;
int numBins = fftSize / 2 + 1; int numBins = fftSize / 2 + 1;
float* windowedSamples = (float*)malloc(fftSize * sizeof(float)); int offset = seg * hopSize;
float* derivWindowedSamples = (float*)malloc(fftSize * sizeof(float)); int samplesToCopy = fftSize;
float complex *complexInput = (float complex*)malloc(fftSize * sizeof(float complex)); if (offset + samplesToCopy > signal->numSamples) {
float complex* fftOutput = (float complex*)malloc(fftSize * sizeof(float complex)); samplesToCopy = signal->numSamples - offset;
memset(sc->windowed, 0, fftSize * sizeof(float));
for (int seg = startSeg; seg < endSeg && seg < result->numSegments; seg++) { memset(sc->derivWindowed, 0, fftSize * sizeof(float));
// Skip if already computed (overview or high-res) } else {
if (result->segments[seg].spectrum != NULL) continue; memcpy(sc->windowed, signal->samples + offset, fftSize * sizeof(float));
memcpy(sc->derivWindowed, signal->samples + offset, fftSize * sizeof(float));
int offset = seg * hopSize;
int samplesToCopy = fftSize;
if (offset + samplesToCopy > signal->numSamples) {
samplesToCopy = signal->numSamples - offset;
memset(windowedSamples, 0, fftSize * sizeof(float));
memset(derivWindowedSamples, 0, fftSize * sizeof(float));
} else {
memcpy(windowedSamples, signal->samples + offset, fftSize * sizeof(float));
memcpy(derivWindowedSamples, signal->samples + offset, fftSize * sizeof(float));
}
// Apply Hann window and derivative window
for (int i = 0; i < fftSize; i++) {
float t = (float)i / (fftSize - 1);
float hann = 0.5f * (1.0f - cosf(2.0f * M_PI * t));
float derivHann = M_PI * sinf(2.0f * M_PI * t);
windowedSamples[i] *= hann;
derivWindowedSamples[i] *= derivHann;
}
// Normal STFT
for (int i = 0; i < fftSize; i++) complexInput[i] = windowedSamples[i] + 0.0f * I;
FFT(complexInput, fftOutput, fftSize, false);
result->segments[seg].numBins = numBins;
result->segments[seg].sampleOffset = offset;
result->segments[seg].sampleCount = samplesToCopy;
result->segments[seg].spectrum = (FrequencyData*)malloc(numBins * sizeof(FrequencyData));
for (int bin = 0; bin < numBins; bin++) {
result->segments[seg].spectrum[bin].frequency = (float)bin * signal->sampleRate / fftSize;
result->segments[seg].spectrum[bin].amplitude = (bin == 0) ? cabsf(fftOutput[bin]) / fftSize : 2.0f * cabsf(fftOutput[bin]) / fftSize;
result->segments[seg].spectrum[bin].phase = cargf(fftOutput[bin]);
}
// Derivative-window STFT for synchrosqueezing
result->segments[seg].derivativeSpectrum = (FrequencyData*)malloc(numBins * sizeof(FrequencyData));
for (int i = 0; i < fftSize; i++) complexInput[i] = derivWindowedSamples[i] + 0.0f * I;
FFT(complexInput, fftOutput, fftSize, false);
for (int bin = 0; bin < numBins; bin++) {
result->segments[seg].derivativeSpectrum[bin].frequency = (float)bin * signal->sampleRate / fftSize;
result->segments[seg].derivativeSpectrum[bin].amplitude = cabsf(fftOutput[bin]) / fftSize;
result->segments[seg].derivativeSpectrum[bin].phase = cargf(fftOutput[bin]);
}
} }
free(windowedSamples); // Hann window h(t) = 0.5*(1 - cos(2πt)); derivative window h'(t) = π*sin(2πt)
free(derivWindowedSamples); for (int i = 0; i < fftSize; i++) {
free(complexInput); float t = (float)i / (fftSize - 1);
free(fftOutput); sc->windowed[i] *= 0.5f * (1.0f - cosf(2.0f * M_PI * t));
sc->derivWindowed[i] *= M_PI * sinf(2.0f * M_PI * t);
}
result->segments[seg].numBins = numBins;
result->segments[seg].sampleOffset = offset;
result->segments[seg].sampleCount = samplesToCopy;
// Normal STFT (V_f)
for (int i = 0; i < fftSize; i++) sc->fftIn[i] = sc->windowed[i] + 0.0f * I;
FFT(sc->fftIn, sc->fftOut, fftSize, false);
result->segments[seg].spectrum = (FrequencyData*)malloc(numBins * sizeof(FrequencyData));
for (int bin = 0; bin < numBins; bin++) {
result->segments[seg].spectrum[bin].frequency = (float)bin * signal->sampleRate / fftSize;
result->segments[seg].spectrum[bin].amplitude = (bin == 0) ? cabsf(sc->fftOut[bin]) / fftSize : 2.0f * cabsf(sc->fftOut[bin]) / fftSize;
result->segments[seg].spectrum[bin].phase = cargf(sc->fftOut[bin]);
}
// Derivative-window STFT (V_fd) for synchrosqueezing
for (int i = 0; i < fftSize; i++) sc->fftIn[i] = sc->derivWindowed[i] + 0.0f * I;
FFT(sc->fftIn, sc->fftOut, fftSize, false);
result->segments[seg].derivativeSpectrum = (FrequencyData*)malloc(numBins * sizeof(FrequencyData));
for (int bin = 0; bin < numBins; bin++) {
result->segments[seg].derivativeSpectrum[bin].frequency = (float)bin * signal->sampleRate / fftSize;
result->segments[seg].derivativeSpectrum[bin].amplitude = cabsf(sc->fftOut[bin]) / fftSize;
result->segments[seg].derivativeSpectrum[bin].phase = cargf(sc->fftOut[bin]);
}
}
// ===== Background high-res computation =====
// Fill segments [startSeg, endSeg) at full resolution, skipping any already
// computed. Returns the next segment index to resume from.
int ComputeNextHighResChunk(AudioSignal* signal, StftResult* result,
int fftSize, int startSeg, int endSeg)
{
SegScratch sc = AllocSegScratch(fftSize);
for (int seg = startSeg; seg < endSeg && seg < result->numSegments; seg++) {
if (result->segments[seg].spectrum != NULL) continue; // already computed
ComputeSegment(signal, result, fftSize, seg, &sc);
}
FreeSegScratch(&sc);
// Return next segment to process
if (endSeg >= result->numSegments) return result->numSegments; if (endSeg >= result->numSegments) return result->numSegments;
return endSeg; return endSeg;
} }
@@ -224,71 +247,13 @@ void ComputeSTFTInit(AudioSignal* signal, StftResult* result, int fftSize)
bool ComputeSTFTIncremental(AudioSignal* signal, StftResult* result, int fftSize, int startSegment) bool ComputeSTFTIncremental(AudioSignal* signal, StftResult* result, int fftSize, int startSegment)
{ {
int hopSize = fftSize / HOP_RATIO; SegScratch sc = AllocSegScratch(fftSize);
int numBins = fftSize / 2 + 1;
float* windowedSamples = (float*)malloc(fftSize * sizeof(float));
float* derivWindowedSamples = (float*)malloc(fftSize * sizeof(float));
float complex *complexInput = (float complex*)malloc(fftSize * sizeof(float complex));
float complex* fftOutput = (float complex*)malloc(fftSize * sizeof(float complex));
for (int seg = startSegment; seg < result->numSegments; seg++) { for (int seg = startSegment; seg < result->numSegments; seg++) {
// Skip segments not aligned with the skip factor (overview mode) if (seg % app.skipFactor != 0) continue; // overview stride
if (seg % app.skipFactor != 0) continue; if (result->segments[seg].spectrum != NULL) continue; // already computed
ComputeSegment(signal, result, fftSize, seg, &sc);
// Skip if already computed as high-res
if (result->segments[seg].spectrum != NULL) continue;
int offset = seg * hopSize;
int samplesToCopy = fftSize;
if (offset + samplesToCopy > signal->numSamples) {
samplesToCopy = signal->numSamples - offset;
memset(windowedSamples, 0, fftSize * sizeof(float));
memset(derivWindowedSamples, 0, fftSize * sizeof(float));
} else {
memcpy(windowedSamples, signal->samples + offset, fftSize * sizeof(float));
memcpy(derivWindowedSamples, signal->samples + offset, fftSize * sizeof(float));
}
// Apply Hann window: h(t) = 0.5 * (1 - cos(2πt))
// And derivative window: h'(t) = π * sin(2πt)
for (int i = 0; i < fftSize; i++) {
float t = (float)i / (fftSize - 1);
float hann = 0.5f * (1.0f - cosf(2.0f * M_PI * t));
float derivHann = M_PI * sinf(2.0f * M_PI * t);
windowedSamples[i] *= hann;
derivWindowedSamples[i] *= derivHann;
}
// Compute normal STFT (V_f)
for (int i = 0; i < fftSize; i++) complexInput[i] = windowedSamples[i] + 0.0f * I;
FFT(complexInput, fftOutput, fftSize, false);
result->segments[seg].numBins = numBins;
result->segments[seg].sampleOffset = offset;
result->segments[seg].sampleCount = samplesToCopy;
result->segments[seg].spectrum = (FrequencyData*)malloc(numBins * sizeof(FrequencyData));
for (int bin = 0; bin < numBins; bin++) {
result->segments[seg].spectrum[bin].frequency = (float)bin * signal->sampleRate / fftSize;
result->segments[seg].spectrum[bin].amplitude = (bin == 0) ? cabsf(fftOutput[bin]) / fftSize : 2.0f * cabsf(fftOutput[bin]) / fftSize;
result->segments[seg].spectrum[bin].phase = cargf(fftOutput[bin]);
}
// Compute derivative-window STFT (V_fd) for synchrosqueezing
result->segments[seg].derivativeSpectrum = (FrequencyData*)malloc(numBins * sizeof(FrequencyData));
for (int i = 0; i < fftSize; i++) complexInput[i] = derivWindowedSamples[i] + 0.0f * I;
FFT(complexInput, fftOutput, fftSize, false);
for (int bin = 0; bin < numBins; bin++) {
result->segments[seg].derivativeSpectrum[bin].frequency = (float)bin * signal->sampleRate / fftSize;
result->segments[seg].derivativeSpectrum[bin].amplitude = cabsf(fftOutput[bin]) / fftSize;
result->segments[seg].derivativeSpectrum[bin].phase = cargf(fftOutput[bin]);
}
} }
FreeSegScratch(&sc);
free(windowedSamples);
free(derivWindowedSamples);
free(complexInput);
free(fftOutput);
return true; return true;
} }