Periodogram viewer

pg-viewer.c

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>

#include "SDL.h"

#include <fftw3.h>

#include <time.h>

struct draw_context
{
  fftwf_plan plan;
  
  float window[2048];
  float *window2;
  float *dft;
  
  float scale;
  float wfunc[2048];
  float dft_interp[854];
  float wv_interp[2048];
  
  SDL_FRect rects[854];
  SDL_FPoint points[2048];

  SDL_AudioDeviceID dev;
};

static inline void generate_hann(float* w, int N)
{
  int i;
  float a;

  for (i = 0; i < N; ++i)
  {
    a = sin (M_PI * (float)i / N);
    
    w[i] = a * a;
  }
}

static void capture_callback (void* userdata, unsigned char* samples, int num_bytes)
{
  struct draw_context *c = userdata;

  int N = num_bytes / sizeof(float);
  int win_off = 2048 - N;

  memmove (c->window, c->window + N, win_off * sizeof(float));
  memcpy (c->window + win_off, samples, num_bytes);
}

static inline void draw (SDL_Renderer *r, struct draw_context *c, int W, int H)
{
  SDL_FPoint point;
  SDL_FRect rect;
  
  int i;
  float x = 0;
  float mag;
  float y;
  float xstep = W / 2048.f;
  float xstep_dft = W / 854.f;
  
  SDL_SetRenderDrawColor(r, 0, 0, 0, 0xFF);
  SDL_RenderClear (r);

  SDL_LockAudioDevice (c->dev);

  for (i = 0; i < 2048; ++i)
      c->window2[i] = c->window[i] * c->wfunc[i] * c->scale;

  fftwf_execute (c->plan);

  for (i = 0; i < 854; ++i)
  {
    if (i == 0)
      mag = c->dft[i];
    else
      mag = hypotf (c->dft[i], c->dft[2048 - i]);

    /* Interpolated magnitude (tau = 0.8) */
    mag = 0.8 * c->dft_interp[i] + 0.2 * mag;
    
    c->dft_interp[i] = mag;

    mag /= 2048.0f;
    mag = 20.0f * log10(mag);

    /* Clamp operation. */
    mag = mag < -70.0f ? -70.0f : (mag > -10.0f ? -10.0f : mag);
    mag = (mag + 70.0f) / 60.0f;
    
    y = mag * H;
    
    rect.x = x;
    rect.w = xstep_dft;
    rect.y = 0.25f * (H - y);
    rect.h = 0.5f * y;
    
    c->rects[i] = rect;
    
    x += xstep_dft;
  }

  x = 0;

  for (i = 0; i < 2048; ++i)
  {
    y = c->window[i];

    point.x = x;
    point.y = H * (0.5f + .25f * y * c->scale + .25f);

    // TODO make a macro
    point.y = point.y < H/2.0 ? H/2.0 : (point.y > H ? H : point.y);
    
    c->points[i] = point;
 
    x += xstep;
  }

  SDL_UnlockAudioDevice (c->dev);
  
  SDL_SetRenderDrawColor (r, 0xFF, 0xFF, 0xFF, 0xFF);
  SDL_RenderDrawLine (r, 0, H/4, W, H/4);
  SDL_RenderFillRectsF (r, c->rects, 854);
  SDL_RenderDrawLinesF (r, c->points, 2048);
  SDL_SetRenderDrawColor (r, 0x00, 0x00, 0xFF, 0xFF);
  SDL_RenderDrawLine (r, 0, H/2, W, H/2);
}

int main(int argc, char ** argv)
{
  SDL_Init (SDL_INIT_AUDIO | SDL_INIT_VIDEO);
  
  SDL_Window *win = NULL;
  SDL_Renderer *r = NULL;
  SDL_Event evt;
  
  SDL_AudioDeviceID dev = 0;
  SDL_AudioSpec want, have;

  struct draw_context *ctx;

  ctx = calloc (sizeof (struct draw_context), 1);
  if (!ctx)
  {
    fputs ("Failed allocating `struct draw_context'", stderr);
    
    goto end;
  }

  ctx->scale = 1;
    
  int running = 1;
  int shown = 1;
  int fullscreen = 0;
  int width = 1280, height = 720;

  want.samples = 800;
  want.freq = 48000;
  want.format = AUDIO_F32;
  want.channels = 1;
  want.callback = capture_callback;
  want.userdata = ctx;

choose:
  int num_devs = SDL_GetNumAudioDevices(1);
  int dev_idx;

  for (int i = 0; i < num_devs; ++i)
	fprintf(stderr, "[%d] %s\n", i, SDL_GetAudioDeviceName(i, 1));
  fprintf(stderr, "> ");

  if (scanf ("%d", &dev_idx) != 1)
	goto choose;

  dev = SDL_OpenAudioDevice (SDL_GetAudioDeviceName(dev_idx, 1), 1, &want, &have, 0);
  
  if (!dev)
  {
    fprintf (stderr, "SDL_OpenAudioDevice(): %s\n", SDL_GetError ());
    
    goto end;
  }
  
  win = SDL_CreateWindow (argv[0],
    SDL_WINDOWPOS_UNDEFINED,
    SDL_WINDOWPOS_UNDEFINED,
    width, height,
    SDL_WINDOW_RESIZABLE);

  if (!win)
  {
    fprintf (stderr, "SDL_CreateWindow(): %s\n", SDL_GetError ());
    
    goto end;
  }

  r = SDL_CreateRenderer (win, -1, SDL_RENDERER_ACCELERATED | SDL_RENDERER_PRESENTVSYNC);

  if (!r)
  {
    fprintf (stderr, "SDL_CreateRenderer(): %s\n", SDL_GetError());

    goto end;  
  }

  ctx->window2 = fftwf_malloc (2048 * sizeof (float));
  ctx->dft = fftwf_malloc (2048 * sizeof (float));
  ctx->plan = fftwf_plan_r2r_1d (2048, ctx->window2, ctx->dft, FFTW_R2HC, 0);

  if (!ctx->plan)
  {
    fprintf (stderr, "Failed creating a fftw plan\n");
    goto end;
  }

  generate_hann (ctx->wfunc, 2048);

  SDL_PauseAudioDevice (dev, 0);

#ifdef NDEBUG
  /* For debugging purposes. */
  struct timespec Tp, Tb, Te;

  unsigned long Tavg = 0;
  unsigned long long iters = 0;
#endif

  while (running)
  {
#ifdef NDEBUG
    clock_gettime (CLOCK_MONOTONIC, &Tb);
#endif
    while (SDL_PollEvent(&evt))
      switch (evt.type)
      {
        case SDL_QUIT:
          running = 0;
        break;
        
        case SDL_KEYDOWN:
          switch (evt.key.keysym.sym)
          {
            case SDLK_ESCAPE:
              running = 0;
            break;

      	    case SDLK_f:
      	      if (!fullscreen)
      	      	SDL_SetWindowFullscreen (win, SDL_WINDOW_FULLSCREEN_DESKTOP);
      	      else
            		SDL_SetWindowFullscreen (win, 0);
      	      fullscreen = !fullscreen;
      	    break;

            case SDLK_UP:
              ctx->scale /= 0.9f;
            break;

            case SDLK_DOWN:
              ctx->scale *= 0.9f;
            break;

            default:
            break;
          }
        break;

        case SDL_WINDOWEVENT:
          switch (evt.window.event)
          {
            case SDL_WINDOWEVENT_HIDDEN:
              shown = 0;
            break;

            case SDL_WINDOWEVENT_SHOWN:
              shown = 1;
              break;

            case SDL_WINDOWEVENT_RESIZED:
              width  = evt.window.data1;
              height = evt.window.data2;
              break;
            default:
            break;
          }
        break;
        default:
        break;
      }

    if (shown)
      draw (r, ctx, width, height);

#ifdef NDEBUG
    clock_gettime (CLOCK_MONOTONIC, &Te);

    Tavg = (Tavg
            + (Te.tv_sec  - Tb.tv_sec) * 1000000000L 
            +  Te.tv_nsec - Tb.tv_nsec) / 2;

    ++iters;
#endif

    SDL_RenderPresent (r);
  }

  SDL_PauseAudioDevice (dev, 1);

#ifdef NDEBUG
  clock_getres (CLOCK_MONOTONIC, &Tp);

  printf ("res = %ld ns, avg_delta = %ld ns over iters = %lld\n",
    Tp.tv_sec * 1000000000L + Tp.tv_nsec, Tavg, iters);
#endif
end:  
  if (r)
    SDL_DestroyRenderer (r);
  
  if (win)
    SDL_DestroyWindow (win);

  if (dev)
    SDL_CloseAudioDevice (dev);

  if (ctx)
  {
    if (ctx->plan)
      fftwf_destroy_plan (ctx->plan);

    if (ctx->dft)
      fftwf_free (ctx->dft);

    if (ctx->window2)
      fftwf_free (ctx->window2);

    free(ctx);
  }
  
  SDL_Quit ();
  return 0;
}