/*
	Spectroscopic Toolkit version 1.96 by Pieter Suurmond, november 8, 2011.

    Copyright (c) 2000-2011 - Pieter Suurmond

    Permission is hereby granted, free of charge, to any person obtaining
    a copy of this software and associated documentation files
    (the "Software"), to deal in the Software without restriction,
    including without limitation the rights to use, copy, modify, merge,
    publish, distribute, sublicense, and/or sell copies of the Software,
    and to permit persons to whom the Software is furnished to do so,
    subject to the following conditions:

    The above copyright notice and this permission notice shall be
    included in all copies or substantial portions of the Software.

    Any person wishing to distribute modifications to the Software is
    requested to send the modifications to the original developer so that
    they can be incorporated into the canonical version.

    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
    MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
    IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR
    ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
    CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
    WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/

#include <stdio.h>          /* Standard headers needed here. */
#include <stdlib.h>
#include <math.h>

#include "ST_integers.h"    /* Integer datatypes (prior to ST_aiff.h). */
#include "ST_aiff.h"        /* To write result to audio file. */
#include "ST_wavelets.h"
#include "ST_engine.h"      /* Private definition of ENGINE struct. */

/*----------------------------------------------------------------------------*/
short addWAVELET(ENGINEp        E,
                 FILE*          MSG,
                 long           RGB,
                 long double    FRQ,
                 double         AMP,
                 double         PAN, 
                 long double    TOP,
                 long double    DUR,
                 short          ENV,
                 short          PHA)
{
    WAVELETp    fresh, ptr, *handle;
    UNSIGNED32  frame, frame2;
    UNSIGNED64  phase;
    double      c1, updown;
    long double realDUR, c2;

    if (E == NULL)
        return 1;
    if (FRQ >= 0.5*(long double)E->samplerate) /* Must be below Nyquist. */
        {
        if (MSG != NULL) fprintf(MSG, "FRQ (%.6Lf Hz) not below Nyquist (%.6f Hz)!\n",
                                      FRQ, 0.5*(double)E->samplerate);
        return 2;
        }
    if (FRQ < 16.0) /* Not below 16 Hz and not negative! */
        {           /* (There are DSP-reasons why FREQ may not be negative). */
        if (MSG != NULL) fprintf(MSG, "FRQ too low (%.6Lf Hz)!\n", FRQ); return 3;
        }
    c1 = fabs(AMP);             /* Contrary to FRQ, AMP may be negative,       */
    if (c1 > 1.0)               /* so if you DO want to invert phase, use AMP. */
        {
        if (MSG != NULL) fprintf(MSG, "AMP too high (%.6f)!\n", AMP); return 4;
        }
    if (c1 < (double)(1.0 / 8388608.0)) /* Store at least 23 bits of magnitude. */
        {
        if (MSG != NULL) fprintf(MSG, "AMP too low (%.6f)!\n", AMP); return 5;
        }
    if (fabs(PAN) > 1.0)                            /* Between -1 and +1. */
        {
        if (MSG != NULL) fprintf(MSG, "PAN out of range (%.4f)!\n", PAN); return 6;
        }      
    if (DUR < (4.0 / (long double)E->samplerate))        /* Minimal 4 samples. */
        {
        if (MSG != NULL) fprintf(MSG, "DUR too low (%.6Lf ms)\n", (1000.0 * DUR)); return 7;
        }
    if (DUR > 3600.0)                               /* Max 1 hour for convenience.   */
        {                                           /* So ACTUAL max. dur=5 hours.   */
        if (MSG != NULL) fprintf(MSG, "DUR too high (%.4Lf s)\n", DUR); return 8;
        }
    if ((ENV < 0) || (ENV >= kNumEnvelopeShapes))   /* 0,1 or 2. */
        {
        if (MSG != NULL) fprintf(MSG, "ENV unknown (%d)!\n", ENV); return 9;
        }
    if ((PHA < -1) || (PHA > 1))                    /* -1, 0 or 1. */
        {
        if (MSG != NULL) fprintf(MSG, "PHA out of range (%d)!\n", PHA); return 10;
        }
    c2 = TOP;                                /* Check starttime. */
    if (ENV == kENV_Gaussian)
        {
        realDUR = 5.0 * DUR;
        c2 -= (2.5 * DUR);
        }
    else /* THUS exponential: */
        {
        if (PHA) /* Illegal ENV/PHA combination? */
            {
            if (MSG != NULL) fprintf(MSG, "Exponential wavelets must have PHA=0!\n"); return 11;
            }
        realDUR = 15.0 * DUR;
        if (ENV == kENV_ExpRise)
            c2 -= realDUR;
        }
    if (c2 < 0.0)                                /* Only positive starttimes. */
        {
        if (MSG != NULL) fprintf(MSG, "TOP too early (starttime = %.4Lf seconds)!\n", c2);
        return 12;
        }
    if (c2 > 86400.0)                            /* Max 24 hours. */
        {
        if (MSG != NULL) fprintf(MSG, "TOP too late (starttime = %.2Lf hours)!\n", c2 / 3600.0);
        return 13;
        }
    if ((RGB < 0) || (RGB > 0xFFFFFF))
        {
        if (MSG != NULL) fprintf(MSG, "Illegal colour-code: %lx!\n", RGB);
        return 14;
        }
    if (!(fresh = malloc(sizeof(struct WAVELET))))
        {
        if (MSG != NULL) fprintf(MSG, "Not enough memory to allocate wavelet!\n"); return 14;
        }
    /*----------------------------------------------------- OK, so far no errors: --------*/
    if (MSG != NULL)
        fprintf(MSG, "Allocating wavelet FRQ=%12.6Lf AMP=%7.4f PAN=%4.1f TOP=%.4Lf DUR=%.4Lf COL=%lx.\n",
                FRQ, AMP, PAN, TOP, DUR, RGB); /* Notice real duration is 5 or 15 times as long! */
    /*---------------------------------------------------------------- Score-colouring: */
    fresh->rgb = RGB;
    /*---------------------------------------------------------------- Envelope timing: */
    fresh->envelopeShape = (SIGNED16)ENV;
    fresh->startFrame    = (UNSIGNED32)(0.5 + ((long double)E->samplerate * c2));
    fresh->numFrames     = (UNSIGNED32)(0.5 + ((long double)E->samplerate * realDUR));
    fresh->envelopeStep  = (UNSIGNED64)(E->two_raise_64 / (long double)fresh->numFrames); /* DOWN! */
    /*
        For the rising exponential, we must be sure that the end of the envelope table
        is exactly reached. Futhermore, the audio-oscillator needs to be aligned.
        Contrary to the audio table, we do not want the envelope table to wrap while
        reaching its' end. Sometimes (when ENV == kENV_ExpRise) we want the
        envelope phase to reach exactly the last sample.
    */
    phase = (UNSIGNED64)(fresh->numFrames - 1) * fresh->envelopeStep;
    if (phase <= (UNSIGNED64)0x8000000000000000ULL) /* CARE: not meant as -max! */
        { if (MSG != NULL) fprintf(MSG, "WARNING: envelope table wrap!\n"); }
    if (ENV == kENV_ExpRise)                  /* No compensation yet for gaussians! */
        { /* Max compensation for reverse envelopes. CARE: this is not meant as -1! */
        fresh->envelopeTime = (UNSIGNED64)0xFFFFFFFFFFFFFFFFULL - phase;
        frame  = fresh->startFrame + fresh->numFrames - 1L; /* Last frame actually calculated and written. */
        frame2 = (long)(0.5+(TOP * (long double)E->samplerate)); /* Specified frame (rounded to nearest frame). */
        if (frame < frame2)
            {
            fresh->startFrame++;    /* Let TOP become equal last calculated frame. */
            frame++;                /* Just for the next check. */
            /* Print starttime correction for debugging:
               if (MSG) fprintf(MSG, "startFrame++ compensation.\n"); */
            }
        if (frame != frame2)
            {
            if (MSG != NULL)
                fprintf(MSG, "WARNING: specified TOP = %ld but last frame = %ld!\n", frame2, frame);
            }
        /* Print envelope phase compensations for debugging:
           if (MSG != NULL) fprintf(MSG, "->envelopeTime INT = %lld\n",   fresh->envelopeTime >> kEnvTableFractBits);
           if (MSG != NULL) fprintf(MSG, "->envelopeTime FRC = %lld\n\n", fresh->envelopeTime & 0x00003FFFFFFFFFFFULL);
        */
        }
    else                                /* All other envelopes start at 0.0. */
        fresh->envelopeTime = 0ULL;     /* (not 0.5 ! to prevent wrap)       */
    /*--------------------------------------------------------------- Audio frequency:  */
    fresh->oscillatorStep = (UNSIGNED64)(0.5 + (FRQ * E->oneOverSamplerate64));
    fresh->oscillatorTime = kSineTableHalfOffset; /* Always positive, already add half. */
    if (ENV == kENV_Gaussian) /* Reach zero phase sine or cosine at TOP (middle of the envelope). */
        {
        if (PHA == -1) /* All oscills in sync at time 0. 64-bit multiplication may safely overflow. */
            fresh->oscillatorTime += (UNSIGNED64)fresh->startFrame * fresh->oscillatorStep;
        else
            {                        /* May need refinement (even or odd numFrames). */
            fresh->oscillatorTime -= ((UNSIGNED64)(fresh->numFrames >> 1) * fresh->oscillatorStep);
            if (PHA == 1) /* Cosine instead of sine. Advance one fourth of the table size. */
                fresh->oscillatorTime += ((UNSIGNED64)kSineTableSize) << (kSineTableFractBits - 2);
            }             /* CARE: we assume kSineTableSize is a muliple of 4 here! */
        }
    else if (ENV == kENV_ExpRise) /* Then PHA may never be set. Last phase actually used. */
        fresh->oscillatorTime -= ((UNSIGNED64)(fresh->numFrames - 1UL) * fresh->oscillatorStep);
    /* Else (ENV == kENV_ExpDecay) and oscillatorTime is already ok. */
    /*---------------------------------------------------------------- Panning -1 to +1: */
    AMP *= 8388608.0;  /* 2^23 = 00000000 10000000 00000000 00000000   24 bits + sign.   */
    if (AMP >= 0.0)    /* Round correctly . */
        updown = +0.5;
    else
        updown = -0.5;
    fresh->amplitudeLeft  = (SIGNED32)(updown + (AMP * sqrt(0.5 * (1.0 - PAN))));
    fresh->amplitudeRight = (SIGNED32)(updown + (AMP * sqrt(0.5 * (1.0 + PAN))));
    /*-------------------------------------------------------------------- Keep up stats: */
    E->WSTAT_number++;  /* Don't depend on this number... (yet,.. we've never needed it). */
    frame = fresh->startFrame + fresh->numFrames;
    if (frame > E->WSTAT_lastFrame)                         /* (init @ 0).   */
        E->WSTAT_lastFrame = frame;                         /* (downrounded) */
    if (fresh->oscillatorStep < E->WSTAT_loOscillStep)      /* (init @ 0xFFFFFFFFFFFFFFFF).*/
        E->WSTAT_loOscillStep = fresh->oscillatorStep;
    if (fresh->oscillatorStep > E->WSTAT_hiOscillStep)      /* (init @ 0).   */
        E->WSTAT_hiOscillStep = fresh->oscillatorStep;
    /*------------------------------------------------------------ Linked list management: */
    /* Maintain ordering on starttime so we don't need to sort the linked list afterwards. */
    if (E->WSTAT_number == 1UL)
        {
        fresh->nxt = NULL;
        E->lastWAVELET = E->firstWAVELET = fresh;
        }
    else if (E->lastWAVELET->startFrame <= fresh->startFrame)
        {
        fresh->nxt = NULL;              /* Append at end. */
        E->lastWAVELET->nxt = fresh;
        E->lastWAVELET = fresh;
        }
    else                                /* Insert. */
        {
        handle = &E->firstWAVELET;
        while (((ptr = *handle)) && (ptr->startFrame < fresh->startFrame))
            handle = &ptr->nxt;
        *handle    = fresh;
        fresh->nxt = ptr;               /* Successive wavelet or NULL. */
        }
    return 0;                           /* return 0 for ok.  */
}

/*---------------------------------------------------------------------------------------*/
short playWAVELETS(ENGINEp E, const char* filename, short bits, FILE* MSG, FILE* releaseMSG)
{
    const short numChannels = 2;  /* Let's keep this fixed and implicit for now: stereo. */
    short       error, error2;
    AUDIOFILEp  A;
    SIGNED64*   interleavAccu;  /* Dynamically allocated memory to mix and buffer audio. */
    UNSIGNED32  bytes, samplesPerBlock;
    UNSIGNED32  framesWritten = 0L;                   /* But we also count locally here. */
                                                       
    error = openAIFFWRITE(&A, filename, E->samplerate, numChannels, bits);
    if (error)
        { if (releaseMSG != NULL) fprintf(releaseMSG, "openAIFFWRITE()=%d!\n", error); return 1; }
    if (releaseMSG != NULL)
        {
        fprintf(releaseMSG, "Rounding audio to %d bits, writing blocks of %ld frames to '%s'.\n",
                            (int)bits, (long)E->framesPerBlock, filename);
        fflush(releaseMSG); /* Handy when filestream is not line-buffered (when not 'stdout'). */
        }
    /*------------------------------- ALLOCATE 64 BIT INTERLEAVED ACCUMULATOR BUFFER: --*/
    samplesPerBlock = E->framesPerBlock * numChannels;
    bytes = samplesPerBlock * sizeof(SIGNED64);
    if (!(interleavAccu = malloc(bytes)))
        {
        if (MSG != NULL)
            fprintf(MSG, "Insufficient memory (requested %ld bytes)!\n", (long)bytes);
        error = 2; goto endFL;          /* Close audio file and return error code 2. */
        }
    while (E->firstWAVELET != NULL)     /* As long as the linked list contains wavelets. */
        {
        WAVELETp    begin, *handle;
        UNSIGNED32  n;
        SIGNED64*   interleavPtr;

        /* First silence the 64-bit (stereo interleaved) output buffer. */
        interleavPtr = interleavAccu;
        n = samplesPerBlock;
        while (n--)
            *interleavPtr++ = 0LL;
        /* Then add up all wavelets within this section of time (a whole buffer).  */
        handle = &E->firstWAVELET;          /* Start at the beginning of the list. */
        while ((begin = *handle) != NULL)
            {
            UNSIGNED64  oTime, eTime;
            SIGNED16*   e;
            
            if (begin->startFrame >= (framesWritten + E->framesPerBlock))
                break;                      /* Starttimes in ascending order, we're    */
                                            /* done now with this section of time.     */           
            n = E->framesPerBlock;          /* Default is the whole accumulator.       */
            interleavPtr = interleavAccu;   /* Wavelet not already started in previous */
            if (begin->startFrame > framesWritten)
                {                           /* block (>= should also work correctly).  */
                UNSIGNED32 skip = begin->startFrame     /* Amount of frames preceding  */
                                - framesWritten;        /* start. Twice for stereo.    */
                interleavPtr += (skip * numChannels);
                n -= skip;
                } /* addWAVELET() already initialised ->oscillatorTime to be in sync.  */
            oTime = begin->oscillatorTime;      /* Load contents in cpu (before loop). */
            eTime = begin->envelopeTime;        /* oTime & eTime maybe register-vars?  */
            e = E->envTables[begin->envelopeShape];  /* Never out of range.            */
            /*------------------------------ Add this wavelet to accu-buffer: ---------*/
            while (n--)                   /* Initially, n is never 0, however, would a */
                {                         /* do-while with pre-dec be more efficient?  */       
                SIGNED32    a;            /* Enveloped but unpanned signal. */
                
                if (!begin->numFrames--)  /* Might be initially 0, is however limited. */
                    {      /* Cut-out this wavelet, we're completely done with it now. */
                    if (releaseMSG != NULL) fprintf(releaseMSG, "\
Release wavelet: start=%.3Lfs dur=%.3Lfs f=%.2LfHz.\n",
                                (long double)begin->startFrame / (long double)E->samplerate,
                                E->oneOverSamplerate64 / (long double)begin->envelopeStep,
                                (long double)begin->oscillatorStep / E->oneOverSamplerate64);
                                            /* Link previous to next (or NULL!). */
                    *handle = begin->nxt;   /* Thus skipping this one. Immediately re- */
                    free(begin);            /* lease memory so list destroys 'itself'. */
                                            /* Behave as if we've just finished the    */
                    goto trrrr;             /* preceding one, but do NOT try to write  */
                    }                       /* times back into memory (they're gone!). */
                /*-------------------------------- Calculate enveloped mono signal: ---*/
                /* 17*16 bit multiplication. Upper 18 bits of oTime and eTime are
                   used to index tables. No interpolation necessary at table-lookups.
                   Unsigned indices so we do not need to mask-off sign bits 
                   (no sign-extension). */
                a = E->sineTable[oTime >> kSineTableFractBits]
                    * (SIGNED32)e[eTime >> kEnvTableFractBits];
                /*----------------------------- Enforce 64 * 64 bit multiplication ----*/
                /*                              by SIGNED64-cast.       16*15*23 gives */
                (*interleavPtr++) -= (SIGNED64)a * begin->amplitudeLeft;   /* 54 bits. */
                (*interleavPtr++) -= (SIGNED64)a * begin->amplitudeRight;
                oTime += begin->oscillatorStep;   /* Advance fractional table-pointers */
                eTime += begin->envelopeStep;     /* (that are hopefully still in CPU  */
                }                                 /* registers), 64-bit unsigned add.  */
            begin->oscillatorTime = oTime; /* Move both audio- and envelope-times back */
            begin->envelopeTime   = eTime; /* in memory again for the next block.      */
            handle = &begin->nxt;          /* Next wavelet in the list. */
trrrr:      n = 13; /* Shameful dummy instruction! */
            }
        /* After calculating a whole buffer, write it to file. Function dumpAIFFWRITE()
           will do the scaling and rounding. It expects the upper 10 bits in the buffer 
           to be used for sign-extension only. Effectively, only 54 bits may be used. */
        error = dumpAIFFWRITE(A, interleavAccu, E->framesPerBlock);
        framesWritten += E->framesPerBlock;  /* Number of frames added to outputfile, */
        if (error)                           /* but only if no write-errors ocurred.  */
            {
            if (MSG != NULL)
                {
                fprintf(MSG, "dumpAIFFWRITE(%s)=%d", filename, error);
                if      (error == 1) fprintf(MSG, " (over- or underflow)");
                else if (error == 2) fprintf(MSG, " (i/o-error)");
                else if (error == 3) fprintf(MSG, " (file becoming huge)");
                fprintf(MSG, "!\n");
                }
            error = 3; break;               /* Immediately quit. */
            }
        }
    free(interleavAccu);                    /* Release memory allocated above. */
    printAIFFWRITE (A, MSG);                /* Both args may be NULL.          */
    if (releaseMSG != MSG)                  /* Not twice to the same stream.   */
        printAIFFWRITE (A, releaseMSG);
endFL:
    error2 = closeAIFFWRITE(&A);
    if (error2)
        {
        if (MSG != NULL) fprintf(MSG, "closeAIFFWRITE(%s)=%d!\n", filename, error2);
        if (!error) error = 4;
        }
    if (error)      /* More descriptive message already reported to MSG. */
        {
        if (releaseMSG != NULL)
            fprintf(releaseMSG, "playWAVELETS(%s)=%d (see logfile)!\n",
                                filename, error);
        }
    return error;
}

/*----------------------------------------------------------------------------*/
void releaseWAVELETS(WAVELETp* handle)   /* Free all wavelets in linked list. */
{  
    WAVELETp  nxt, ptr;

    if (handle == NULL)
        return;         /* Added in version 1.94, april 2011. */
    ptr = *handle;
    while (ptr != NULL)
        {
        nxt = ptr->nxt;
        free(ptr);
        ptr = nxt;
        }
    *handle = NULL;
}

/*----------------------------------------------------------------------------*/
static SIGNED32* newSineTable(FILE* MSG)      /* Only used within this file.  */
{
    SIGNED32        *C;

    C = malloc(kSineTableSize * sizeof(SIGNED32));
    if (C != NULL)
        {
        long double timestep, twoPi;
        UNSIGNED32  t;                        /* 32 bit to index tables.      */
        SIGNED32    *sinoid;
        
        sinoid = C;
        twoPi = 8.0 * atan(1.0);
        timestep = twoPi / (long double)kSineTableSize;
        for (t = 0; t < kSineTableSize; t++)
            {                                 /* 17 signed bits are needed. */
            long double wave = 32768.0 * sinl((long double)t * timestep);
            if (wave >= 0.0)
                wave += 0.5;
            else
                wave -= 0.5;
            *sinoid++ = (SIGNED32)wave;                 /* 32 bits needed. */
            }
        if (MSG != NULL)
            fprintf(MSG, "Created sine-table (audio 32 bit, %ld samples)\n",
                         (long)kSineTableSize);
        }
    else if (MSG != NULL)
        fprintf(MSG, "Not enough memory for sine-table (audio)!\n");
    return C;
}

/*----------------------------------------------------------------------------*/
static SIGNED16* newGaussTable(FILE* MSG)      /* Only used within this file. */
{
    SIGNED16        *G, env, *content, *contentMirror;      /* malloc inits G */
    double          timestep, gt;
    UNSIGNED32      t;                             /* 32 bit to index tables. */
    /*------------------------------------ GUASSIAN ENVELOPE -----------------*/
    /*  y=e^-(x^2)   start at x = -sqrt(25/2) --> y = 3.7267*10^-6
                       end at x = +sqrt(25/2) (about 0.12 at 15 bit unsigned.)
        Inflection points at: x = +sqrt(1/2)  equals 0.4 of the total tablesize.
                              x = -sqrt(1/2)         0.6 tablesize.
    Amplitude at TOP (where x=0) is exactly 1.0. At the inflection points, 
    amplitude is e^-0.5 which is approximately 0.6065. DUR is the effective
    duration, which is the time between the inflection points. The actually
    calculated duration is 5 times as long. One of the -6 dB point lies at:
            0.5 = e ^ (-xh^2) 
            ln(0.5) = -xh^2 
            -ln(0.5) = xh^2  
            xh = sqrt(-ln(0.5)) = 0.832554611 
    */
    G = malloc(kEnvTableSize * sizeof(SIGNED16));
    if (G != NULL)
        {
        timestep = sqrt(50.0) / (double)kEnvTableSize;
        content = G + (kEnvTableSize >> 1);               /* Start of upper half. */
        contentMirror = content;    
        for (t=0; t<(kEnvTableSize >> 1); t++)
            {
            gt = (double)t * timestep;
            gt = -1.0 * gt * gt;
            env = (SIGNED16)((-32768.0 * exp(gt)) - 0.5); /* Stored negative!!! */
            *content++ = env;                             /* Symmetrical around middle  */
            *contentMirror-- = env;                       /* Index 0 IS used so there   */
            }                                             /* will be time-err of 1 samp */
        *contentMirror = 0;
        if (MSG != NULL)
            fprintf(MSG, "Created symmetrical gauss-table (subaudio 16 bit, %ld samples)\n",
                         (long)kEnvTableSize);
        }
    else
        if (MSG != NULL) fprintf(MSG, "Not enough memory for gauss-table (envelope)\n");
    return(G);
}

/*----------------------------------------------------------------------------------------*/
static SIGNED16* newExpTable(short rise, FILE* MSG) /* Only used within this file. Decay- */
{                                                   /* ing if argument rise is zero,      */
    SIGNED16        *G, *content;                   /* rising if rise is nonzero.         */
    double          timestep;
    UNSIGNED32      t, tt;                        /* 32 bits to index tables. */
    const char*     riseTxt[2] = { "rising", "decaying" };
    /*-------------------------------------- EXPONENTIAL DECAY ---------------*/
    /*  y=e^-x      start at x = 0       --> y = 1
                      end at x = ln(0.5) --> y = 0.5
    Amplitude at TOP (where x=0) is exactly 1.0.
    We call DUR the effective duration, which is the time between the top (1.0) and
    the -6 dB point (0.5). The actually calculated duration is 15 times as long.
    */
    G = malloc(kEnvTableSize * sizeof(SIGNED16));
    if (G != NULL)
        {
        content = G;
        timestep = 15.0 * log(0.5) / (double)kEnvTableSize;   /* 'log()' is natural log! */
        for (t=0; t<kEnvTableSize; t++)
            {                                                      /* Stored negative!!! */
            if (rise) tt = kEnvTableSize - 1 - t; else tt = t;
            *content++ = (SIGNED16)((-32768.0 * exp((double)tt * timestep)) - 0.5);
            }
        if (MSG != NULL)
            fprintf(MSG, "Created %s exp-table (subaudio 16 bit, %ld samples)\n",
                         riseTxt[rise], (long)kEnvTableSize);
        }
    else
        if (MSG != NULL) fprintf(MSG, "Not enough memory for exp-table (envelope)\n");
    return G;
}

/*----------------------------------------------------------------------------*/
ENGINEp newENGINE(long sr, long blocksize, FILE* msg)
{
    ENGINEp E = NULL;

    if ((sr < 64L) || (sr > 144000L))                         /* Be tolerant. */
        { if (msg != NULL) fprintf(msg, "Invalid samplerate (%ld Hz) for wavelet-engine!\n", sr); }
    else if ((blocksize < 1L) || (blocksize > 1048576L)) /* Not more than 1 million frames.*/
        { if (msg != NULL) fprintf(msg, "Invalid blocksize (%ld sampleframes)!\n", blocksize); }
    /*------------------------------------ First a small datasize-check: -----------------*/
    else if ((sizeof(SIGNED16)<2)   || (sizeof(SIGNED32)<4)   || (sizeof(SIGNED64)<8)   ||
             (sizeof(UNSIGNED16)<2) || (sizeof(UNSIGNED32)<4) || (sizeof(UNSIGNED64)<8) ||
             (sizeof(float)<4)      || (sizeof(double)<8)     || (sizeof(char)!=1)      ||
             (sizeof(long double)<12) || 
             E) /* Surpress 'constant expression'-warning with E. */
        { if (msg) fprintf(msg, "Datatype(s) incompatible with platform/compiler!\n"); }
        /* Which should better be tested @ compile time instead of runtime. */
    /*------------------------------------ Allocate & init DSP-engine itself: ------*/
    else
        {
        if (msg != NULL)
            fprintf(msg, "Selected audio samplerate = %ld Hertz.\n", sr);
        E = malloc(sizeof(struct ENGINE));
        if (E != NULL)
            {
            E->two_raise_64 = powl(2.0, 64.0);  /* Constant. */
            
            E->samplerate = sr;                         /* Used by addWAVELET(), play-  */
            E->oneOverSamplerate64 = E->two_raise_64    /* WAVELETS(), scoreWAVELETS(). */
                                   / (long double)sr;
            E->framesPerBlock = blocksize;
            
            E->lastWAVELET = E->firstWAVELET = (WAVELETp)NULL; /* Start with no wavelets. */
            
            E->WSTAT_number = 0L;                   /* Incremented by addWAVELET(). */
            E->WSTAT_lastFrame = 0L;                /* Updated by addWAVELET(). */
            E->WSTAT_loOscillStep = (UNSIGNED64)0xFFFFFFFFFFFFFFFFULL; /* NOT -1! */
            E->WSTAT_hiOscillStep = (UNSIGNED64)0x0000000000000000ULL;
            /*------------------------------------- ALLOCATE AND FILL WAVETABLES: */
            E->sineTable                = newSineTable(msg);
            E->envTables[kENV_Gaussian] = newGaussTable(msg);
            E->envTables[kENV_ExpDecay] = newExpTable(0, msg);
            E->envTables[kENV_ExpRise]  = newExpTable(1, msg);
            if ((E->sineTable == NULL) ||                       /* All 4 must be there. */
                (E->envTables[kENV_Gaussian] == NULL) ||
                (E->envTables[kENV_ExpDecay] == NULL) ||
                (E->envTables[kENV_ExpRise] == NULL))
                {
                releaseENGINE(&E);                      /* Also sets E to NULL. */
                if (msg != NULL)
                    fprintf(msg, "Not enough memory to allocate wavetables!\n");
                }
            }
        else if (msg != NULL)
            fprintf(msg, "Not enough memory to allocate a new wavelet engine!\n");
        }
    return E;
}

/*----------------------------------------------------------------------------*/
long getSamplerateENGINE(ENGINEp E)
{
    return E->samplerate;
}

/*----------------------------------------------------------------------------*/
void releaseENGINE(ENGINEp* E)
{
    short    i;

    if ((E != NULL) && (*E != NULL))
        {
        for (i=0; i<kNumEnvelopeShapes; i++)    /* Release all envelope   */
            if ((*E)->envTables[i] != NULL)     /* tables that are there. */
                free((*E)->envTables[i]);
        if ((*E)->sineTable != NULL)
            free((*E)->sineTable);
        releaseWAVELETS(&((*E)->firstWAVELET)); /* Cleanup linked list.   */
        free(*E);                               /* Release ENGINE itself. */
        *E = (ENGINEp)NULL;                     /* Clear reference.       */
        }
}