From bbd476c3896c83cbeccc6cfed78e27220fca4d15 Mon Sep 17 00:00:00 2001 From: Tyler Date: Mon, 30 Mar 2026 19:40:20 -0700 Subject: [PATCH] Remove reassignment code - back to standard spectrogram Co-authored-by: Qwen-Coder --- src/spectrogram.c | 126 +++++++++------------------------------------- 1 file changed, 24 insertions(+), 102 deletions(-) diff --git a/src/spectrogram.c b/src/spectrogram.c index c1bbdbf..16c8d33 100644 --- a/src/spectrogram.c +++ b/src/spectrogram.c @@ -34,9 +34,6 @@ #define MAX_SAMPLE_RATE 48000 #define LOUDNESS_FLOOR_DB -80.0f -// Reassignment method for sharper time-frequency localization -#define USE_REASSIGNMENT 1 - // Colormap types typedef enum { COLORMAP_GRAYS = 0, @@ -270,34 +267,44 @@ static void ComputeSTFT(AudioSignal* signal, StftResult* result, int fftSize) int hopSize = fftSize / HOP_RATIO; // 75% overlap int numSegments = (signal->numSamples - fftSize) / hopSize + 1; if (numSegments <= 0) numSegments = 1; - + result->numSegments = numSegments; result->segments = (StftSegment*)malloc(numSegments * sizeof(StftSegment)); result->sampleRate = signal->sampleRate; result->totalSamples = signal->numSamples; result->useHannWindow = true; - + int numBins = fftSize / 2 + 1; float* windowedSamples = (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 = 0; seg < numSegments; seg++) { int offset = seg * hopSize; int samplesToCopy = fftSize; if (offset + samplesToCopy > signal->numSamples) { samplesToCopy = signal->numSamples - offset; memset(windowedSamples, 0, fftSize * sizeof(float)); + } else { + memcpy(windowedSamples, signal->samples + offset, fftSize * sizeof(float)); } - memcpy(windowedSamples, signal->samples + offset, samplesToCopy * sizeof(float)); - ApplyHannWindow(windowedSamples, fftSize); + + // Apply Hann window: h(t) = 0.5 * (1 - cos(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)); + windowedSamples[i] *= hann; + } + + // Compute 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; @@ -309,7 +316,9 @@ static void ComputeSTFT(AudioSignal* signal, StftResult* result, int fftSize) static void FreeSTFT(StftResult* result) { - for (int i = 0; i < result->numSegments; i++) free(result->segments[i].spectrum); + for (int i = 0; i < result->numSegments; i++) { + free(result->segments[i].spectrum); + } free(result->segments); result->segments = NULL; result->numSegments = 0; @@ -362,7 +371,7 @@ static void FreeSignal(AudioSignal* signal) } // ============================================================================ -// Spectrogram Generation with Reassignment +// Spectrogram Generation // ============================================================================ static void GenerateSpectrogramTexture(StftResult* stft, Image* image, Texture2D* texture) @@ -372,11 +381,6 @@ static void GenerateSpectrogramTexture(StftResult* stft, Image* image, Texture2D int height = stft->segments[0].numBins; *image = GenImageColor(width, height, BLACK); Color* pixels = (Color*)image->data; - - // Initialize to black - for (int i = 0; i < width * height; i++) { - pixels[i] = BLACK; - } // Find max amplitude for normalization float maxAmplitude = 0.0001f; @@ -385,95 +389,16 @@ static void GenerateSpectrogramTexture(StftResult* stft, Image* image, Texture2D if (stft->segments[seg].spectrum[bin].amplitude > maxAmplitude) maxAmplitude = stft->segments[seg].spectrum[bin].amplitude; - float hopSize = (float)(stft->segments[0].sampleOffset > 0 ? - stft->segments[1].sampleOffset - stft->segments[0].sampleOffset : - stft->segments[0].sampleCount); - float freqPerBin = (float)stft->sampleRate / (height * 2 - 2); - - // Create a floating-point accumulation buffer for reassignment - float* accumBuffer = (float*)calloc(width * height, sizeof(float)); - for (int seg = 0; seg < width; seg++) { - float segmentTime = (seg * hopSize) / (float)stft->sampleRate; - for (int bin = 0; bin < height; bin++) { float amplitude = stft->segments[seg].spectrum[bin].amplitude; - float phase = stft->segments[seg].spectrum[bin].phase; - - if (amplitude < 0.0001f) continue; - - float reassignedFreq = bin * freqPerBin; - float reassignedTime = segmentTime; - -#if USE_REASSIGNMENT - // ===== Reassignment Method ===== - // Estimate instantaneous frequency from phase derivative over time - if (seg > 0 && seg < width - 1) { - float prevPhase = stft->segments[seg-1].spectrum[bin].phase; - float nextPhase = stft->segments[seg+1].spectrum[bin].phase; - - // Phase difference (unwrapped) - float phaseDiff = nextPhase - prevPhase; - // Unwrap to [-pi, pi] - while (phaseDiff > M_PI) phaseDiff -= 2.0f * M_PI; - while (phaseDiff < -M_PI) phaseDiff += 2.0f * M_PI; - - // Instantaneous frequency deviation - float expectedPhaseShift = 2.0f * M_PI * bin * hopSize / (height * 2 - 2); - float phaseDeviation = phaseDiff - expectedPhaseShift; - float instantFreqDev = phaseDeviation * stft->sampleRate / (2.0f * M_PI * hopSize); - - reassignedFreq = bin * freqPerBin + instantFreqDev; - } - - // Estimate group delay from phase derivative over frequency - if (bin > 0 && bin < height - 1) { - float prevPhaseAdj = stft->segments[seg].spectrum[bin-1].phase; - float nextPhaseAdj = stft->segments[seg].spectrum[bin+1].phase; - - float phaseGrad = nextPhaseAdj - prevPhaseAdj; - // Unwrap - while (phaseGrad > M_PI) phaseGrad -= 2.0f * M_PI; - while (phaseGrad < -M_PI) phaseGrad += 2.0f * M_PI; - - // Group delay (time correction) - float groupDelay = -phaseGrad / (2.0f * M_PI * freqPerBin); - reassignedTime = segmentTime + groupDelay / (float)stft->sampleRate; - } - - // Clamp to valid range - if (reassignedFreq < 0) reassignedFreq = 0; - if (reassignedFreq >= stft->sampleRate / 2.0f) reassignedFreq = stft->sampleRate / 2.0f - 1; -#endif - - // Map reassigned coordinates to pixel indices - int reassignedBin = (int)(reassignedFreq / freqPerBin); - int reassignedSeg = (int)((reassignedTime * stft->sampleRate) / hopSize); - - // Clamp to texture bounds - if (reassignedBin < 0) reassignedBin = 0; - if (reassignedBin >= height) reassignedBin = height - 1; - if (reassignedSeg < 0) reassignedSeg = 0; - if (reassignedSeg >= width) reassignedSeg = width - 1; - - // Accumulate amplitude at reassigned location - int pixelIndex = (height - 1 - reassignedBin) * width + reassignedSeg; - accumBuffer[pixelIndex] += amplitude; - } - } - - // Convert accumulation buffer to colors - for (int i = 0; i < width * height; i++) { - if (accumBuffer[i] > 0.0001f) { - float db = AmplitudeToDecibels(accumBuffer[i]); + float db = AmplitudeToDecibels(amplitude); float normalized = (db - app.amplitudeFloorDb) / (app.amplitudeCeilingDb - app.amplitudeFloorDb); normalized = Clamp(normalized, 0.0f, 1.0f); - pixels[i] = GetColormapColor(normalized, app.colormap); + int pixelIndex = (height - 1 - bin) * width + seg; + pixels[pixelIndex] = GetColormapColor(normalized, app.colormap); } } - - free(accumBuffer); - if (texture->id != 0) UnloadTexture(*texture); *texture = LoadTextureFromImage(*image); SetTextureFilter(*texture, TEXTURE_FILTER_BILINEAR); @@ -893,10 +818,7 @@ static void DrawInfo(Rectangle bounds) DrawText(TextFormat("Sample Rate: %d Hz", app.signal.sampleRate), 10, y, fontSize, LIGHTGRAY); y += 20; DrawText(TextFormat("Duration: %.2f sec", app.signal.duration), 10, y, fontSize, LIGHTGRAY); y += 20; DrawText(TextFormat("View: %.1f%%-%.1f%% (%.2f sec)", app.viewStart*100, app.viewEnd*100, (app.viewEnd-app.viewStart)*app.signal.duration), 10, y, fontSize, LIGHTGRAY); y += 20; - DrawText(TextFormat("FFT: %d (%.1f Hz/bin, 75%% overlap)", app.fftSize, (float)app.signal.sampleRate / app.fftSize), 10, y, fontSize, LIGHTGRAY); y += 20; -#if USE_REASSIGNMENT - DrawText("Reassignment: ON (sharp)", 10, y, fontSize, (Color){ 80, 255, 80, 255 }); y += 20; -#endif + DrawText(TextFormat("FFT: %d (%.1f Hz/bin)", app.fftSize, (float)app.signal.sampleRate / app.fftSize), 10, y, fontSize, LIGHTGRAY); y += 20; DrawText(TextFormat("Max Freq: %.1f kHz", (float)app.signal.sampleRate / 2000.0f), 10, y, fontSize, LIGHTGRAY); y += 20; y += 10;