src/bend.c

   1 /* libSoX effect: Pitch Bend   (c) 2008 robs@users.sourceforge.net
   2  *
   3  * This library is free software; you can redistribute it and/or modify it
   4  * under the terms of the GNU Lesser General Public License as published by
   5  * the Free Software Foundation; either version 2.1 of the License, or (at
   6  * your option) any later version.
   7  *
   8  * This library is distributed in the hope that it will be useful, but
   9  * WITHOUT ANY WARRANTY; without even the implied warranty of
  10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser
  11  * General Public License for more details.
  12  *
  13  * You should have received a copy of the GNU Lesser General Public License
  14  * along with this library; if not, write to the Free Software Foundation,
  15  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  16  */
  17
  18 /* Portions based on http://www.dspdimension.com/download smbPitchShift.cpp:
  19  *
  20  * COPYRIGHT 1999-2006 Stephan M. Bernsee <smb [AT] dspdimension [DOT] com>
  21  *
  22  *             The Wide Open License (WOL)
  23  *
  24  * Permission to use, copy, modify, distribute and sell this software and its
  25  * documentation for any purpose is hereby granted without fee, provided that
  26  * the above copyright notice and this license appear in all source copies.
  27  * THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT EXPRESS OR IMPLIED WARRANTY OF
  28  * ANY KIND. See http://www.dspguru.com/wol.htm for more information.
  29  */
  30
  31 #ifdef NDEBUG /* Enable assert always. */
  32 #undef NDEBUG /* Must undef above assert.h or other that might include it. */
  33 #endif
  34
  35 #include "sox_i.h"
  36 #include <assert.h>
  37
  38 #define MAX_FRAME_LENGTH 8192
  39
  40 typedef struct {
  41   unsigned nbends;       /* Number of bends requested */
  42   struct {
  43     char *str;           /* Command-line argument to parse for this bend */
  44     uint64_t start;      /* Start bending when in_pos equals this */
  45     double cents;
  46     uint64_t duration;   /* Number of samples to bend */
  47   } *bends;
  48
  49   unsigned frame_rate;
  50   size_t in_pos;         /* Number of samples read from the input stream */
  51   unsigned bends_pos;    /* Number of bends completed so far */
  52
  53   double shift;
  54
  55   float gInFIFO[MAX_FRAME_LENGTH];
  56   float gOutFIFO[MAX_FRAME_LENGTH];
  57   double gFFTworksp[2 * MAX_FRAME_LENGTH];
  58   float gLastPhase[MAX_FRAME_LENGTH / 2 + 1];
  59   float gSumPhase[MAX_FRAME_LENGTH / 2 + 1];
  60   float gOutputAccum[2 * MAX_FRAME_LENGTH];
  61   float gAnaFreq[MAX_FRAME_LENGTH];
  62   float gAnaMagn[MAX_FRAME_LENGTH];
  63   float gSynFreq[MAX_FRAME_LENGTH];
  64   float gSynMagn[MAX_FRAME_LENGTH];
  65   long gRover;
  66   int fftFrameSize, ovsamp;
  67 } priv_t;
  68
  69 static int parse(sox_effect_t * effp, char **argv, sox_rate_t rate)
  70 {
  71   priv_t *p = (priv_t *) effp->priv;
  72   size_t i;
  73   char const *next;
  74   uint64_t last_seen = 0;
  75   const uint64_t in_length = argv ? 0 :
  76     (effp->in_signal.length != SOX_UNKNOWN_LEN ?
  77      effp->in_signal.length / effp->in_signal.channels : SOX_UNKNOWN_LEN);
  78
  79   for (i = 0; i < p->nbends; ++i) {
  80     if (argv)   /* 1st parse only */
  81       p->bends[i].str = lsx_strdup(argv[i]);
  82
  83     next = lsx_parseposition(rate, p->bends[i].str,
  84              argv ? NULL : &p->bends[i].start, last_seen, in_length, '+');
  85     last_seen = p->bends[i].start;
  86     if (next == NULL || *next != ',')
  87       break;
  88
  89     p->bends[i].cents = strtod(next + 1, (char **)&next);
  90     if (p->bends[i].cents == 0 || *next != ',')
  91       break;
  92
  93     next = lsx_parseposition(rate, next + 1,
  94              argv ? NULL : &p->bends[i].duration, last_seen, in_length, '+');
  95     last_seen = p->bends[i].duration;
  96     if (next == NULL || *next != '\0')
  97       break;
  98
  99     /* sanity checks */
 100     if (!argv && p->bends[i].duration < p->bends[i].start) {
 101       lsx_fail("Bend %" PRIuPTR " has negative width", i+1);
 102       break;
 103     }
 104     if (!argv && i && p->bends[i].start < p->bends[i-1].start) {
 105       lsx_fail("Bend %" PRIuPTR " overlaps with previous one", i+1);
 106       break;
 107     }
 108
 109     p->bends[i].duration -= p->bends[i].start;
 110   }
 111   if (i < p->nbends)
 112     return lsx_usage(effp);
 113   return SOX_SUCCESS;
 114 }
 115
 116 static int create(sox_effect_t * effp, int argc, char **argv)
 117 {
 118   priv_t *p = (priv_t *) effp->priv;
 119   char const * opts = "f:o:";
 120   int c;
 121   lsx_getopt_t optstate;
 122   lsx_getopt_init(argc, argv, opts, NULL, lsx_getopt_flag_none, 1, &optstate);
 123
 124   p->frame_rate = 25;
 125   p->ovsamp = 16;
 126   while ((c = lsx_getopt(&optstate)) != -1) switch (c) {
 127     GETOPT_NUMERIC(optstate, 'f', frame_rate, 10 , 80)
 128     GETOPT_NUMERIC(optstate, 'o', ovsamp,  4 , 32)
 129     default: lsx_fail("unknown option `-%c'", optstate.opt); return lsx_usage(effp);
 130   }
 131   argc -= optstate.ind, argv += optstate.ind;
 132
 133   p->nbends = argc;
 134   p->bends = lsx_calloc(p->nbends, sizeof(*p->bends));
 135   return parse(effp, argv, 0.);     /* No rate yet; parse with dummy */
 136 }
 137
 138 static int start(sox_effect_t * effp)
 139 {
 140   priv_t *p = (priv_t *) effp->priv;
 141   unsigned i;
 142
 143   int n = effp->in_signal.rate / p->frame_rate + .5;
 144   for (p->fftFrameSize = 2; n > 2; p->fftFrameSize <<= 1, n >>= 1);
 145   assert(p->fftFrameSize <= MAX_FRAME_LENGTH);
 146   p->shift = 1;
 147   parse(effp, 0, effp->in_signal.rate); /* Re-parse now rate is known */
 148   p->in_pos = p->bends_pos = 0;
 149   for (i = 0; i < p->nbends; ++i)
 150     if (p->bends[i].duration)
 151       return SOX_SUCCESS;
 152   return SOX_EFF_NULL;
 153 }
 154
 155 static int flow(sox_effect_t * effp, const sox_sample_t * ibuf,
 156                 sox_sample_t * obuf, size_t * isamp, size_t * osamp)
 157 {
 158   priv_t *p = (priv_t *) effp->priv;
 159   size_t i, len = *isamp = *osamp = min(*isamp, *osamp);
 160   double magn, phase, tmp, window, real, imag;
 161   double freqPerBin, expct;
 162   long k, qpd, index, inFifoLatency, stepSize, fftFrameSize2;
 163   float pitchShift = p->shift;
 164
 165   /* set up some handy variables */
 166   fftFrameSize2 = p->fftFrameSize / 2;
 167   stepSize = p->fftFrameSize / p->ovsamp;
 168   freqPerBin = effp->in_signal.rate / p->fftFrameSize;
 169   expct = 2. * M_PI * (double) stepSize / (double) p->fftFrameSize;
 170   inFifoLatency = p->fftFrameSize - stepSize;
 171   if (!p->gRover)
 172     p->gRover = inFifoLatency;
 173
 174   /* main processing loop */
 175   for (i = 0; i < len; i++) {
 176     SOX_SAMPLE_LOCALS;
 177     ++p->in_pos;
 178
 179     /* As long as we have not yet collected enough data just read in */
 180     p->gInFIFO[p->gRover] = SOX_SAMPLE_TO_FLOAT_32BIT(ibuf[i], effp->clips);
 181     obuf[i] = SOX_FLOAT_32BIT_TO_SAMPLE(
 182         p->gOutFIFO[p->gRover - inFifoLatency], effp->clips);
 183     p->gRover++;
 184
 185     /* now we have enough data for processing */
 186     if (p->gRover >= p->fftFrameSize) {
 187       if (p->bends_pos != p->nbends && p->in_pos >=
 188           p->bends[p->bends_pos].start + p->bends[p->bends_pos].duration) {
 189         pitchShift = p->shift *= pow(2., p->bends[p->bends_pos].cents / 1200);
 190         ++p->bends_pos;
 191       }
 192       if (p->bends_pos != p->nbends && p->in_pos >= p->bends[p->bends_pos].start) {
 193         double progress = (double)(p->in_pos - p->bends[p->bends_pos].start) /
 194              p->bends[p->bends_pos].duration;
 195         progress = 1 - cos(M_PI * progress);
 196         progress *= p->bends[p->bends_pos].cents * (.5 / 1200);
 197         pitchShift = p->shift * pow(2., progress);
 198       }
 199
 200       p->gRover = inFifoLatency;
 201
 202       /* do windowing and re,im interleave */
 203       for (k = 0; k < p->fftFrameSize; k++) {
 204         window = -.5 * cos(2 * M_PI * k / (double) p->fftFrameSize) + .5;
 205         p->gFFTworksp[2 * k] = p->gInFIFO[k] * window;
 206         p->gFFTworksp[2 * k + 1] = 0.;
 207       }
 208
 209       /* ***************** ANALYSIS ******************* */
 210       lsx_safe_cdft(2 * p->fftFrameSize, 1, p->gFFTworksp);
 211
 212       /* this is the analysis step */
 213       for (k = 0; k <= fftFrameSize2; k++) {
 214         /* de-interlace FFT buffer */
 215         real = p->gFFTworksp[2 * k];
 216         imag = - p->gFFTworksp[2 * k + 1];
 217
 218         /* compute magnitude and phase */
 219         magn = 2. * sqrt(real * real + imag * imag);
 220         phase = atan2(imag, real);
 221
 222         /* compute phase difference */
 223         tmp = phase - p->gLastPhase[k];
 224         p->gLastPhase[k] = phase;
 225
 226         tmp -= (double) k *expct; /* subtract expected phase difference */
 227
 228         /* map delta phase into +/- Pi interval */
 229         qpd = tmp / M_PI;
 230         if (qpd >= 0)
 231           qpd += qpd & 1;
 232         else qpd -= qpd & 1;
 233         tmp -= M_PI * (double) qpd;
 234
 235         /* get deviation from bin frequency from the +/- Pi interval */
 236         tmp = p->ovsamp * tmp / (2. * M_PI);
 237
 238         /* compute the k-th partials' true frequency */
 239         tmp = (double) k *freqPerBin + tmp * freqPerBin;
 240
 241         /* store magnitude and true frequency in analysis arrays */
 242         p->gAnaMagn[k] = magn;
 243         p->gAnaFreq[k] = tmp;
 244
 245       }
 246
 247       /* this does the actual pitch shifting */
 248       memset(p->gSynMagn, 0, p->fftFrameSize * sizeof(float));
 249       memset(p->gSynFreq, 0, p->fftFrameSize * sizeof(float));
 250       for (k = 0; k <= fftFrameSize2; k++) {
 251         index = k * pitchShift;
 252         if (index <= fftFrameSize2) {
 253           p->gSynMagn[index] += p->gAnaMagn[k];
 254           p->gSynFreq[index] = p->gAnaFreq[k] * pitchShift;
 255         }
 256       }
 257
 258       for (k = 0; k <= fftFrameSize2; k++) { /* SYNTHESIS */
 259         /* get magnitude and true frequency from synthesis arrays */
 260         magn = p->gSynMagn[k], tmp = p->gSynFreq[k];
 261         tmp -= (double) k *freqPerBin; /* subtract bin mid frequency */
 262         tmp /= freqPerBin; /* get bin deviation from freq deviation */
 263         tmp = 2. * M_PI * tmp / p->ovsamp; /* take p->ovsamp into account */
 264         tmp += (double) k *expct; /* add the overlap phase advance back in */
 265         p->gSumPhase[k] += tmp; /* accumulate delta phase to get bin phase */
 266         phase = p->gSumPhase[k];
 267         /* get real and imag part and re-interleave */
 268         p->gFFTworksp[2 * k] = magn * cos(phase);
 269         p->gFFTworksp[2 * k + 1] = - magn * sin(phase);
 270       }
 271
 272       for (k = p->fftFrameSize + 2; k < 2 * p->fftFrameSize; k++)
 273         p->gFFTworksp[k] = 0.; /* zero negative frequencies */
 274
 275       lsx_safe_cdft(2 * p->fftFrameSize, -1, p->gFFTworksp);
 276
 277       /* do windowing and add to output accumulator */
 278       for (k = 0; k < p->fftFrameSize; k++) {
 279         window =
 280             -.5 * cos(2. * M_PI * (double) k / (double) p->fftFrameSize) + .5;
 281         p->gOutputAccum[k] +=
 282             2. * window * p->gFFTworksp[2 * k] / (fftFrameSize2 * p->ovsamp);
 283       }
 284       for (k = 0; k < stepSize; k++)
 285         p->gOutFIFO[k] = p->gOutputAccum[k];
 286
 287       memmove(p->gOutputAccum, /* shift accumulator */
 288           p->gOutputAccum + stepSize, p->fftFrameSize * sizeof(float));
 289
 290       for (k = 0; k < inFifoLatency; k++) /* move input FIFO */
 291         p->gInFIFO[k] = p->gInFIFO[k + stepSize];
 292     }
 293   }
 294   return SOX_SUCCESS;
 295 }
 296
 297 static int stop(sox_effect_t * effp)
 298 {
 299   priv_t *p = (priv_t *) effp->priv;
 300
 301   if (p->bends_pos != p->nbends)
 302     lsx_warn("Input audio too short; bends not applied: %u",
 303         p->nbends - p->bends_pos);
 304   return SOX_SUCCESS;
 305 }
 306
 307 static int lsx_kill(sox_effect_t * effp)
 308 {
 309   priv_t *p = (priv_t *) effp->priv;
 310   unsigned i;
 311
 312   for (i = 0; i < p->nbends; ++i)
 313     free(p->bends[i].str);
 314   free(p->bends);
 315   return SOX_SUCCESS;
 316 }
 317
 318 sox_effect_handler_t const *lsx_bend_effect_fn(void)
 319 {
 320   static sox_effect_handler_t handler = {
 321     "bend", "[-f frame-rate(25)] [-o over-sample(16)] {start,cents,end}",
 322     0, create, start, flow, 0, stop, lsx_kill, sizeof(priv_t)
 323   };
 324   return &handler;
 325 }