/*  synaescope.cpp
 *  Copyright (C) 1999,2002 Richard Boulton <richard@tartarus.org>
 *
 *  Much code copied from Synaesthesia - a program to display sound
 *  graphically, by Paul Francis Harrison <pfh@yoyo.cc.monash.edu.au>
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
 */

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include "synaescope.h"

#include <pthread.h>
#include <dirent.h>
#include <sys/stat.h>
#include <sys/time.h>
#include <time.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <string.h>
#include <assert.h>

#define SCOPE_BG_RED 0
#define SCOPE_BG_GREEN 0
#define SCOPE_BG_BLUE 0

#define FFT_BUFFER_SIZE_LOG 9
#define FFT_BUFFER_SIZE (1 << FFT_BUFFER_SIZE_LOG)

#define syn_width 320
#define syn_height 200
#define brightMin 200
#define brightMax 2000
#define brightDec 10
#define brightInc 6
#define brTotTargetLow 5000
#define brTotTargetHigh 15000

static int autobrightness = 1;  /* Whether to use automatic brightness adjust */
static unsigned int brightFactor = 400;
static unsigned char output[syn_width * syn_height * 2];
static guint32 display[syn_width * syn_height];
static gint16 pcmt_l[FFT_BUFFER_SIZE];
static gint16 pcmt_r[FFT_BUFFER_SIZE];
static gint16 pcm_l[FFT_BUFFER_SIZE];
static gint16 pcm_r[FFT_BUFFER_SIZE];
static double fftout_l[FFT_BUFFER_SIZE];
static double fftout_r[FFT_BUFFER_SIZE];
static double fftmult[FFT_BUFFER_SIZE / 2 + 1];
static double corr_l[FFT_BUFFER_SIZE];
static double corr_r[FFT_BUFFER_SIZE];
static int clarity[FFT_BUFFER_SIZE];    /* Surround sound */
static double cosTable[FFT_BUFFER_SIZE];
static double negSinTable[FFT_BUFFER_SIZE];
static int bitReverse[FFT_BUFFER_SIZE];
static int scaleDown[256];

static void synaes_fft (double *x, double *y);
static void synaescope_coreGo (void);

#define SYNAESCOPE_DOLOOP() \
while (running) { \
    gint bar; \
    guint val; \
    gint val2; \
    unsigned char *outptr = output; \
	int w; \
\
    synaescope_coreGo(); \
\
    outptr = output; \
    for (w=0; w < syn_width * syn_height; w++) { \
	bits[w] = colEq[(outptr[0] >> 4) + (outptr[1] & 0xf0)]; \
	outptr += 2; \
    } \
\
    GDK_THREADS_ENTER(); \
    gdk_draw_image(win,gc,image,0,0,0,0,-1,-1); \
    gdk_flush(); \
    GDK_THREADS_LEAVE(); \
    dosleep(SCOPE_SLEEP); \
}

static inline void
addPixel (unsigned char *output, int x, int y, int br1, int br2)
{
  unsigned char *p;

  if (x < 0 || x >= syn_width || y < 0 || y >= syn_height)
    return;

  p = output + x * 2 + y * syn_width * 2;
  if (p[0] < 255 - br1)
    p[0] += br1;
  else
    p[0] = 255;
  if (p[1] < 255 - br2)
    p[1] += br2;
  else
    p[1] = 255;
}

static inline void
addPixelFast (unsigned char *p, int br1, int br2)
{
  if (p[0] < 255 - br1)
    p[0] += br1;
  else
    p[0] = 255;
  if (p[1] < 255 - br2)
    p[1] += br2;
  else
    p[1] = 255;
}

static void
synaescope_coreGo (void)
{
  int i, j;
  register unsigned long *ptr;
  register unsigned long *end;
  int heightFactor;
  int actualHeight;
  int heightAdd;
  double brightFactor2;
  long int brtot;

  memcpy (pcm_l, pcmt_l, sizeof (pcm_l));
  memcpy (pcm_r, pcmt_r, sizeof (pcm_r));

  for (i = 0; i < FFT_BUFFER_SIZE; i++) {
    fftout_l[i] = pcm_l[i];
    fftout_r[i] = pcm_r[i];
  }

  synaes_fft (fftout_l, fftout_r);

  for (i = 0 + 1; i < FFT_BUFFER_SIZE; i++) {
    double x1 = fftout_l[bitReverse[i]];
    double y1 = fftout_r[bitReverse[i]];
    double x2 = fftout_l[bitReverse[FFT_BUFFER_SIZE - i]];
    double y2 = fftout_r[bitReverse[FFT_BUFFER_SIZE - i]];
    double aa, bb;

    corr_l[i] = sqrt (aa = (x1 + x2) * (x1 + x2) + (y1 - y2) * (y1 - y2));
    corr_r[i] = sqrt (bb = (x1 - x2) * (x1 - x2) + (y1 + y2) * (y1 + y2));
    clarity[i] = (int) (
        ((x1 + x2) * (x1 - x2) + (y1 + y2) * (y1 - y2)) / (aa + bb) * 256);
  }

  /* Asger Alstrupt's optimized 32 bit fade */
  /* (alstrup@diku.dk) */
  ptr = (unsigned long *) output;
  end = (unsigned long *) (output + syn_width * syn_height * 2);
  do {
    /*Bytewize version was: *(ptr++) -= *ptr+(*ptr>>1)>>4; */
    if (*ptr) {
      if (*ptr & 0xf0f0f0f0) {
        *ptr = *ptr - ((*ptr & 0xf0f0f0f0) >> 4) - ((*ptr & 0xe0e0e0e0) >> 5);
      } else {
        *ptr = (*ptr * 14 >> 4) & 0x0f0f0f0f;
        /*Should be 29/32 to be consistent. Who cares. This is totally */
        /* hacked anyway.  */
        /*unsigned char *subptr = (unsigned char*)(ptr++); */
        /*subptr[0] = (int)subptr[0] * 29 / 32; */
        /*subptr[1] = (int)subptr[0] * 29 / 32; */
        /*subptr[2] = (int)subptr[0] * 29 / 32; */
        /*subptr[3] = (int)subptr[0] * 29 / 32; */
      }
    }
    ptr++;
  } while (ptr < end);

  heightFactor = FFT_BUFFER_SIZE / 2 / syn_height + 1;
  actualHeight = FFT_BUFFER_SIZE / 2 / heightFactor;
  heightAdd = (syn_height + actualHeight) >> 1;

  /* Correct for window size */
  brightFactor2 = (brightFactor / 65536.0 / FFT_BUFFER_SIZE) *
      sqrt (actualHeight * syn_width / (320.0 * 200.0));

  brtot = 0;
  for (i = 1; i < FFT_BUFFER_SIZE / 2; i++) {
    /*int h = (int)( corr_r[i]*280 / (corr_l[i]+corr_r[i]+0.0001)+20 ); */
    if (corr_l[i] > 0 || corr_r[i] > 0) {
      int h = (int) (corr_r[i] * syn_width / (corr_l[i] + corr_r[i]));

/*      int h = (int)( syn_width - 1 ); */
      int br1, br2, br = (int) ((corr_l[i] + corr_r[i]) * i * brightFactor2);
      int px = h, py = heightAdd - i / heightFactor;

      brtot += br;
      br1 = br * (clarity[i] + 128) >> 8;
      br2 = br * (128 - clarity[i]) >> 8;
      if (br1 < 0)
        br1 = 0;
      else if (br1 > 255)
        br1 = 255;
      if (br2 < 0)
        br2 = 0;
      else if (br2 > 255)
        br2 = 255;
      /*unsigned char *p = output+ h*2+(164-((i<<8)>>FFT_BUFFER_SIZE_LOG))*(syn_width*2);  */

      if (px < 30 || py < 30 || px > syn_width - 30 || py > syn_height - 30) {
        addPixel (output, px, py, br1, br2);
        for (j = 1; br1 > 0 || br2 > 0;
            j++, br1 = scaleDown[br1], br2 = scaleDown[br2]) {
          addPixel (output, px + j, py, br1, br2);
          addPixel (output, px, py + j, br1, br2);
          addPixel (output, px - j, py, br1, br2);
          addPixel (output, px, py - j, br1, br2);
        }
      } else {
        unsigned char *p = output + px * 2 + py * syn_width * 2, *p1 = p, *p2 =
            p, *p3 = p, *p4 = p;
        addPixelFast (p, br1, br2);
        for (; br1 > 0 || br2 > 0; br1 = scaleDown[br1], br2 = scaleDown[br2]) {
          p1 += 2;
          addPixelFast (p1, br1, br2);
          p2 -= 2;
          addPixelFast (p2, br1, br2);
          p3 += syn_width * 2;
          addPixelFast (p3, br1, br2);
          p4 -= syn_width * 2;
          addPixelFast (p4, br1, br2);
        }
      }
    }
  }

  /* Apply autoscaling: makes quiet bits brighter, and loud bits
   * darker, but still keeps loud bits brighter than quiet bits. */
  if (brtot != 0 && autobrightness) {
    long int brTotTarget = brTotTargetHigh;

    if (brightMax != brightMin) {
      brTotTarget -= ((brTotTargetHigh - brTotTargetLow) *
          (brightFactor - brightMin)) / (brightMax - brightMin);
    }
    if (brtot < brTotTarget) {
      brightFactor += brightInc;
      if (brightFactor > brightMax)
        brightFactor = brightMax;
    } else {
      brightFactor -= brightDec;
      if (brightFactor < brightMin)
        brightFactor = brightMin;
    }
    /* printf("brtot: %ld\tbrightFactor: %d\tbrTotTarget: %d\n",
       brtot, brightFactor, brTotTarget); */
  }
}

#define BOUND(x) ((x) > 255 ? 255 : (x))
#define PEAKIFY(x) BOUND((x) - (x)*(255-(x))/255/2)

static void
synaescope32 ()
{
  unsigned char *outptr;
  guint32 colEq[256];
  int i;
  guint32 bg_color;

  for (i = 0; i < 256; i++) {
    int red = PEAKIFY ((i & 15 * 16));
    int green = PEAKIFY ((i & 15) * 16 + (i & 15 * 16) / 4);
    int blue = PEAKIFY ((i & 15) * 16);

    colEq[i] = (red << 16) + (green << 8) + blue;
  }
  bg_color = (SCOPE_BG_RED << 16) + (SCOPE_BG_GREEN << 8) + SCOPE_BG_BLUE;

  synaescope_coreGo ();

  outptr = output;
  for (i = 0; i < syn_width * syn_height; i++) {
    display[i] = colEq[(outptr[0] >> 4) + (outptr[1] & 0xf0)];
    outptr += 2;
  }
}


#if 0
static void
synaescope16 (void *data)
{
  guint16 *bits;
  guint16 colEq[256];
  int i;
  GdkWindow *win;
  GdkColormap *c;
  GdkVisual *v;
  GdkGC *gc;
  GdkColor bg_color;

  win = (GdkWindow *) data;
  GDK_THREADS_ENTER ();
  c = gdk_colormap_get_system ();
  gc = gdk_gc_new (win);
  v = gdk_window_get_visual (win);

  for (i = 0; i < 256; i++) {
    GdkColor color;

    color.red = PEAKIFY ((i & 15 * 16)) << 8;
    color.green = PEAKIFY ((i & 15) * 16 + (i & 15 * 16) / 4) << 8;
    color.blue = PEAKIFY ((i & 15) * 16) << 8;
    gdk_color_alloc (c, &color);
    colEq[i] = color.pixel;
  }

  /* Create render image */
  if (image) {
    gdk_image_destroy (image);
    image = NULL;
  }
  image = gdk_image_new (GDK_IMAGE_FASTEST, v, syn_width, syn_height);
  bg_color.red = SCOPE_BG_RED << 8;
  bg_color.green = SCOPE_BG_GREEN << 8;
  bg_color.blue = SCOPE_BG_BLUE << 8;
  gdk_color_alloc (c, &bg_color);
  GDK_THREADS_LEAVE ();

  assert (image);
  assert (image->bpp == 2);

  bits = (guint16 *) image->mem;

  running = 1;

  SYNAESCOPE_DOLOOP ();
}


static void
synaescope8 (void *data)
{
  unsigned char *outptr;
  guint8 *bits;
  guint8 colEq[256];
  int i;
  GdkWindow *win;
  GdkColormap *c;
  GdkVisual *v;
  GdkGC *gc;
  GdkColor bg_color;

  win = (GdkWindow *) data;
  GDK_THREADS_ENTER ();
  c = gdk_colormap_get_system ();
  gc = gdk_gc_new (win);
  v = gdk_window_get_visual (win);

  for (i = 0; i < 64; i++) {
    GdkColor color;

    color.red = PEAKIFY ((i & 7 * 8) * 4) << 8;
    color.green = PEAKIFY ((i & 7) * 32 + (i & 7 * 8) * 2) << 8;
    color.blue = PEAKIFY ((i & 7) * 32) << 8;
    gdk_color_alloc (c, &color);
    colEq[i * 4] = color.pixel;
    colEq[i * 4 + 1] = color.pixel;
    colEq[i * 4 + 2] = color.pixel;
    colEq[i * 4 + 3] = color.pixel;
  }

  /* Create render image */
  if (image) {
    gdk_image_destroy (image);
    image = NULL;
  }
  image = gdk_image_new (GDK_IMAGE_FASTEST, v, syn_width, syn_height);
  bg_color.red = SCOPE_BG_RED << 8;
  bg_color.green = SCOPE_BG_GREEN << 8;
  bg_color.blue = SCOPE_BG_BLUE << 8;
  gdk_color_alloc (c, &bg_color);
  GDK_THREADS_LEAVE ();

  assert (image);
  assert (image->bpp == 1);

  bits = (guint8 *) image->mem;

  running = 1;

  SYNAESCOPE_DOLOOP ();
}
#endif

#if 0
static void
run_synaescope (void *data)
{
  switch (depth) {
    case 8:
      synaescope8 (win);
      break;
    case 16:
      synaescope16 (win);
      break;
    case 24:
    case 32:
      synaescope32 (win);
      break;

  }
}

static void
start_synaescope (void *data)
{
  init_synaescope_window ();
}

#endif


static int
bitReverser (int i)
{
  int sum = 0;
  int j;

  for (j = 0; j < FFT_BUFFER_SIZE_LOG; j++) {
    sum = (i & 1) + sum * 2;
    i >>= 1;
  }

  return sum;
}

static void
init_synaescope ()
{
  int i;

  for (i = 0; i <= FFT_BUFFER_SIZE / 2 + 1; i++) {
    double mult = (double) 128 / ((FFT_BUFFER_SIZE * 16384) ^ 2);

    /* Result now guaranteed (well, almost) to be in range 0..128 */

    /* Low values represent more frequencies, and thus get more */
    /* intensity - this helps correct for that. */
    mult *= log (i + 1) / log (2);

    mult *= 3;                  /* Adhoc parameter, looks about right for me. */

    fftmult[i] = mult;
  }

  for (i = 0; i < FFT_BUFFER_SIZE; i++) {
    negSinTable[i] = -sin (M_PI * 2 / FFT_BUFFER_SIZE * i);
    cosTable[i] = cos (M_PI * 2 / FFT_BUFFER_SIZE * i);
    bitReverse[i] = bitReverser (i);
  }

  for (i = 0; i < 256; i++)
    scaleDown[i] = i * 200 >> 8;

  memset (output, 0, syn_width * syn_height * 2);
}

static void
synaes_fft (double *x, double *y)
{
  int n2 = FFT_BUFFER_SIZE;
  int n1;
  int twoToTheK;
  int j;

  for (twoToTheK = 1; twoToTheK < FFT_BUFFER_SIZE; twoToTheK *= 2) {
    n1 = n2;
    n2 /= 2;
    for (j = 0; j < n2; j++) {
      double c = cosTable[j * twoToTheK & (FFT_BUFFER_SIZE - 1)];
      double s = negSinTable[j * twoToTheK & (FFT_BUFFER_SIZE - 1)];
      int i;

      for (i = j; i < FFT_BUFFER_SIZE; i += n1) {
        int l = i + n2;
        double xt = x[i] - x[l];
        double yt = y[i] - y[l];

        x[i] = (x[i] + x[l]);
        y[i] = (y[i] + y[l]);
        x[l] = xt * c - yt * s;
        y[l] = xt * s + yt * c;
      }
    }
  }
}

static void
synaescope_set_data (gint16 data[2][512])
{
  int i;
  gint16 *newset_l = pcmt_l;
  gint16 *newset_r = pcmt_r;

  for (i = 0; i < FFT_BUFFER_SIZE; i++) {
    newset_l[i] = data[0][i];
    newset_r[i] = data[1][i];
  }
}

void
synaesthesia_init (guint32 resx, guint32 resy)
{
  init_synaescope ();
}

guint32 *
synaesthesia_update (gint16 data[2][512])
{
  synaescope_set_data (data);
  synaescope32 ();
  return display;
}

void
synaesthesia_close ()
{
}