/*
 * concert_a.c - Generate a pure 441 Hz tone for 5 seconds.
 *
 * Concert A is usually pitched somewhere between 440 Hz (USA) 
 *    and 442 Hz (Europe).
 *
 *		     Written by Scott Porter, Gints Klimanis 
 *			Silicon Graphics Computer Systems
 *				 1994
 */
#include <dmedia/audio.h>
#include <math.h>
#include <stdio.h>

/* A semitone ratio - 2^(1/12) */
#define SEMITONE 1.95500086538685 
#define MIDI2Freq(x) ((double 440)*(SEMITONE^(x-69.0)))

long GetAudioHardwareInputRate(void);
long GetAudioHardwareOutputRate(void);
void ComputeSinusoidPlus(short *out, int length, double normalFrequency,
		     double phase, double amplitude);
void ComputeSinusoid(short *out, int length, double normalFrequency,
		     double phase, double amplitude);
void GrowHarmony(double frequency, double length_sec, double samplingRate,
		 ALport audioPort);
void AdSyn(double frequency, double length_sec, double samplingRate,
		 ALport audioPort);



#define NPARTIALS 32

main()
{

  int	    i,j;
  short	    buf[60000];
  ALconfig  config;
  ALport    audioPort;
  double    samplingRate;
  int	    length;
  double    amplitude, frequency, phase, harmonic_frequency;
  
  ALseterrorhandler(0);

  /* query audio hardware for output sampling */
  /* if undefined, default to 44100 Hertz */
  samplingRate = GetAudioHardwareOutputRate();
  if (samplingRate == AL_RATE_UNDEFINED)
    samplingRate = 44100;
  length = samplingRate;
  
  config = ALnewconfig();
  ALsetqueuesize(config, samplingRate);
  ALsetchannels(config, AL_MONO);
  
  audioPort = ALopenport("outport", "w", config);
  if (!audioPort) 
  {
    fprintf(stderr, "couldn't open port\n");
    exit(1);
  }

  AdSyn(150, 1, samplingRate, audioPort); 
//  GrowHarmony(100.0, 0.25, samplingRate, audioPort);

  
  /* allow audio to drain from port */
  while (ALgetfilled(audioPort) > 0) 
    sginap(1);

  ALfreeconfig(config);
  ALcloseport(audioPort);
}



void AdSyn(double frequency, double length_sec, double samplingRate,
		 ALport audioPort)
{
  double    coef;
  int	    i,j;
  short	    buf[60000];
  ALconfig  config;
  int	    length;
  double    amplitude, phase, harmonic_frequency, harmonic_amplitude;
  double    a[NPARTIALS], f[NPARTIALS];

  phase = M_PI/2;
  amplitude = 32767/NPARTIALS;
  length = samplingRate*length_sec;
 
  ComputeSinusoid(buf, length, frequency/samplingRate, phase, amplitude);  
  for (i=1;i<NPARTIALS;i++)
  {
    coef = exp((double) -i/4.0);
    printf("\n Coeff: %f; i: %d",coef, i);
    f[i] = frequency*(i+1);
    a[i] = amplitude*coef;
    ComputeSinusoidPlus(buf, length, f[i]/samplingRate, phase, a[i]);
  }
  ALwritesamps(audioPort, buf, length);


} 

void GrowHarmony(double frequency, double length_sec, double samplingRate,
		 ALport audioPort)
{		 
  int	    i,j;
  short	    buf[60000];
  ALconfig  config;
  int	    length;
  double    amplitude, phase, harmonic_frequency;

 /* generate quarter of a second of sine wave */
  amplitude = 32767/NPARTIALS;	/* for wave in range [-32767 .. 32767] */
  phase = M_PI/2;		/* Radians */
  length = samplingRate*length_sec;
  
  /* write 5 seconds to audio hardware */
  for (i = 1; i <= NPARTIALS; i++)
  {
    ComputeSinusoid(buf, length, frequency/samplingRate, phase, amplitude);
    for (j=2, harmonic_frequency=2*frequency; j<=i; j++)
    {
      printf("\nComputing partial %d",j);
      harmonic_frequency += frequency;
      ComputeSinusoidPlus(buf, length, harmonic_frequency/samplingRate, phase, amplitude);
    }
    printf("\nplaying %d partials",i);
    ALwritesamps(audioPort, buf, length);

  }
} 


long GetAudioHardwareInputRate(void)
{
  long buffer[6];

  /* acquire state variables of audio hardware */
  buffer[0] = AL_INPUT_RATE;
  buffer[2] = AL_INPUT_SOURCE;
  buffer[4] = AL_DIGITAL_INPUT_RATE;
  ALgetparams(AL_DEFAULT_DEVICE, buffer, 6);
  
  /* for input sources microphone or line and input rate not AES word clock */
  if	((buffer[3] != AL_INPUT_DIGITAL)&&(buffer[1] > 0))
    return (buffer[1]);

  /* for input rate AES word clock and machine has ability to read digital 
     input sampling rate, return AL_DIGITAL_INPUT_RATE */
  else if (ALgetdefault(AL_DEFAULT_DEVICE, AL_DIGITAL_INPUT_RATE) >= 0) 
    return (buffer[5]);

  /* could not determine sampling rate */
  else
    return (AL_RATE_UNDEFINED);
}   /* ---- end GetAudioHardwareInputRate() ---- */



/* ******************************************************************
 * GetAudioHardwareOutputRate	acquire audio hardware output sampling rate
 * ****************************************************************** */
long GetAudioHardwareOutputRate(void)
{
long	buffer[4];

buffer[0] = AL_OUTPUT_RATE;
buffer[2] = AL_DIGITAL_INPUT_RATE;
ALgetparams(AL_DEFAULT_DEVICE, buffer, 4);

/* return output rate (Hertz) */
if	(buffer[1] > 0) 
    return (buffer[1]);

/* when output rate is input rate, get input rate */
else if (AL_RATE_INPUTRATE == buffer[1])
    return (GetAudioHardwareInputRate());

/* for input rate AES word clock and machine has ability to read digital 
    input sampling rate, return AL_DIGITAL_INPUT_RATE */
else if (ALgetdefault(AL_DEFAULT_DEVICE, AL_DIGITAL_INPUT_RATE) >= 0) 
    return (buffer[3]);

/* could not determine sampling rate */
else
    return (AL_RATE_UNDEFINED);
} /* ---- end GetAudioHardwareOutputRate() ---- */





  
/* ******************************************************************
 * ComputeSinusoidPlus:    use waveguide oscillator to compute sinusoid
 *           if there is something already there, it will add to it.
 * ****************************************************************** */
void ComputeSinusoidPlus(short *out, int length, double normalFrequency,
			 double phase, double amplitude)
/*
    out			ptr to output 
    length		length of output
    normalFrequency	frequency/samplingFrequency
    phase		in Radians.  0=cosine wave, PI/2= sine wave
    amplitude		maximum&minimum of sinusoid
*/
{
int	i;
double	sum1, sum2;
double	z1, z2;
double	amplitudeCoefficient, tuningCoefficient;

/* 
 * -------- compute transformer normalized waveguide sinusoid oscillator 
 */
/*
more efficient version of:

argInc = frequency*2*M_PI/((float) samplingRate);
arg = phase;
for (i = 0; i < length; i++, arg += argInc) 
    buf[i] = (short) (amplitude*cos(arg));
*/

/* ring in z2:  z2 output is ring-value amplitude cosine wave  
		z1 output is low amplitude -sine 

   ring in z1:  z1 output is ring-value amplitude cosine wave
		z2 output is huge amplitude sine wave.  	

   to start at arbitrary initial phase:
	z1 = -amplitudeCoefficient*amplitude*sin(phase)
	z2 = amplitude*cos(phase)
*/
tuningCoefficient = cos(2*M_PI*normalFrequency);
amplitudeCoefficient = sqrt(((1-tuningCoefficient)/(1+tuningCoefficient)));
z2 = amplitude*cos(phase);
z1 = -amplitudeCoefficient*amplitude*sin(phase);
for (i = 0; i < length; i++)
    {
    sum2 = (z1 + z2)*tuningCoefficient;
    sum1 = sum2 - z2;
    z2 = z1 + sum2;
    z1 = sum1;

    out[i] += (short) z2;
    }

}   /* ---- end ComputeSinusoidPlus() ---- */




void ComputeSinusoid(short *out, int length, double normalFrequency,
			 double phase, double amplitude)
/*
    out			ptr to output 
    length		length of output
    normalFrequency	frequency/samplingFrequency
    phase		in Radians.  0=cosine wave, PI/2= sine wave
    amplitude		maximum&minimum of sinusoid
*/
{
int	i;
double	sum1, sum2;
double	z1, z2;
double	amplitudeCoefficient, tuningCoefficient;

/* 
 * -------- compute transformer normalized waveguide sinusoid oscillator 
 */
/*
more efficient version of:

argInc = frequency*2*M_PI/((float) samplingRate);
arg = phase;
for (i = 0; i < length; i++, arg += argInc) 
    buf[i] = (short) (amplitude*cos(arg));
*/

/* ring in z2:  z2 output is ring-value amplitude cosine wave  
		z1 output is low amplitude -sine 

   ring in z1:  z1 output is ring-value amplitude cosine wave
		z2 output is huge amplitude sine wave.  	

   to start at arbitrary initial phase:
	z1 = -amplitudeCoefficient*amplitude*sin(phase)
	z2 = amplitude*cos(phase)
*/
tuningCoefficient = cos(2*M_PI*normalFrequency);
amplitudeCoefficient = sqrt(((1-tuningCoefficient)/(1+tuningCoefficient)));
z2 = amplitude*cos(phase);
z1 = -amplitudeCoefficient*amplitude*sin(phase);
for (i = 0; i < length; i++)
    {
    sum2 = (z1 + z2)*tuningCoefficient;
    sum1 = sum2 - z2;
    z2 = z1 + sum2;
    z1 = sum1;

    out[i] = (short) z2;
    }

}   /* ---- end ComputeSinusoid() ---- */