diff --git a/src/spectrogram.c b/src/spectrogram.c index 9567e0c..e3f62af 100644 --- a/src/spectrogram.c +++ b/src/spectrogram.c @@ -27,9 +27,7 @@ // Configuration // ============================================================================ -#define FFT_SIZE_DEFAULT 2048 -#define FFT_SIZE_MAX 2048 -#define FFT_SIZE_MIN 128 +#define FFT_SIZE 128 #define HOP_RATIO 4 // FFT_SIZE / HOP_SIZE = 4 means 75% overlap #define MAX_SAMPLE_RATE 48000 #define LOUDNESS_FLOOR_DB -80.0f @@ -116,8 +114,6 @@ typedef struct { float amplitudeCeilingDb; ColormapType colormap; bool showGrid; - int fftSize; // Current FFT size (128-2048) - bool useSynchrosqueezing; // Enable synchrosqueezing for sharper display // File browser state bool showFileBrowser; @@ -466,96 +462,81 @@ static void GenerateSpectrogramTexture(StftResult* stft, Image* image, Texture2D if (stft->segments[seg].spectrum[bin].amplitude > maxAmplitude) maxAmplitude = stft->segments[seg].spectrum[bin].amplitude; - if (app.useSynchrosqueezing) { - // ===== SYNCHROSQUEEZING ===== - // Reassign energy to true frequencies using derivative STFT - - // Accumulation buffer for reassigned energy - float* accumBuffer = (float*)calloc(width * height, sizeof(float)); - - // Noise threshold: only reassign bins with significant energy - // This prevents noise from being reassigned to random frequencies - float noiseThreshold = maxAmplitude * 0.01f; // 1% of max amplitude - - for (int seg = 0; seg < width; seg++) { - for (int bin = 0; bin < height; bin++) { - FrequencyData* V_f = &stft->segments[seg].spectrum[bin]; - FrequencyData* V_fd = &stft->segments[seg].derivativeSpectrum[bin]; - - float amplitude = V_f->amplitude; - - // Skip noise bins - if (amplitude < noiseThreshold) continue; - - // Compute instantaneous frequency using synchrosqueezing formula: - // ω̂ = bin_freq + Re[V_fd / (i * V_f)] - // Complex division: (a+bi)/(c+di) = ((ac+bd) + (bc-ad)i) / (c²+d²) - // We need Re[(a+bi) / (i*(c+di))] = Re[(a+bi) / (-d+ci)] = (ad+bc)/(c²+d²) - - float V_f_real = amplitude * cosf(V_f->phase); - float V_f_imag = amplitude * sinf(V_f->phase); - float V_fd_real = V_fd->amplitude * cosf(V_fd->phase); - float V_fd_imag = V_fd->amplitude * sinf(V_fd->phase); - - float denom = V_f_real * V_f_real + V_f_imag * V_f_imag; - - float trueFreq = V_f->frequency; // Default to bin frequency - - if (denom > 1e-10f) { - // Re[V_fd / (i * V_f)] = (-V_fd_real * V_f_imag + V_fd_imag * V_f_real) / denom - // Note the MINUS sign on the first term - float correction = (-V_fd_real * V_f_imag + V_fd_imag * V_f_real) / denom; - trueFreq = V_f->frequency + correction; - } - - // Clamp to valid range - if (trueFreq < 0) trueFreq = 0; - if (trueFreq >= stft->sampleRate / 2.0f) trueFreq = stft->sampleRate / 2.0f - 1; - - // Map to bin coordinate - float targetBinF = trueFreq / freqPerBin; - if (targetBinF < 0) targetBinF = 0; - if (targetBinF >= height) targetBinF = height - 0.001f; - - // Bilinear splatting to neighboring bins - int bin0 = (int)targetBinF; - int bin1 = bin0 + 1; - if (bin1 >= height) bin1 = height - 1; - - float frac = targetBinF - bin0; - int idx0 = (height - 1 - bin0) * width + seg; - int idx1 = (height - 1 - bin1) * width + seg; - - accumBuffer[idx0] += amplitude * (1 - frac); - accumBuffer[idx1] += amplitude * frac; - } - } - - // Convert accumulation buffer to colors - for (int i = 0; i < width * height; i++) { - if (accumBuffer[i] > 0.0001f) { - float db = AmplitudeToDecibels(accumBuffer[i]); - float normalized = (db - app.amplitudeFloorDb) / (app.amplitudeCeilingDb - app.amplitudeFloorDb); - normalized = Clamp(normalized, 0.0f, 1.0f); - pixels[i] = GetColormapColor(normalized, app.colormap); - } - } - - free(accumBuffer); - } else { - // ===== STANDARD STFT (no synchrosqueezing) ===== - for (int seg = 0; seg < width; seg++) { - for (int bin = 0; bin < height; bin++) { - float amplitude = stft->segments[seg].spectrum[bin].amplitude; - float db = AmplitudeToDecibels(amplitude); - float normalized = (db - app.amplitudeFloorDb) / (app.amplitudeCeilingDb - app.amplitudeFloorDb); - normalized = Clamp(normalized, 0.0f, 1.0f); - int pixelIndex = (height - 1 - bin) * width + seg; - pixels[pixelIndex] = GetColormapColor(normalized, app.colormap); + // ===== SYNCHROSQUEEZING ===== + // Reassign energy to true frequencies using derivative STFT + + // Accumulation buffer for reassigned energy + float* accumBuffer = (float*)calloc(width * height, sizeof(float)); + + // Noise threshold: only reassign bins with significant energy + float noiseThreshold = maxAmplitude * 0.01f; // 1% of max amplitude + + for (int seg = 0; seg < width; seg++) { + for (int bin = 0; bin < height; bin++) { + FrequencyData* V_f = &stft->segments[seg].spectrum[bin]; + FrequencyData* V_fd = &stft->segments[seg].derivativeSpectrum[bin]; + + float amplitude = V_f->amplitude; + + // Skip noise bins + if (amplitude < noiseThreshold) continue; + + // Compute instantaneous frequency using synchrosqueezing formula: + // ω̂ = bin_freq + Re[V_fd / (i * V_f)] + // Complex division: (a+bi)/(c+di) = ((ac+bd) + (bc-ad)i) / (c²+d²) + // We need Re[(a+bi) / (i*(c+di))] = Re[(a+bi) / (-d+ci)] = (ad+bc)/(c²+d²) + + float V_f_real = amplitude * cosf(V_f->phase); + float V_f_imag = amplitude * sinf(V_f->phase); + float V_fd_real = V_fd->amplitude * cosf(V_fd->phase); + float V_fd_imag = V_fd->amplitude * sinf(V_fd->phase); + + float denom = V_f_real * V_f_real + V_f_imag * V_f_imag; + + float trueFreq = V_f->frequency; // Default to bin frequency + + if (denom > 1e-10f) { + // Re[V_fd / (i * V_f)] = (-V_fd_real * V_f_imag + V_fd_imag * V_f_real) / denom + // Note the MINUS sign on the first term + float correction = (-V_fd_real * V_f_imag + V_fd_imag * V_f_real) / denom; + trueFreq = V_f->frequency + correction; } + + // Clamp to valid range + if (trueFreq < 0) trueFreq = 0; + if (trueFreq >= stft->sampleRate / 2.0f) trueFreq = stft->sampleRate / 2.0f - 1; + + // Map to bin coordinate + float targetBinF = trueFreq / freqPerBin; + if (targetBinF < 0) targetBinF = 0; + if (targetBinF >= height) targetBinF = height - 0.001f; + + // Bilinear splatting to neighboring bins + int bin0 = (int)targetBinF; + int bin1 = bin0 + 1; + if (bin1 >= height) bin1 = height - 1; + + float frac = targetBinF - bin0; + int idx0 = (height - 1 - bin0) * width + seg; + int idx1 = (height - 1 - bin1) * width + seg; + + accumBuffer[idx0] += amplitude * (1 - frac); + accumBuffer[idx1] += amplitude * frac; } } + // Convert accumulation buffer to colors + for (int i = 0; i < width * height; i++) { + if (accumBuffer[i] > 0.0001f) { + float db = AmplitudeToDecibels(accumBuffer[i]); + float normalized = (db - app.amplitudeFloorDb) / (app.amplitudeCeilingDb - app.amplitudeFloorDb); + normalized = Clamp(normalized, 0.0f, 1.0f); + pixels[i] = GetColormapColor(normalized, app.colormap); + } + } + + free(accumBuffer); + if (texture->id != 0) UnloadTexture(*texture); *texture = LoadTextureFromImage(*image); SetTextureFilter(*texture, TEXTURE_FILTER_BILINEAR); @@ -988,24 +969,7 @@ static void DrawSidebar(void) // Title DrawText("Spectrogram Controls", x, y, 14, WHITE); y += 25; - - // FFT Size slider - DrawText(TextFormat("FFT Size: %d (%.1f Hz/bin)", app.fftSize, (float)app.signal.sampleRate / app.fftSize), x, y, fontSize, LIGHTGRAY); y += 18; - Rectangle fftSlider = { x, y, sidebarWidth - 10, 20 }; - float fftValue = (float)(app.fftSize - FFT_SIZE_MIN) / (FFT_SIZE_MAX - FFT_SIZE_MIN); - DrawSlider(fftSlider, fftValue); - if (UpdateSlider(fftSlider, &fftValue)) { - int newFFT = FFT_SIZE_MIN + (int)(fftValue * (FFT_SIZE_MAX - FFT_SIZE_MIN)); - // Round to power of 2 - while ((newFFT & (newFFT - 1)) != 0) newFFT++; - if (newFFT != app.fftSize && newFFT >= FFT_SIZE_MIN && newFFT <= FFT_SIZE_MAX) { - app.fftSize = newFFT; - app.stftComputed = false; - needsRegen = true; - } - } - y += 25; - + // dB Floor slider DrawText(TextFormat("dB Floor: %.1f", app.amplitudeFloorDb), x, y, fontSize, LIGHTGRAY); y += 18; Rectangle dbSlider = { x, y, sidebarWidth - 10, 20 }; @@ -1016,7 +980,7 @@ static void DrawSidebar(void) needsRegen = true; // Trigger immediate redraw } y += 25; - + // Colormap dropdown DrawText("Colormap:", x, y, fontSize, LIGHTGRAY); y += 18; const char* colormapNames[] = { "Grays", "Inferno", "Viridis", "Plasma", "Hot", "Cool" }; @@ -1029,11 +993,11 @@ static void DrawSidebar(void) DrawRectangleRec(cmapButton, (Color){ 50, 50, 60, 255 }); DrawRectangleLinesEx(cmapButton, 1, GRAY); DrawText(colormapNames[app.colormap], cmapButton.x + 10, cmapButton.y + 6, fontSize, WHITE); - DrawTexturePro(colormapTexture, (Rectangle){ 0, 0, 256, 1 }, + DrawTexturePro(colormapTexture, (Rectangle){ 0, 0, 256, 1 }, (Rectangle){ cmapButton.x + cmapButton.width - 60, cmapButton.y + 5, 50, 15 }, (Vector2){ 0, 0 }, 0.0f, WHITE); y += 30; - + // Grid toggle Rectangle gridCheck = { x, y, 18, 18 }; if (CheckCollisionPointRec(GetMousePosition(), gridCheck) && IsMouseButtonPressed(MOUSE_LEFT_BUTTON)) { @@ -1042,17 +1006,7 @@ static void DrawSidebar(void) DrawRectangleRec(gridCheck, app.showGrid ? BLUE : DARKGRAY); DrawRectangleLinesEx(gridCheck, 1, WHITE); DrawText("Show Grid", x + 25, y + 2, fontSize, LIGHTGRAY); y += 25; - - // Synchrosqueezing toggle - Rectangle sqCheck = { x, y, 18, 18 }; - if (CheckCollisionPointRec(GetMousePosition(), sqCheck) && IsMouseButtonPressed(MOUSE_LEFT_BUTTON)) { - app.useSynchrosqueezing = !app.useSynchrosqueezing; - needsRegen = true; - } - DrawRectangleRec(sqCheck, app.useSynchrosqueezing ? BLUE : DARKGRAY); - DrawRectangleLinesEx(sqCheck, 1, WHITE); - DrawText("Synchrosqueezing (sharp)", x + 25, y + 2, fontSize, LIGHTGRAY); y += 25; - + // File loading DrawText("File:", x, y, fontSize, LIGHTGRAY); y += 18; Rectangle fileButton = { x, y, sidebarWidth - 10, 25 }; @@ -1111,7 +1065,7 @@ static void DrawSidebar(void) if (app.loaded) { DrawText(TextFormat("Sample Rate: %d Hz", app.signal.sampleRate), x, y, fontSize, GRAY); y += 16; DrawText(TextFormat("Duration: %.2f sec", app.signal.duration), x, y, fontSize, GRAY); y += 16; - DrawText(TextFormat("Max Freq: %.1f kHz", (float)app.signal.sampleRate / 2000.0f), x, y, fontSize, GRAY); y += 16; + DrawText(TextFormat("FFT: %d, %.1f Hz/bin (synchrosqueezed)", FFT_SIZE, (float)app.signal.sampleRate / FFT_SIZE), x, y, fontSize, GRAY); y += 16; } else { DrawText("No file loaded", x, y, fontSize, GRAY); y += 16; } @@ -1177,10 +1131,8 @@ int main(int argc, char* argv[]) app.cachedVisibleStart = -1; app.cachedVisibleEnd = -1; app.visibleTextureValid = false; - app.fftSize = FFT_SIZE_DEFAULT; app.isPlaying = false; app.playbackFinished = false; - app.useSynchrosqueezing = true; // Enabled by default GenerateColormapTexture(); ScanDirectory(GetWorkingDirectory()); @@ -1396,7 +1348,7 @@ int main(int argc, char* argv[]) if (app.loaded && !app.stftComputed) { TraceLog(LOG_INFO, "Computing STFT..."); double startTime = GetTime(); - ComputeSTFT(&app.signal, &app.stft, app.fftSize); + ComputeSTFT(&app.signal, &app.stft, FFT_SIZE); TraceLog(LOG_INFO, "STFT computed in %.2f sec (%d segments)", GetTime() - startTime, app.stft.numSegments); GenerateSpectrogramTexture(&app.stft, &app.spectrogramImage, &app.spectrogramTexture); app.stftComputed = true;