1 /*
2     paulstretch.c:
3 
4     This is an implementation of the paulstretch algorithm by
5     Paul Nasca Octavian, based off of the Numpy/Scipy python
6     implementation written by the same author.
7 
8     Copyright (C) 2016 by Paul Batchelor
9 
10     This file is part of Csound.
11 
12     The Csound Library is free software; you can redistribute it
13     and/or modify it under the terms of the GNU Lesser General Public
14     License as published by the Free Software Foundation; either
15     version 2.1 of the License, or (at your option) any later version.
16 
17     Csound is distributed in the hope that it will be useful,
18     but WITHOUT ANY WARRANTY; without even the implied warranty of
19     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
20     GNU Lesser General Public License for more details.
21 
22     You should have received a copy of the GNU Lesser General Public
23     License along with Csound; if not, write to the Free Software
24     Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
25     02110-1301 USA
26 */
27 #include <stdlib.h>
28 #include <complex.h>
29 
30 #ifdef _MSC_VER
31 #define _USE_MATH_DEFINES
32 #endif
33 #include <math.h>
34 #include <math.h>
35 #include "csoundCore.h"
36 #include "interlocks.h"
37 #include "H/fftlib.h"
38 
39 #ifdef ANDROID
40 float crealf(_Complex float);
41 float cimagf(_Complex float);
42 #endif
43 
44 typedef struct {
45     OPDS h;
46     MYFLT *out, *stretch, *winsize, *ifn;
47     MYFLT start_pos, displace_pos;
48     MYFLT *window;
49     MYFLT *old_windowed_buf;
50     MYFLT *hinv_buf;
51     MYFLT *output;
52     FUNC *ft;
53     uint32_t windowsize;
54     uint32_t half_windowsize;
55     MYFLT *tmp;
56     uint32_t counter;
57     AUXCH m_window;
58     AUXCH m_old_windowed_buf;
59     AUXCH m_hinv_buf;
60     AUXCH m_output;
61     AUXCH m_tmp;
62 } PAULSTRETCH;
63 
compute_block(CSOUND * csound,PAULSTRETCH * p)64 static void compute_block(CSOUND *csound, PAULSTRETCH *p)
65 {
66     uint32_t istart_pos = floor(p->start_pos);
67     uint32_t pos;
68     uint32_t i;
69     uint32_t windowsize = p->windowsize;
70     uint32_t half_windowsize = p->half_windowsize;
71     MYFLT *hinv_buf = p->hinv_buf;
72     MYFLT *old_windowed_buf= p->old_windowed_buf;
73     MYFLT *tbl = p->ft->ftable;
74     MYFLT *window = p->window;
75     MYFLT *output= p->output;
76     MYFLT *tmp = p->tmp;
77     for (i = 0; i < windowsize; i++) {
78       pos = istart_pos + i;
79       if (LIKELY(pos < p->ft->flen)) {
80         tmp[i] = tbl[pos] * window[i];
81       } else {
82         tmp[i] = FL(0.0);
83       }
84     }
85     /* re-order bins and take FFT */
86     tmp[p->windowsize] = tmp[1];
87     tmp[p->windowsize + 1] = FL(0.0);
88     csoundRealFFTnp2(csound, tmp, p->windowsize);
89 
90     /* randomize phase */
91     for (i = 0; i < windowsize + 2; i += 2) {
92       MYFLT mag = HYPOT(tmp[i], tmp[i + 1]);
93       // Android 5.1 does not seem to have cexpf ...
94       // complex ph = cexpf(I * ((MYFLT)rand() / RAND_MAX) * 2 * PI);
95       // so ...
96       MYFLT  x = (((MYFLT)rand() / RAND_MAX) * 2 * PI);
97 #ifdef MSVC
98       // TODO - Double check this is equivalent to non-windows complex definition
99           _Fcomplex ph = { cos(x), sin(x) };
100 #else
101       complex double ph =  cos(x) + I*sin(x);
102 #endif
103       tmp[i] = mag * (MYFLT)crealf(ph);
104       tmp[i + 1] = mag * (MYFLT)cimagf(ph);
105     }
106 
107     /* re-order bins and take inverse FFT */
108     tmp[1] = tmp[p->windowsize];
109     csoundInverseRealFFTnp2(csound, tmp, p->windowsize);
110 
111     /* apply window and overlap */
112     for (i = 0; i < windowsize; i++) {
113       tmp[i] *= window[i];
114       if (i < half_windowsize) {
115         output[i] = (MYFLT)(tmp[i] + old_windowed_buf[half_windowsize + i]);
116         output[i] *= hinv_buf[i];
117       }
118       old_windowed_buf[i] = tmp[i];
119     }
120     p->start_pos += p->displace_pos;
121 }
122 
ps_init(CSOUND * csound,PAULSTRETCH * p)123 static int32_t ps_init(CSOUND* csound, PAULSTRETCH *p)
124 {
125     FUNC *ftp = csound->FTnp2Find(csound, p->ifn);
126     uint32_t i = 0;
127     uint32_t size;
128 
129     if (ftp == NULL)
130       return csound->InitError(csound, Str("paulstretch: table not found"));
131 
132     p->ft = ftp;
133     p->windowsize = (uint32_t)FLOOR((CS_ESR * *p->winsize));
134     if (p->windowsize < 16) {
135       p->windowsize = 16;
136     }
137     p->half_windowsize = p->windowsize / 2;
138     p->displace_pos = (p->windowsize * FL(0.5)) / *p->stretch;
139 
140     size = sizeof(MYFLT) * p->windowsize;
141     csound->AuxAlloc(csound, size, &(p->m_window));
142     p->window = p->m_window.auxp;
143 
144     csound->AuxAlloc(csound, size, &p->m_old_windowed_buf);
145     p->old_windowed_buf = p->m_old_windowed_buf.auxp;
146 
147     csound->AuxAlloc(csound,
148                      (size_t)(sizeof(MYFLT) * p->half_windowsize), &p->m_hinv_buf);
149     p->hinv_buf = p->m_hinv_buf.auxp;
150 
151     csound->AuxAlloc(csound,
152                      (size_t)(sizeof(MYFLT) * p->half_windowsize), &p->m_output);
153     p->output = p->m_output.auxp;
154 
155     csound->AuxAlloc(csound, size + 2 * sizeof(MYFLT), &p->m_tmp);
156     p->tmp = p->m_tmp.auxp;
157 
158     /* Create Hann window */
159     for (i = 0; i < p->windowsize; i++) {
160       p->window[i] = FL(0.5) - COS(i * TWOPI_F / (p->windowsize - 1)) * FL(0.5);
161     }
162     /* create inverse Hann window */
163     {
164       MYFLT hinv_sqrt2 = (1 + SQRT(FL(0.5))) * FL(0.5);
165       for (i = 0; i < p->half_windowsize; i++) {
166         p->hinv_buf[i] =
167           hinv_sqrt2 - (FL(1.0) - hinv_sqrt2) *
168           COS(i * TWOPI_F / p->half_windowsize);
169       }
170     }
171     p->start_pos = FL(0.0);
172     p->counter = 0;
173 
174     return OK;
175 }
176 
paulstretch_perf(CSOUND * csound,PAULSTRETCH * p)177 static int32_t paulstretch_perf(CSOUND* csound, PAULSTRETCH *p)
178 {
179     uint32_t offset = p->h.insdshead->ksmps_offset;
180     uint32_t early  = p->h.insdshead->ksmps_no_end;
181     uint32_t n, nsmps = CS_KSMPS;
182     MYFLT *out = p->out;
183 
184     if (UNLIKELY(offset)) {
185       memset(p->out, '\0', offset*sizeof(MYFLT));
186       memset(p->out, '\0', offset*sizeof(MYFLT));
187     }
188     if (UNLIKELY(early)) {
189       nsmps -= early;
190       memset(&p->out[nsmps], '\0', early*sizeof(MYFLT));
191     }
192 
193     for (n = offset; n < nsmps; n++) {
194       if (p->counter == 0) {
195         compute_block(csound, p);
196       }
197       out[n] = p->output[p->counter];
198       p->counter = (p->counter + 1) % p->half_windowsize;
199     }
200     return OK;
201 }
202 
203 static OENTRY paulstretch_localops[] = {
204   { "paulstretch", (int32_t) sizeof(PAULSTRETCH), TR, 3, "a", "iii",
205     (int32_t (*)(CSOUND *, void *)) ps_init,
206     (int32_t (*)(CSOUND *, void *)) paulstretch_perf}
207 };
208 
209 LINKAGE_BUILTIN(paulstretch_localops)
210