Logo Search packages:      
Sourcecode: csound version File versions

reverbsc.c

/*
    reverbsc.c:

    8 delay line FDN reverb, with feedback matrix based upon
    physical modeling scattering junction of 8 lossless waveguides
    of equal characteristic impedance. Based on Julius O. Smith III,
    "A New Approach to Digital Reverberation using Closed Waveguide
    Networks," Proceedings of the International Computer Music
    Conference 1985, p. 47-53 (also available as a seperate
    publication from CCRMA), as well as some more recent papers by
    Smith and others.

    Csound orchestra version coded by Sean Costello, October 1999

    C implementation (C) 2005 Istvan Varga

    This file is part of Csound.

    The Csound Library is free software; you can redistribute it
    and/or modify it under the terms of the GNU Lesser General Public
    License as published by the Free Software Foundation; either
    version 2.1 of the License, or (at your option) any later version.

    Csound 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 Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public
    License along with Csound; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
    02111-1307 USA
*/

#include "csdl.h"
#include <math.h>

#define DEFAULT_SRATE   44100.0
#define MIN_SRATE       5000.0
#define MAX_SRATE       1000000.0
#define MAX_PITCHMOD    20.0
#define DELAYPOS_SHIFT  28
#define DELAYPOS_SCALE  0x10000000
#define DELAYPOS_MASK   0x0FFFFFFF

/* reverbParams[n][0] = delay time (in seconds)                     */
/* reverbParams[n][1] = random variation in delay time (in seconds) */
/* reverbParams[n][2] = random variation frequency (in 1/sec)       */
/* reverbParams[n][3] = random seed (0 - 32767)                     */

static const double reverbParams[8][4] = {
    { (2473.0 / DEFAULT_SRATE), 0.0010, 3.100,  1966.0 },
    { (2767.0 / DEFAULT_SRATE), 0.0011, 3.500, 29491.0 },
    { (3217.0 / DEFAULT_SRATE), 0.0017, 1.110, 22937.0 },
    { (3557.0 / DEFAULT_SRATE), 0.0006, 3.973,  9830.0 },
    { (3907.0 / DEFAULT_SRATE), 0.0010, 2.341, 20643.0 },
    { (4127.0 / DEFAULT_SRATE), 0.0011, 1.897, 22937.0 },
    { (2143.0 / DEFAULT_SRATE), 0.0017, 0.891, 29491.0 },
    { (1933.0 / DEFAULT_SRATE), 0.0006, 3.221, 14417.0 }
};

static const double outputGain  = 0.35;
static const double jpScale     = 0.25;

00065 typedef struct {
    int         writePos;
    int         bufferSize;
    int         readPos;
    int         readPosFrac;
    int         readPosFrac_inc;
    int         dummy;
    int         seedVal;
    int         randLine_cnt;
    double      filterState;
    MYFLT       buf[1];
} delayLine;

00078 typedef struct {
    OPDS        h;
    MYFLT       *aoutL, *aoutR, *ainL, *ainR, *kFeedBack, *kLPFreq;
    MYFLT       *iSampleRate, *iPitchMod, *iSkipInit;
    double      sampleRate;
    double      dampFact;
    MYFLT       prv_LPFreq;
    int         initDone;
    delayLine   *delayLines[8];
    AUXCH       auxData;
} SC_REVERB;

static int delay_line_max_samples(SC_REVERB *p, int n)
{
    double  maxDel;

    maxDel = reverbParams[n][0];
    maxDel += (reverbParams[n][1] * (double) *(p->iPitchMod) * 1.125);
    return (int) (maxDel * p->sampleRate + 16.5);
}

static int delay_line_bytes_alloc(SC_REVERB *p, int n)
{
    int nBytes;

    nBytes = (int) sizeof(delayLine) - (int) sizeof(MYFLT);
    nBytes += (delay_line_max_samples(p, n) * (int) sizeof(MYFLT));
    nBytes = (nBytes + 15) & (~15);
    return nBytes;
}

static void next_random_lineseg(SC_REVERB *p, delayLine *lp, int n)
{
    double  prvDel, nxtDel, phs_incVal;

    /* update random seed */
    if (lp->seedVal < 0)
      lp->seedVal += 0x10000;
    lp->seedVal = (lp->seedVal * 15625 + 1) & 0xFFFF;
    if (lp->seedVal >= 0x8000)
      lp->seedVal -= 0x10000;
    /* length of next segment in samples */
    lp->randLine_cnt = (int) ((p->sampleRate / reverbParams[n][2]) + 0.5);
    prvDel = (double) lp->writePos;
    prvDel -= ((double) lp->readPos
               + ((double) lp->readPosFrac / (double) DELAYPOS_SCALE));
    while (prvDel < 0.0)
      prvDel += (double) lp->bufferSize;
    prvDel = prvDel / p->sampleRate;    /* previous delay time in seconds */
    nxtDel = (double) lp->seedVal * reverbParams[n][1] / 32768.0;
    /* next delay time in seconds */
    nxtDel = reverbParams[n][0] + (nxtDel * (double) *(p->iPitchMod));
    /* calculate phase increment per sample */
    phs_incVal = (prvDel - nxtDel) / (double) lp->randLine_cnt;
    phs_incVal = phs_incVal * p->sampleRate + 1.0;
    lp->readPosFrac_inc = (int) (phs_incVal * DELAYPOS_SCALE + 0.5);
}

static void init_delay_line(SC_REVERB *p, delayLine *lp, int n)
{
    double  readPos;
    /* int     i; */

    /* calculate length of delay line */
    lp->bufferSize = delay_line_max_samples(p, n);
    lp->dummy = 0;
    lp->writePos = 0;
    /* set random seed */
    lp->seedVal = (int) (reverbParams[n][3] + 0.5);
    /* set initial delay time */
    readPos = (double) lp->seedVal * reverbParams[n][1] / 32768;
    readPos = reverbParams[n][0] + (readPos * (double) *(p->iPitchMod));
    readPos = (double) lp->bufferSize - (readPos * p->sampleRate);
    lp->readPos = (int) readPos;
    readPos = (readPos - (double) lp->readPos) * (double) DELAYPOS_SCALE;
    lp->readPosFrac = (int) (readPos + 0.5);
    /* initialise first random line segment */
    next_random_lineseg(p, lp, n);
    /* clear delay line to zero */
    lp->filterState = 0.0;
    memset(lp->buf, 0, sizeof(MYFLT)*lp->bufferSize);
    /* for (i = 0; i < lp->bufferSize; i++) */
    /*   lp->buf[i] = FL(0.0); */
}

static int sc_reverb_init(CSOUND *csound, SC_REVERB *p)
{
    int i;
    int nBytes;

    /* check for valid parameters */
    if (*(p->iSampleRate) <= FL(0.0))
      p->sampleRate = (double) csound->esr;
    else
      p->sampleRate = (double) *(p->iSampleRate);
    if (UNLIKELY(p->sampleRate < MIN_SRATE || p->sampleRate > MAX_SRATE)) {
      return csound->InitError(csound,
                               Str("reverbsc: sample rate is out of range"));
    }
    if (UNLIKELY(*(p->iPitchMod) < FL(0.0) || *(p->iPitchMod) > (MYFLT) MAX_PITCHMOD)) {
      return csound->InitError(csound,
                               Str("reverbsc: invalid pitch modulation factor"));
    }
    /* calculate the number of bytes to allocate */
    nBytes = 0;
    for (i = 0; i < 8; i++)
      nBytes += delay_line_bytes_alloc(p, i);
    if (nBytes != p->auxData.size)
      csound->AuxAlloc(csound, (size_t) nBytes, &(p->auxData));
    else if (p->initDone && *(p->iSkipInit) != FL(0.0))
      return OK;    /* skip initialisation if requested */
    /* set up delay lines */
    nBytes = 0;
    for (i = 0; i < 8; i++) {
      p->delayLines[i] = (delayLine*) ((unsigned char*) (p->auxData.auxp)
                                       + (int) nBytes);
      init_delay_line(p, p->delayLines[i], i);
      nBytes += delay_line_bytes_alloc(p, i);
    }
    p->dampFact = 1.0;
    p->prv_LPFreq = FL(0.0);
    p->initDone = 1;

    return OK;
}

static int sc_reverb_perf(CSOUND *csound, SC_REVERB *p)
{
    double    ainL, ainR, aoutL, aoutR;
    double    vm1, v0, v1, v2, am1, a0, a1, a2, frac;
    delayLine *lp;
    int       i, n, nsmps = csound->ksmps, readPos;
    int       bufferSize; /* Local copy */
    double    dampFact = p->dampFact;

    if (p->initDone <= 0) goto err1;
    /* calculate tone filter coefficient if frequency changed */
    if (*(p->kLPFreq) != p->prv_LPFreq) {
      p->prv_LPFreq = *(p->kLPFreq);
      dampFact = 2.0 - cos(p->prv_LPFreq * TWOPI / p->sampleRate);
      dampFact = p->dampFact = dampFact - sqrt(dampFact * dampFact - 1.0);
    }
    /* update delay lines */
    for (i = 0; i < nsmps; i++) {
      /* calculate "resultant junction pressure" and mix to input signals */
      ainL = aoutL = aoutR = 0.0;
      for (n = 0; n < 8; n++)
        ainL += p->delayLines[n]->filterState;
      ainL *= jpScale;
      ainR = ainL + (double) p->ainR[i];
      ainL = ainL + (double) p->ainL[i];
      /* loop through all delay lines */
      for (n = 0; n < 8; n++) {
        lp = p->delayLines[n];
        bufferSize = lp->bufferSize;
        /* send input signal and feedback to delay line */
        lp->buf[lp->writePos] = (MYFLT) ((n & 1 ? ainR : ainL)
                                         - lp->filterState);
        if (UNLIKELY(++lp->writePos >= bufferSize))
          lp->writePos -= bufferSize;
        /* read from delay line with cubic interpolation */
        if (lp->readPosFrac >= DELAYPOS_SCALE) {
          lp->readPos += (lp->readPosFrac >> DELAYPOS_SHIFT);
          lp->readPosFrac &= DELAYPOS_MASK;
        }
        if (UNLIKELY(lp->readPos >= bufferSize))
          lp->readPos -= bufferSize;
        readPos = lp->readPos;
        frac = (double) lp->readPosFrac * (1.0 / (double) DELAYPOS_SCALE);
        /* calculate interpolation coefficients */
        a2 = frac * frac; a2 -= 1.0; a2 *= (1.0 / 6.0);
        a1 = frac; a1 += 1.0; a1 *= 0.5; am1 = a1 - 1.0;
        a0 = 3.0 * a2; a1 -= a0; am1 -= a2; a0 -= frac;
        /* read four samples for interpolation */
        if (readPos > 0 && readPos < (bufferSize - 2)) {
          vm1 = (double) (lp->buf[readPos - 1]);
          v0  = (double) (lp->buf[readPos]);
          v1  = (double) (lp->buf[readPos + 1]);
          v2  = (double) (lp->buf[readPos + 2]);
        }
        else {
          /* at buffer wrap-around, need to check index */
          if (--readPos < 0) readPos += bufferSize;
          vm1 = (double) lp->buf[readPos];
          if (++readPos >= bufferSize) readPos -= bufferSize;
          v0 = (double) lp->buf[readPos];
          if (++readPos >= bufferSize) readPos -= bufferSize;
          v1 = (double) lp->buf[readPos];
          if (++readPos >= bufferSize) readPos -= bufferSize;
          v2 = (double) lp->buf[readPos];
        }
        v0 = (am1 * vm1 + a0 * v0 + a1 * v1 + a2 * v2) * frac + v0;
        /* update buffer read position */
        lp->readPosFrac += lp->readPosFrac_inc;
        /* apply feedback gain and lowpass filter */
        v0 *= (double) *(p->kFeedBack);
        v0 = (lp->filterState - v0) * dampFact + v0;
        lp->filterState = v0;
        /* mix to output */
        if (n & 1)
          aoutR += v0;
        else
          aoutL += v0;
        /* start next random line segment if current one has reached endpoint */
        if (--(lp->randLine_cnt) <= 0)
          next_random_lineseg(p, lp, n);
      }
      p->aoutL[i] = (MYFLT) (aoutL * outputGain);
      p->aoutR[i] = (MYFLT) (aoutR * outputGain);
    }

    return OK;
 err1:
    return csound->PerfError(csound, Str("reverbsc: not initialised"));
}

/* module interface functions */

int reverbsc_init_(CSOUND *csound)
{
    return csound->AppendOpcode(csound, "reverbsc",
                                (int) sizeof(SC_REVERB), 5, "aa", "aakkjpo",
                                (int (*)(CSOUND *, void *)) sc_reverb_init,
                                (int (*)(CSOUND *, void *)) NULL,
                                (int (*)(CSOUND *, void *)) sc_reverb_perf);
}


Generated by  Doxygen 1.6.0   Back to index