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