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

Remove reassignment code - back to standard spectrogram

Browse Source 
Co-authored-by: Qwen-Coder <qwen-coder@alibabacloud.com>
This commit is contained in:
Tyler
2026-03-30 19:40:20 -07:00
parent c00c665226
commit bbd476c389
1 changed files with 24 additions and 102 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/spectrogram.c
+24 -102
View File
@@ -34,9 +34,6 @@
#define MAX_SAMPLE_RATE 48000 #define MAX_SAMPLE_RATE 48000
#define LOUDNESS_FLOOR_DB -80.0f #define LOUDNESS_FLOOR_DB -80.0f
// Reassignment method for sharper time-frequency localization
#define USE_REASSIGNMENT 1
// Colormap types // Colormap types
typedef enum { typedef enum {
COLORMAP_GRAYS = 0, COLORMAP_GRAYS = 0,
@@ -270,34 +267,44 @@ static void ComputeSTFT(AudioSignal* signal, StftResult* result, int fftSize)
int hopSize = fftSize / HOP_RATIO; // 75% overlap int hopSize = fftSize / HOP_RATIO; // 75% overlap
int numSegments = (signal->numSamples - fftSize) / hopSize + 1; int numSegments = (signal->numSamples - fftSize) / hopSize + 1;
if (numSegments <= 0) numSegments = 1; if (numSegments <= 0) numSegments = 1;
result->numSegments = numSegments; result->numSegments = numSegments;
result->segments = (StftSegment*)malloc(numSegments * sizeof(StftSegment)); result->segments = (StftSegment*)malloc(numSegments * sizeof(StftSegment));
result->sampleRate = signal->sampleRate; result->sampleRate = signal->sampleRate;
result->totalSamples = signal->numSamples; result->totalSamples = signal->numSamples;
result->useHannWindow = true; result->useHannWindow = true;
int numBins = fftSize / 2 + 1; int numBins = fftSize / 2 + 1;
float* windowedSamples = (float*)malloc(fftSize * sizeof(float)); float* windowedSamples = (float*)malloc(fftSize * sizeof(float));
float complex *complexInput = (float complex*)malloc(fftSize * sizeof(float complex)); float complex *complexInput = (float complex*)malloc(fftSize * sizeof(float complex));
float complex* fftOutput = (float complex*)malloc(fftSize * sizeof(float complex)); float complex* fftOutput = (float complex*)malloc(fftSize * sizeof(float complex));
for (int seg = 0; seg < numSegments; seg++) { for (int seg = 0; seg < numSegments; seg++) {
int offset = seg * hopSize; int offset = seg * hopSize;
int samplesToCopy = fftSize; int samplesToCopy = fftSize;
if (offset + samplesToCopy > signal->numSamples) { if (offset + samplesToCopy > signal->numSamples) {
samplesToCopy = signal->numSamples - offset; samplesToCopy = signal->numSamples - offset;
memset(windowedSamples, 0, fftSize * sizeof(float)); 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; for (int i = 0; i < fftSize; i++) complexInput[i] = windowedSamples[i] + 0.0f * I;
FFT(complexInput, fftOutput, fftSize, false); FFT(complexInput, fftOutput, fftSize, false);
result->segments[seg].numBins = numBins; result->segments[seg].numBins = numBins;
result->segments[seg].sampleOffset = offset; result->segments[seg].sampleOffset = offset;
result->segments[seg].sampleCount = samplesToCopy; result->segments[seg].sampleCount = samplesToCopy;
result->segments[seg].spectrum = (FrequencyData*)malloc(numBins * sizeof(FrequencyData)); result->segments[seg].spectrum = (FrequencyData*)malloc(numBins * sizeof(FrequencyData));
for (int bin = 0; bin < numBins; bin++) { for (int bin = 0; bin < numBins; bin++) {
result->segments[seg].spectrum[bin].frequency = (float)bin * signal->sampleRate / fftSize; 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].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) 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); free(result->segments);
result->segments = NULL; result->segments = NULL;
result->numSegments = 0; 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) 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; int height = stft->segments[0].numBins;
*image = GenImageColor(width, height, BLACK); *image = GenImageColor(width, height, BLACK);
Color* pixels = (Color*)image->data; Color* pixels = (Color*)image->data;
// Initialize to black
for (int i = 0; i < width * height; i++) {
pixels[i] = BLACK;
}
// Find max amplitude for normalization // Find max amplitude for normalization
float maxAmplitude = 0.0001f; float maxAmplitude = 0.0001f;
@@ -385,95 +389,16 @@ static void GenerateSpectrogramTexture(StftResult* stft, Image* image, Texture2D
if (stft->segments[seg].spectrum[bin].amplitude > maxAmplitude) if (stft->segments[seg].spectrum[bin].amplitude > maxAmplitude)
maxAmplitude = stft->segments[seg].spectrum[bin].amplitude; 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++) { for (int seg = 0; seg < width; seg++) {
float segmentTime = (seg * hopSize) / (float)stft->sampleRate;
for (int bin = 0; bin < height; bin++) { for (int bin = 0; bin < height; bin++) {
float amplitude = stft->segments[seg].spectrum[bin].amplitude; float amplitude = stft->segments[seg].spectrum[bin].amplitude;
float phase = stft->segments[seg].spectrum[bin].phase; float db = AmplitudeToDecibels(amplitude);
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 normalized = (db - app.amplitudeFloorDb) / (app.amplitudeCeilingDb - app.amplitudeFloorDb); float normalized = (db - app.amplitudeFloorDb) / (app.amplitudeCeilingDb - app.amplitudeFloorDb);
normalized = Clamp(normalized, 0.0f, 1.0f); 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); if (texture->id != 0) UnloadTexture(*texture);
*texture = LoadTextureFromImage(*image); *texture = LoadTextureFromImage(*image);
SetTextureFilter(*texture, TEXTURE_FILTER_BILINEAR); 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("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("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("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; DrawText(TextFormat("FFT: %d (%.1f Hz/bin)", 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("Max Freq: %.1f kHz", (float)app.signal.sampleRate / 2000.0f), 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; y += 10;