diff --git a/src/stft.c b/src/stft.c index 14bca63..ad963ea 100644 --- a/src/stft.c +++ b/src/stft.c @@ -135,74 +135,97 @@ void SaveToCache(void) app.fftSize, app.stft.numSegments); } -// ===== Background high-res computation ===== -int ComputeNextHighResChunk(AudioSignal* signal, StftResult* result, - int fftSize, int startSeg, int endSeg) +// ===== Per-segment STFT (shared by the overview and high-res passes) ===== + +// 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 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 = startSeg; seg < endSeg && seg < result->numSegments; seg++) { - // Skip if already computed (overview or 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 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]); - } + int offset = seg * hopSize; + int samplesToCopy = fftSize; + if (offset + samplesToCopy > signal->numSamples) { + samplesToCopy = signal->numSamples - offset; + memset(sc->windowed, 0, fftSize * sizeof(float)); + memset(sc->derivWindowed, 0, fftSize * sizeof(float)); + } else { + memcpy(sc->windowed, signal->samples + offset, fftSize * sizeof(float)); + memcpy(sc->derivWindowed, signal->samples + offset, fftSize * sizeof(float)); } - free(windowedSamples); - free(derivWindowedSamples); - free(complexInput); - free(fftOutput); + // Hann window h(t) = 0.5*(1 - cos(2πt)); derivative window h'(t) = π*sin(2πt) + for (int i = 0; i < fftSize; i++) { + float t = (float)i / (fftSize - 1); + 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; return endSeg; } @@ -224,71 +247,13 @@ void ComputeSTFTInit(AudioSignal* signal, StftResult* result, int fftSize) bool ComputeSTFTIncremental(AudioSignal* signal, StftResult* result, int fftSize, int startSegment) { - int hopSize = fftSize / HOP_RATIO; - 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)); - + SegScratch sc = AllocSegScratch(fftSize); for (int seg = startSegment; seg < result->numSegments; seg++) { - // Skip segments not aligned with the skip factor (overview mode) - if (seg % app.skipFactor != 0) continue; - - // 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]); - } + if (seg % app.skipFactor != 0) continue; // overview stride + if (result->segments[seg].spectrum != NULL) continue; // already computed + ComputeSegment(signal, result, fftSize, seg, &sc); } - - free(windowedSamples); - free(derivWindowedSamples); - free(complexInput); - free(fftOutput); + FreeSegScratch(&sc); return true; }