Logo Search packages:      
Sourcecode: csound version File versions  Download package

ugens5.c

/*
    ugens5.c:

    Copyright (C) 1991 Barry Vercoe, John ffitch, Gabriel Maldonado

    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 "csoundCore.h"         /*                      UGENS5.C        */
#include "ugens5.h"
#include <math.h>
#include "oload.h"

/*
 * LPC storage slots
 */

#define MAX_LPC_SLOT 20

int porset(CSOUND *csound, PORT *p)
{
    p->c2 = pow(0.5, (double)csound->onedkr / *p->ihtim);
    p->c1 = 1.0 - p->c2;
    if (*p->isig >= FL(0.0))
      p->yt1 = (double)(*p->isig);
    return OK;
}

int port(CSOUND *csound, PORT *p)
{
    p->yt1 = p->c1 * (double)*p->ksig + p->c2 * p->yt1;
    *p->kr = (MYFLT)p->yt1;
    return OK;
}

int tonset(CSOUND *csound, TONE *p)
{
    double b;
    p->prvhp = (double)*p->khp;
    b = 2.0 - cos((double)(p->prvhp * csound->tpidsr));
    p->c2 = b - sqrt(b * b - 1.0);
    p->c1 = 1.0 - p->c2;

    if (!(*p->istor))
      p->yt1 = 0.0;
    return OK;
}

int tone(CSOUND *csound, TONE *p)
{
    MYFLT       *ar, *asig;
    int         n, nsmps = csound->ksmps;
    double      c1 = p->c1, c2 = p->c2;
    double      yt1 = p->yt1;

    if (*p->khp != (MYFLT)p->prvhp) {
      double b;
      p->prvhp = (double)*p->khp;
      b = 2.0 - cos((double)(p->prvhp * csound->tpidsr));
      p->c2 = c2 = b - sqrt(b * b - 1.0);
      p->c1 = c1 = 1.0 - c2;
    }
    ar = p->ar;
    asig = p->asig;
    for (n=0; n<nsmps; n++) {
      yt1 = c1 * (double)(asig[n]) + c2 * yt1;
      ar[n] = (MYFLT)yt1;
    }
    p->yt1 = yt1;
    return OK;
}

int tonsetx(CSOUND *csound, TONEX *p)
{                   /* From Gabriel Maldonado, modified for arbitrary order */
    {
      double b;
      p->prvhp = *p->khp;
      b = 2.0 - cos((double)(*p->khp * csound->tpidsr));
      p->c2 = b - sqrt(b * b - 1.0);
      p->c1 = 1.0 - p->c2;
    }
    if ((p->loop = (int) (*p->ord + FL(0.5))) < 1) p->loop = 4;
    if (!*p->istor && (p->aux.auxp == NULL ||
                       (int)(p->loop*sizeof(double)) > p->aux.size))
      csound->AuxAlloc(csound, (long)(p->loop*sizeof(double)), &p->aux);
    p->yt1 = (double*)p->aux.auxp;
    if (!(*p->istor)) {
      memset(p->yt1, 0, p->loop*sizeof(double)); /* Punning zero and 0.0 */
    }
    return OK;
}

int tonex(CSOUND *csound, TONEX *p)      /* From Gabriel Maldonado, modified */
{
    int j;
    int nsmps;
    MYFLT *asig, *ar;
    double c1, c2, *yt1;
    if (*p->khp != p->prvhp) {
      double b;
      p->prvhp = (double)*p->khp;
      b = 2.0 - cos(p->prvhp * (double)csound->tpidsr);
      p->c2 = b - sqrt(b * b - 1.0);
      p->c1 = 1.0 - p->c2;
    }
    c1 = p->c1;
    c2 = p->c2;
    yt1= p->yt1;
    asig = p->asig;
    nsmps = csound->ksmps;
    ar = p->ar;
    for (j=0; j< p->loop; j++) {
      int n;
      for (n=0; n<nsmps; n++) {
        double x = c1 * asig[n] + c2 * *yt1;
        *yt1 = x;
        ar[n] = (MYFLT)x;
      }
      asig = p->ar;
      yt1++;
    }
    return OK;
}

int atone(CSOUND *csound, TONE *p)
{
    MYFLT       *ar, *asig;
    int n, nsmps = csound->ksmps;
    double      c2 = p->c2, yt1 = p->yt1;

    if (*p->khp != p->prvhp) {
      double b;
      p->prvhp = *p->khp;
      b = 2.0 - cos((double)(*p->khp * csound->tpidsr));
      p->c2 = c2 = b - sqrt(b * b - 1.0);
/*      p->c1 = c1 = 1.0 - c2; */
    }
    ar = p->ar;
    asig = p->asig;
    for (n=0; n<nsmps; n++) {
      double sig = (double)asig[n];
      double x = yt1 = c2 * (yt1 + sig);
      ar[n] = (MYFLT)x;
      yt1 -= sig;               /* yt1 contains yt1-xt1 */
    }
    p->yt1 = yt1;
    return OK;
}

int atonex(CSOUND *csound, TONEX *p)      /* Gabriel Maldonado, modified */
{
    MYFLT       *ar = p->ar, *asig;
    double      c2, *yt1;
    int         n, nsmps=csound->ksmps, j;

    if (*p->khp != p->prvhp) {
      double b;
      p->prvhp = *p->khp;
      b = 2.0 - cos((double)(*p->khp * csound->tpidsr));
      p->c2 = b - sqrt(b * b - 1.0);
      /*p->c1 = 1. - p->c2;*/
    }

    c2 = p->c2;
    yt1=p->yt1;
    asig = p->asig;
    for (j=0; j< p->loop; j++) {
      for (n=0; n<nsmps; n++) {
        double sig = asig[n];
        double x = c2 * (yt1[j] + sig);
        yt1[j] = x - sig;            /* yt1 contains yt1-xt1 */
        ar[n] = (MYFLT)x;
      }
    }
    return OK;
}

int rsnset(CSOUND *csound, RESON *p)
{
    int scale;
    p->scale = scale = (int)*p->iscl;
    if (scale && scale != 1 && scale != 2) {
      return csound->InitError(csound, Str("illegal reson iscl value, %f"),
                                       *p->iscl);
    }
    p->prvcf = p->prvbw = -100.0;
    if (!(*p->istor))
      p->yt1 = p->yt2 = 0.0;
    return OK;
}

int reson(CSOUND *csound, RESON *p)
{
    int flag = 0, n, nsmps = csound->ksmps;
    MYFLT       *ar, *asig;
    double      c3p1, c3t4, omc3, c2sqr;
    double      yt1, yt2, c1 = p->c1, c2 = p->c2, c3 = p->c3;

    if (*p->kcf != (MYFLT)p->prvcf) {
      p->prvcf = (double)*p->kcf;
      p->cosf = cos(p->prvcf * (double)(csound->tpidsr));
      flag = 1;                 /* Mark as changed */
    }
    if (*p->kbw != (MYFLT)p->prvbw) {
      p->prvbw = (double)*p->kbw;
      c3 = p->c3 = exp(p->prvbw * (double)(csound->mtpdsr));
      flag = 1;                /* Mark as changed */
    }
    if (flag) {
      c3p1 = c3 + 1.0;
      c3t4 = c3 * 4.0;
      omc3 = 1.0 - c3;
      c2 = p->c2 = c3t4 * p->cosf / c3p1;               /* -B, so + below */
      c2sqr = c2 * c2;
      if (p->scale == 1)
        c1 = p->c1 = omc3 * sqrt(1.0 - c2sqr / c3t4);
      else if (p->scale == 2)
        c1 = p->c1 = sqrt((c3p1*c3p1-c2sqr) * omc3/c3p1);
      else c1 = p->c1 = 1.0;
    }
    asig = p->asig;
    ar = p->ar;
    yt1 = p->yt1; yt2 = p->yt2;
    for (n=0; n<nsmps; n++) {
      double yt0 = c1 * (double)(asig[n]) + c2 * yt1 - c3 * yt2;
      ar[n] = (MYFLT)yt0;
      yt2 = yt1;
      yt1 = yt0;
    }
    p->yt1 = yt1; p->yt2 = yt2; /* Write back for next cycle */
    return OK;
}

int rsnsetx(CSOUND *csound, RESONX *p)
{                               /* Gabriel Maldonado, modifies for arb order */
    int scale;
    p->scale = scale = (int) *p->iscl;
    if ((p->loop = (int) (*p->ord + FL(0.5))) < 1)
      p->loop = 4; /* default value */
    if (!*p->istor && (p->aux.auxp == NULL ||
                       (int)(p->loop*2*sizeof(double)) > p->aux.size))
      csound->AuxAlloc(csound, (long)(p->loop*2*sizeof(double)), &p->aux);
    p->yt1 = (double*)p->aux.auxp; p->yt2 = (double*)p->aux.auxp + p->loop;
    if (scale && scale != 1 && scale != 2) {
      return csound->InitError(csound, Str("illegal reson iscl value, %f"),
                                       *p->iscl);
    }
    p->prvcf = p->prvbw = -100.0;

    if (!(*p->istor)) {
      memset(p->yt1, 0, p->loop*sizeof(double));
      memset(p->yt2, 0, p->loop*sizeof(double));
    }
    return OK;
}

int resonx(CSOUND *csound, RESONX *p)   /* Gabriel Maldonado, modified  */
{
    int flag = 0, nsmps = csound->ksmps, j;
    MYFLT       *ar, *asig;
    double       c3p1, c3t4, omc3, c2sqr;
    double       *yt1, *yt2, c1,c2,c3;

    if (*p->kcf != (MYFLT)p->prvcf) {
      p->prvcf = (double)*p->kcf;
      p->cosf = cos(p->prvcf * (double)(csound->tpidsr));
      flag = 1;
    }
    if (*p->kbw != (MYFLT)p->prvbw) {
      p->prvbw = (double)*p->kbw;
      p->c3 = exp(p->prvbw * (double)(csound->mtpdsr));
      flag = 1;
    }
    if (flag) {
      c3p1 = p->c3 + 1.0;
      c3t4 = p->c3 * 4.0;
      omc3 = 1.0 - p->c3;
      p->c2 = c3t4 * p->cosf / c3p1;            /* -B, so + below */
      c2sqr = p->c2 * p->c2;
      if (p->scale == 1)
        p->c1 = omc3 * sqrt(1.0 - (c2sqr / c3t4));
      else if (p->scale == 2)
        p->c1 = sqrt((c3p1*c3p1-c2sqr) * omc3/c3p1);
      else p->c1 = 1.0;
    }

    ar   = p->ar;
    c1   = p->c1;
    c2   = p->c2;
    c3   = p->c3;
    yt1  = p->yt1;
    yt2  = p->yt2;
    ar = p->ar;
    asig = p->asig;
    for (j=0; j< p->loop; j++) {
      int n;
      for (n=0; n<nsmps; n++) {
        double x =
          c1 * (double)asig[n] + c2 * yt1[j] - c3 * yt2[j];
        yt2[j] = yt1[j];
        ar[n] = (MYFLT)x;
        yt1[j] = x;
      }
    }
    return OK;
}

int areson(CSOUND *csound, RESON *p)
{
    int flag = 0, nsmps = csound->ksmps;
    MYFLT       *ar, *asig;
    double      c3p1, c3t4, omc3, c2sqr, D = 2.0; /* 1/RMS = root2 (rand) */
                                                   /*      or 1/.5  (sine) */
    double      yt1, yt2, c1, c2, c3;

    if (*p->kcf != (MYFLT)p->prvcf) {
      p->prvcf = (double)*p->kcf;
      p->cosf = cos(p->prvcf * (double)(csound->tpidsr));
      flag = 1;
    }
    if (*p->kbw != (MYFLT)p->prvbw) {
      p->prvbw = (double)*p->kbw;
      p->c3 = exp(p->prvbw * (double)(csound->mtpdsr));
      flag = 1;
    }
    if (flag) {
      c3p1 = p->c3 + 1.0;
      c3t4 = p->c3 * 4.0;
      omc3 = 1.0 - p->c3;
      p->c2 = c3t4 * p->cosf / c3p1;
      c2sqr = p->c2 * p->c2;
      if (p->scale == 1)                        /* i.e. 1 - A(reson) */
        p->c1 = 1.0 - omc3 * sqrt(1.0 - c2sqr / c3t4);
      else if (p->scale == 2)                 /* i.e. D - A(reson) */
        p->c1 = D - sqrt((c3p1*c3p1-c2sqr)*omc3/c3p1);
      else p->c1 = 0.0;                        /* cannot tell        */
    }
    asig = p->asig;
    ar = p->ar;
    c1 = p->c1; c2 = p->c2; c3 = p->c3; yt1 = p->yt1; yt2 = p->yt2;
    if (p->scale == 1 || p->scale == 0) {
      int n;
      for (n=0; n<nsmps; n++) {
        double sig = (double)asig[n];
        double ans = c1 * sig + c2 * yt1 - c3 * yt2;
        yt2 = yt1;
        yt1 = ans - sig;  /* yt1 contains yt1-xt1 */
        ar[n] = (MYFLT)ans;
      }
    }
    else if (p->scale == 2) {
      int n;
      for (n=0; n<nsmps; n++) {
        double sig = (double)asig[n];
        double ans = c1 * sig + c2 * yt1 - c3 * yt2;
        yt2 = yt1;
        yt1 = ans - D * sig;      /* yt1 contains yt1-D*xt1 */
        ar[n] = (MYFLT)ans;
      }
    }
    p->yt1 = yt1; p->yt2 = yt2;
    return OK;
}

/*
 *
 * LPREAD opcode : initialisation phase
 *
 *
 */

int lprdset(CSOUND *csound, LPREAD *p)
{
    LPHEADER *lph;
    MEMFIL   *mfp;
    long     magic;
    long     totvals;  /* NB - presumes sizeof(MYFLT) == sizeof(long) !! */
    char     lpfilname[MAXNAME];

    /* Store adress of opcode for other lpXXXX init to point to */
    if (csound->lprdaddr == NULL ||
        csound->currentLPCSlot >= csound->max_lpc_slot) {
      csound->max_lpc_slot = csound->currentLPCSlot + MAX_LPC_SLOT;
      csound->lprdaddr = mrealloc(csound,
                                  csound->lprdaddr,
                                  csound->max_lpc_slot * sizeof(LPREAD*));
    }
    ((LPREAD**) csound->lprdaddr)[csound->currentLPCSlot] = p;

    /* Build file name */
    csound->strarg2name(csound, lpfilname, p->ifilno, "lp.", p->XSTRCODE);

    /* Do not reload existing file ? */
    if ((mfp = p->mfp) != NULL && strcmp(mfp->filename, lpfilname) == 0)
      goto lpend;                             /* rtn if file prv known */
    /* Load analysis in memory file */
    /* else read file  */
    if ((mfp = ldmemfile2(csound, lpfilname, CSFTYPE_LPC)) == NULL) {
      return csound->InitError(csound, Str("LPREAD cannot load %s"), lpfilname);
    }
    /* Store memory file location in opcode */
    p->mfp = mfp;                                   /*  & record facts   */
    /* Take a peek to the header if exisiting. Else take input arguments */
    lph = (LPHEADER *) mfp->beginp;

    magic=lph->lpmagic;
    if ((magic==LP_MAGIC)||(magic==LP_MAGIC2)) {
      p->storePoles = (magic==LP_MAGIC2);

      csound->Message(csound, Str("Using %s type of file.\n"),
                      p->storePoles?Str("pole"):Str("filter coefficient"));
      /* Store header length */
      p->headlongs = lph->headersize/sizeof(long);
      /* Check if input values where available */
      if (*p->inpoles || *p->ifrmrate) {
        csound->Warning(csound, Str("lpheader overriding inputs"));
      }
      /* Check orc/analysis sample rate compatibility */
      if (lph->srate != csound->esr) {
        csound->Warning(csound, Str("lpfile srate != orch sr"));
      }
      p->npoles = lph->npoles;                /* note npoles, etc. */
      /* Store header info in opcode */
      p->nvals = lph->nvals;
      p->framrat16 = lph->framrate * FL(65536.0);/* scaled framno cvt */
    }
    else if (BYTREVL(lph->lpmagic) == LP_MAGIC) {   /* Header reversed:  */
      return csound->InitError(csound, Str("file %s bytes are in wrong order"),
                                       lpfilname);
    }
    else {                                    /* No Header on file:*/
      p->headlongs = 0;
      p->npoles = (long)*p->inpoles;          /*  data from inargs */
      p->nvals = p->npoles + 4;
      p->framrat16 = *p->ifrmrate * FL(65536.0);
      if (!p->npoles || !p->framrat16) {
        return csound->InitError(csound,
                                 Str("insufficient args and no file header"));
      }
    }
    /* Check  pole number */
    if (p->npoles > MAXPOLES) {
      return csound->InitError(csound, Str("npoles > MAXPOLES"));
    }
    /* Look for total frame data size (file size - header) */
    totvals = (mfp->length/sizeof(long)) - p->headlongs;   /* see NB above!! */
    /* Store the size of a frame in integer */
    p->lastfram16 = (((totvals - p->nvals) / p->nvals) << 16) - 1;
    if (csound->oparms->odebug)
      csound->Message(csound, Str(
                 "npoles %ld, nvals %ld, totvals %ld, lastfram16 = %lx\n"),
             p->npoles, p->nvals, totvals, p->lastfram16);
 lpend:
    p->lastmsg = 0;
    return OK;
}

#ifdef TRACE_POLES
static void
 DumpPolesF(int poleCount, MYFLT *part1, MYFLT *part2, int isMagn, char *where)
{
    int i;

    csound->Message(csound, "%s\n", where);
    for (i=0; i<poleCount; i++) {
      if (isMagn)
        csound->Message(csound, Str("magnitude: %f   Phase: %f\n"),
                                part1[i], part2[i]);
      else
        csound->Message(csound, Str("Real: %f   Imag: %f\n"),
                                part1[i], part2[i]);
    }
}
#endif

static void SortPoles(int poleCount, MYFLT *poleMagn, MYFLT *polePhas)
{
    int i, j;
    MYFLT diff, fTemp;
    int shouldSwap;

    /*  DumpPolesF(poleCount, poleMagn, polePhas, 1, "Before sort"); */

    for (i=1; i<poleCount; i++) {
      for (j=0; j<i; j++) {

        shouldSwap = 0;

        diff = (MYFLT)(fabs(polePhas[j])-fabs(polePhas[i]));
        if (diff>1.0e-10)
          shouldSwap = 1;
        else if (diff>-1.0e-10) {
          diff = poleMagn[j]-poleMagn[i];

          if (diff>1.0e-10)
            shouldSwap = 1;
          else if (diff>-1.0e-10)
            {
              if (polePhas[j]>polePhas[i])
                shouldSwap = 1;
            }
        }
        if (shouldSwap) {
          fTemp = poleMagn[i];
          poleMagn[i] = poleMagn[j];
          poleMagn[j] = fTemp;

          fTemp = polePhas[i];
          polePhas[i] = polePhas[j];
          polePhas[j] = fTemp;
        }
      }
    }
/*  DumpPolesF(poleCount, poleMagn, polePhas, 1, "After sort"); */
}

static int DoPoleInterpolation(int poleCount,
                               MYFLT *pm1, MYFLT *pp1,
                               MYFLT *pm2, MYFLT *pp2,
                               MYFLT factor, MYFLT *outMagn, MYFLT *outPhas)
{
    int i;

    if (poleCount%2!=0) {
/*    printf (Str("Cannot handle uneven pole count yet \n")); */
      return (0);
    }

    for (i=0; i<poleCount; i++) {
      if (fabs(fabs(pp1[i])-PI)<1.0e-5) {
        pm1[i] = -pm1[i];
        pp1[i] = FL(0.0);
      }

      if (fabs(fabs(pp2[i])-PI)<1.0e-5) {
        pm2[i] = -pm2[i];
        pp2[i] = FL(0.0);
      }
    }

    /* Sort poles according to abs(phase) */

    SortPoles(poleCount, pm1, pp1);
    SortPoles(poleCount, pm2, pp2);

        /*      DumpPolesF(poleCount, pm1, pp1, 1, "Sorted poles 1"); */
        /*      DumpPolesF(poleCount, pm2, pp2, 1, "Sorted poles 2"); */

        /*      printf ("factor=%f\n", factor); */

    for (i=0; i<poleCount; i++) {
      outMagn[i] = pm1[i]+(pm2[i]-pm1[i])*factor;
      outPhas[i] = pp1[i]+(pp2[i]-pp1[i])*factor;
    }

/*     DumpPolesF(poleCount, outMagn, outPhas, 1, "Interpolated poles"); */

    return(1);
}

static inline void InvertPoles(int count, double *real, double *imag)
{
    int    i;
    double pr,pi,mag;

    for (i=0; i<count; i++) {
      pr = real[i];
      pi = imag[i];
      mag = pr*pr+pi*pi;
      real[i] = pr/mag;
      imag[i] = -pi/mag;
    }
}

/*
 *
 * Resynthetize filter coefficients from poles values
 *
 */

static inline void
    synthetize(int    poleCount,
               double *poleReal, double *poleImag,
               double *polyReal, double *polyImag)
{
    int    j, k;
    double pr, pi, cr, ci;

    polyReal[0] = 1;
    polyImag[0] = 0;

    for (j=0; j<poleCount; j++) {
      polyReal[j+1] = 1;
      polyImag[j+1] = 0;

      pr = poleReal[j];
      pi = poleImag[j];

      for (k=j; k>=0; k--) {
        cr = polyReal[k];
        ci = polyImag[k];

        polyReal[k] = -(cr*pr-ci*pi);
        polyImag[k] = -(ci*pr+cr*pi);

        if (k>0) {
            polyReal[k] += polyReal[k-1];
            polyImag[k] += polyImag[k-1];
        }
      }
    }

    /* Makes it 1+a1.x+...+anXn */

    pr = polyReal[0];
    for (j=0; j<=poleCount; j++)
      polyReal[j] /= pr;
}

/*
 *
 * LPREAD k/a time access. This will setup current pole values
 *
 */

int lpread(CSOUND *csound, LPREAD *p)
{
    MYFLT   *bp, *np, *cp;
    long    nn, framphase;
    MYFLT   fract;
    int     i, status;
    MYFLT   poleMagn1[MAXPOLES], polePhas1[MAXPOLES];
    MYFLT   poleMagn2[MAXPOLES], polePhas2[MAXPOLES];
    MYFLT   interMagn[MAXPOLES], interPhas[MAXPOLES];

    if (p->mfp==NULL) {
      return csound->PerfError(csound, Str("lpread: not initialised"));
    }
    /* Locate frame position range */
    if ((framphase = (long)(*p->ktimpt*p->framrat16)) < 0) { /* for kfram reqd*/
      return csound->PerfError(csound, Str("lpread timpnt < 0"));
    }
    if (framphase > p->lastfram16) {                /* not past last one */
      framphase = p->lastfram16;
      if (!p->lastmsg) {
        p->lastmsg = 1;
        csound->Warning(csound, Str("lpread ktimpnt truncated to last frame"));
      }
    }
    /* Locate frames bounding current time */
    nn = (framphase >> 16) * p->nvals + p->headlongs;   /* see NB above!! */
    bp = (MYFLT *)p->mfp->beginp + nn;          /* locate begin this frame */
    np = bp + p->nvals;                         /* & interp betw adj frams */
    fract = (framphase & 0x0FFFFL) / FL(65536.0);
    /* Interpolate freq/amplpitude and store in opcode */
    *p->krmr = *bp + (*np - *bp) * fract;   bp++;   np++; /* for 4 rslts */
    *p->krmo = *bp + (*np - *bp) * fract;   bp++;   np++;
    *p->kerr = *bp + (*np - *bp) * fract;   bp++;   np++;
    *p->kcps = *bp + (*np - *bp) * fract;   bp++;   np++;

   /* Interpolate filter or poles coef values and store in opcode */

    cp = p->kcoefs;      /* This is where the coefs get stored */
    if (p->storePoles) {
      for (i=0; i<p->npoles; i++) {
        poleMagn1[i] = *bp++;
        polePhas1[i] = *bp++;
        poleMagn2[i] = *np++;
        polePhas2[i] = *np++;
      }

      status =
        DoPoleInterpolation(p->npoles,poleMagn1,polePhas1,poleMagn2,
                            polePhas2,fract,interMagn,interPhas);
      if (!status) {
        return csound->PerfError(csound, Str("Interpolation failed"));
      }
      for (i=0; i<p->npoles; i++) {
        *cp++ = interMagn[i];
        *cp++ = interPhas[i];
      }
    }
    else {
      nn = p->npoles;
      do {
        *cp = *bp + (*np - *bp) * fract;
        cp++; bp++; np++;
      }
      while (--nn);
    }
/*  if (csound->oparms->odebug) {
      csound->Message(csound, "phase:%lx fract:%6.2f rmsr:%6.2f rmso:%6.2f kerr:%6.2f kcps:%6.2f\n",
             framphase,fract,*p->krmr,*p->krmo,*p->kerr,*p->kcps);
      cp = p->kcoefs;
      nn = p->npoles;
      do {
        csound->Message(csound, " %6.2f",*cp++);
      } while (--nn);
      csound->Message(csound, "\n");
    }  */
    return OK;
}

/*
 *
 * LPRESON: initialisation time
 *
 *
 */
int lprsnset(CSOUND *csound, LPRESON *p)
{
    LPREAD *q;

   /* connect to previously loaded lpc analysis */
   /* get adr lpread struct */
    p->lpread = q = ((LPREAD**) csound->lprdaddr)[csound->currentLPCSlot];

   /* Initialize pointer to circulat buffer (for filtering) */
    p->circjp = p->circbuf;
    p->jp2lim = p->circbuf + (q->npoles << 1);  /* npoles det circbuflim */
    return OK;
}

/*
 *
 * LPRESON: k & a time access. Will actually filter the signal
 *                  Uses a circular buffer to store previous signal values.
 */

int lpreson(CSOUND *csound, LPRESON *p)
{
    LPREAD *q = p->lpread;
    int     nn, nsmps = csound->ksmps;
    MYFLT   *coefp, *pastp, *jp, *jp2, *rslt = p->ar, *asig = p->asig;
    MYFLT   x;
    double  poleReal[MAXPOLES], poleImag[MAXPOLES];
    double  polyReal[MAXPOLES+1], polyImag[MAXPOLES+1];
    int     i;
    double  pm,pp;

    jp = p->circjp;
    jp2 = jp + q->npoles;

    /* If we where using poles, we have to compute filter coefs now */
    if (q->storePoles) {
      coefp = q->kcoefs;
      for (i=0; i<q->npoles; i++) {
        pm = *coefp++;
        pp = *coefp++;
        if (fabs(pm)>0.999999)
          pm = 1/pm;
        poleReal[i] = pm*cos(pp);
        poleImag[i] = pm*sin(pp);
      }
      /*     DumpPoles(q->npoles,poleReal,poleImag,0,"About to filter"); */
      InvertPoles(q->npoles,poleReal,poleImag);
      synthetize(q->npoles,poleReal,poleImag,polyReal,polyImag);
      coefp = q->kcoefs;
      for (i=0; i<q->npoles; i++) {
        coefp[i] = -(MYFLT)polyReal[q->npoles-i]; /* MR_WHY - somthing with the atan2 ? */
#ifdef _DEBUG
/*                      if (polyImag[i]>1.0e-10) */
/*                      { */
/*                              printf ("bad polymag: %f\n",polyImag[i]); */
/*                      } */
#endif
      }
    }

    /* For each sample */
    do {
      /* Compute Xn = Yn + CkXn-k */

#ifdef TRACE_FILTER
      csound->Message(csound, "Asig=%f\n", *asig);
#endif
      x = *asig++;
      coefp = q->kcoefs;              /* using lpread interp coefs */
      pastp = jp;
      nn = q->npoles;
      do {
#ifdef TRACE_FILTER
        csound->Message(csound, "\t%f,%f\n", *coefp, *pastp);
#endif
        x += *coefp++ * *pastp++;
      } while (--nn);
#ifdef TRACE_FILTER
      csound->Message(csound, "result=%f\n", x);
#endif
      /* Store result signal in circular and output buffers */

      *jp++ = *jp2++ = x;
      *rslt++ = x;

      /* Check if end of buffer reached */
      if (jp2 >= p->jp2lim) {
        jp2 = jp;
        jp = p->circbuf;
      }
    } while (--nsmps);
    p->circjp = jp;
    return OK;
}

/*
 *
 * LPFRESON : Initialisation time
 *
 */
int lpfrsnset(CSOUND *csound, LPFRESON *p)
{

   /* Connect to previously loaded analysis file */

    if (((LPREAD**) csound->lprdaddr)[csound->currentLPCSlot]->storePoles) {
      return csound->InitError(csound, Str("Pole file not supported "
                                           "for this opcode !"));
    }
    p->lpread = ((LPREAD**) csound->lprdaddr)[csound->currentLPCSlot];
    p->prvratio = FL(1.0);
    p->d = FL(0.0);
    p->prvout = FL(0.0);
    return OK;
}

/*
 *
 * LPFRESON : k & a time : actually filters the data
 *
 */
int lpfreson(CSOUND *csound, LPFRESON *p)
{
    LPREAD  *q = p->lpread;
    int     nn, nsmps = csound->ksmps;
    MYFLT   *coefp, *pastp, *pastp1, *rslt = p->ar, *asig = p->asig;
    MYFLT   x, temp1, temp2, ampscale, cq;

    if (*p->kfrqratio != p->prvratio) {             /* for new freqratio */
      if (*p->kfrqratio <= FL(0.0)) {
        return csound->PerfError(csound, Str("illegal frqratio, %5.2f"),
                                         *p->kfrqratio);
      }                                             /*      calculate d  */
      p->d = (*p->kfrqratio - FL(1.0)) / (*p->kfrqratio + FL(1.0));
      p->prvratio = *p->kfrqratio;
    }
    if (p->d != FL(0.0)) {                          /* for non-zero d,   */
      coefp = q->kcoefs;
      nn = q->npoles - 1;
      do {
        temp1 = p->d * *coefp++;                    /*    shift formants */
        *coefp += temp1;
      }
      while (--nn);
      ampscale = FL(1.0) / (FL(1.0) - p->d * *coefp); /*    & reset scales */
      cq = (FL(1.0) - p->d * p->d) * ampscale;
    }
    else {
      cq = FL(1.0);
      ampscale = FL(1.0);
    }
    x = p->prvout;
    do {
      nn = q->npoles - 1;
      pastp  = pastp1 = p->past + nn;
      temp1 = *pastp;
      *pastp = cq * x - p->d * *pastp;
      pastp--;
      do {
        temp2 = *pastp;
        *pastp = (*pastp1 - *pastp) * p->d + temp1;
        pastp--;   pastp1--;
        temp1 = temp2;
      }
      while (--nn);
      x = *asig++;
      pastp = p->past;
      coefp = q->kcoefs;
      nn = q->npoles;
      do  x += *coefp++ * *pastp++;
      while (--nn);
      *rslt++ = x * ampscale;
    }
    while (--nsmps);
    p->prvout = x;
    return OK;
}

int rmsset(CSOUND *csound, RMS *p)
{
    double   b;

    b = 2.0 - cos((double)(*p->ihp * csound->tpidsr));
    p->c2 = b - sqrt(b*b - 1.0);
    p->c1 = 1.0 - p->c2;
    if (!*p->istor)
      p->prvq = 0.0;
    return OK;
}

int gainset(CSOUND *csound, GAIN *p)
{
    double   b;

    b = 2.0 - cos((double)(*p->ihp * csound->tpidsr));
    p->c2 = b - sqrt(b*b - 1.0);
    p->c1 = 1.0 - p->c2;
    if (!*p->istor)
      p->prvq = p->prva = 0.0;
    return OK;
}

int balnset(CSOUND *csound, BALANCE *p)
{
    double   b;

    b = 2.0 - cos((double)(*p->ihp * csound->tpidsr));
    p->c2 = b - sqrt(b*b - 1.0);
    p->c1 = 1.0 - p->c2;
    if (!*p->istor)
      p->prvq = p->prvr = p->prva = 0.0;
    return OK;
}

int rms(CSOUND *csound, RMS *p)
{
    int     n, nsmps = csound->ksmps;
    MYFLT   *asig;
    double  q;
    double  c1 = p->c1, c2 = p->c2;

    q = p->prvq;
    asig = p->asig;
    for (n=0; n<nsmps; n++) {
      double as = (double)asig[n];
      q = c1 * as * as + c2 * q;
    }
    p->prvq = q;
    *p->kr = (MYFLT) sqrt(q);
    return OK;
}

int gain(CSOUND *csound, GAIN *p)
{
    int     nsmps = csound->ksmps;
    MYFLT   *ar, *asig;
    double  q, a, m, diff, inc;
    double  c1 = p->c1, c2 = p->c2;
    int     n;

    q = p->prvq;
    asig = p->asig;
    for (n = 0; n < nsmps; n++) {
      double as = (double)asig[n];
      q = c1 * as * as + c2 * q;
    }
    p->prvq = q;
    if (q > 0.0)
      a = *p->krms / sqrt(q);
    else
      a = *p->krms;
    ar = p->ar;
    if ((diff = a - p->prva) != 0.0) {
      m = p->prva;
      inc = diff * (double)csound->onedksmps;
      for (n = 0; n < nsmps; n++) {
        ar[n] = asig[n] * m;
        m += inc;
      }
      p->prva = a;
    }
    else {
      for (n = 0; n < nsmps; n++) {
        ar[n] = asig[n] * a;
      }
    }
    return OK;
}

int balance(CSOUND *csound, BALANCE *p)
{
    int     n, nsmps = csound->ksmps;
    MYFLT   *ar, *asig, *csig;
    double  q, r, a, m, diff, inc;
    double  c1 = p->c1, c2 = p->c2;

    q = p->prvq;
    r = p->prvr;
    asig = p->asig;
    csig = p->csig;
    for (n = 0; n < nsmps; n++) {
      double as = (double)asig[n];
      double cs = (double)csig[n];
      q = c1 * as * as + c2 * q;
      r = c1 * cs * cs + c2 * r;
    }
    p->prvq = q;
    p->prvr = r;
    if (q != 0.0)
      a = sqrt(r/q);
    else
      a = sqrt(r);
    ar = p->ar;
    if ((diff = a - p->prva) != 0.0) {
      m = p->prva;
      inc = diff * (double)csound->onedksmps;
      for (n = 0; n < nsmps; n++) {
        ar[n] = asig[n] * m;
        m += inc;
      }
      p->prva = a;
    }
    else {
      for (n = 0; n < nsmps; n++) {
        ar[n] = asig[n] * a;
      }
    }
    return OK;
}

/*
 *   Set current lpc slot
 */
int lpslotset(CSOUND *csound, LPSLOT *p)
{
    int n;

    n = (int) *(p->islotnum);
    if (n < 0)
      return csound->InitError(csound, Str("lpslot number should be positive"));
    else {
      if (n >= csound->max_lpc_slot) {
        csound->max_lpc_slot = n + MAX_LPC_SLOT;
        csound->lprdaddr = mrealloc(csound,
                                    csound->lprdaddr,
                                    csound->max_lpc_slot * sizeof(LPREAD**));
      }
      csound->currentLPCSlot = n;
    }
    return OK;
}

int lpitpset(CSOUND *csound, LPINTERPOL *p)
{

    if ((unsigned int) ((int) *(p->islot1))
          >= (unsigned int) csound->max_lpc_slot ||
        (unsigned int) ((int) *(p->islot2))
          >= (unsigned int) csound->max_lpc_slot)
      return csound->InitError(csound, Str("LPC slot is not allocated"));
  /* Get lpread pointers */
    p->lp1 = ((LPREAD**) csound->lprdaddr)[(int) *(p->islot1)];
    p->lp2 = ((LPREAD**) csound->lprdaddr)[(int) *(p->islot2)];

  /* Check if workable */

    if ((!p->lp1->storePoles) || (!p->lp2->storePoles)) {
      return csound->InitError(csound, Str("lpinterpol works only "
                                           "with poles files.."));
    }
    if (p->lp1->npoles != p->lp2->npoles) {
      return csound->InitError(csound, Str("The poles files "
                                           "have different pole count"));
    }

#if 0                   /* This is incorrect C */
    if (&p->kcoefs-p != &p->lp1->kcoefs-p->lp1)
      return csound->InitError(csound, Str("padding error"));
#endif

    p->npoles = p->lp1->npoles;
    p->storePoles = 1;
    ((LPREAD**) csound->lprdaddr)[csound->currentLPCSlot] = (LPREAD*) p;
    return OK;
}

int lpinterpol(CSOUND *csound, LPINTERPOL *p)
{
    int     i,status;
    MYFLT   *cp,*cp1,*cp2;
    MYFLT   poleMagn1[MAXPOLES], polePhas1[MAXPOLES];
    MYFLT   poleMagn2[MAXPOLES], polePhas2[MAXPOLES];
    MYFLT   interMagn[MAXPOLES], interPhas[MAXPOLES];

    /* RWD: guessing this... */
    if (p->lp1==NULL || p->lp2==NULL) {
      return csound->PerfError(csound, Str("lpinterpol: not initialised"));
    }
    cp1 =  p->lp1->kcoefs;
    cp2 =  p->lp2->kcoefs;

    for (i=0; i<p->npoles; i++) {
      poleMagn1[i] = *cp1++;
      polePhas1[i] = *cp1++;
      poleMagn2[i] = *cp2++;
      polePhas2[i] = *cp2++;
    }

    status = DoPoleInterpolation(p->npoles,poleMagn1,polePhas1,poleMagn2,
                                     polePhas2,*p->kmix,interMagn,interPhas);
    if (!status) {
      return csound->PerfError(csound, Str("Interpolation failed"));
    }

    cp = p->kcoefs;      /* This is where the coefs get stored */
    for (i=0; i<p->npoles; i++) {
      *cp++ = interMagn[i];
      *cp++ = interPhas[i];
    }
    return OK;
}


Generated by  Doxygen 1.6.0   Back to index