1;;; versions of the Moore-Klingbeil-Trevisani-Edwards phase-vocoder
2
3(provide 'snd-pvoc.scm)
4
5(define make-pvocoder
6  (let ((+documentation+ "(make-pvocoder fftsize overlap interp analyze edit synthesize) makes a new (Scheme-based, not CLM) phase-vocoder generator"))
7
8    (lambda* (fftsize overlap interp analyze edit synthesize)
9      (let ((N (or fftsize 512)))
10	(let ((N2 (floor (/ N 2)))
11	      (D (floor (/ N (or overlap 4)))))
12
13	  ;; basic: fftsize overlap
14	  ;;  everything else via closures (interp in particular)
15	  ;;  pv: counter ("output" here)
16	  ;;      interp
17	  ;;      fftsize ("N"), hop ("D")
18	  ;;      in-counter ("filptr")
19	  ;;      hamming window scaled
20	  ;;      slot for in-coming data ("in-data") (created on first call)
21	  ;;      float-vectors: ampinc amp freqinc phaseinc phase lastphase
22	  ;;      funcs: analysize, edit, resynthesize
23
24	  (list
25	   interp                        ;output
26	   interp                        ;interp
27	   0                             ;filptr
28	   N                             ;N
29	   (let ((window (make-fft-window hamming-window fftsize)))
30	     (float-vector-scale! window (/ 2.0 (* 0.54 fftsize))) ;den = hamming window integrated
31	     window)                     ; window
32	   D                             ;D
33	   #f                            ;in-data (created in pvocoder gen)
34	   (make-float-vector fftsize)            ;ampinc
35	   (make-float-vector fftsize)            ;freqs
36	   (make-float-vector N2)                 ;amps
37	   (make-float-vector N2)                 ;phaseinc
38	   (make-float-vector N2)                 ;phases
39	   (make-float-vector N2)                 ;lastphaseinc
40	   analyze
41	   edit
42	   synthesize))))))
43
44;;; pvocoder generator:
45;;     input data func
46;;     analysis func with fallback
47;;     editing func with fallback
48;;     resynthesis func with fallback
49
50(define pvocoder
51  (let ((+documentation+ "(pvocoder pv input) is the phase-vocoder generator associated with make-pvocoder")
52	;; pvocoder list accessors
53	(pvoc-output (lambda (pv) (pv 0)))
54	(set-pvoc-output (lambda (pv val) (set! (pv 0) val)))
55	(pvoc-interp (lambda (pv) (pv 1)))
56	(pvoc-filptr (lambda (pv) (pv 2)))
57	(set-pvoc-filptr (lambda (pv val) (set! (pv 2) val)))
58	(pvoc-N (lambda (pv) (pv 3)))
59	(pvoc-window (lambda (pv) (pv 4)))
60	(pvoc-D (lambda (pv) (pv 5)))
61	(pvoc-in-data (lambda (pv) (pv 6)))
62	(set-pvoc-in-data (lambda (pv val) (set! (pv 6) val)))
63	(pvoc-ampinc (lambda (pv) (pv 7)))
64	(pvoc-freqs (lambda (pv) (pv 8)))
65	(pvoc-amps (lambda (pv) (pv 9)))
66	(pvoc-phaseinc (lambda (pv) (pv 10)))
67	(pvoc-phases (lambda (pv) (pv 11)))
68	(pvoc-lastphase (lambda (pv) (pv 12)))
69	(pvoc-analyze (lambda (pv) (pv 13)))
70	(pvoc-edit (lambda (pv) (pv 14)))
71	(pvoc-synthesize (lambda (pv) (pv 15)))
72	(sine-bank (lambda* (amps phases size)
73		     (let ((len (or size (length amps)))
74			   (sum 0.0))
75		       (do ((i 0 (+ i 1)))
76			   ((= i len))
77			 (set! sum (+ sum (* (amps i) (sin (phases i))))))
78		       sum)))
79	(pi2 (* 2.0 pi)))
80
81    (lambda (pv input)
82      (when (>= (pvoc-output pv) (pvoc-interp pv))
83	;; get next block of amp/phase info
84	(let ((N (pvoc-N pv))
85	      (D (pvoc-D pv))
86	      (amps (pvoc-ampinc pv))
87	      (freqs (pvoc-freqs pv))
88	      (filptr (pvoc-filptr pv)))
89
90	  (if (pvoc-analyze pv)
91	      ((pvoc-analyze pv) pv input)
92	      ;; if no analysis func:
93	      (begin
94		(fill! freqs 0.0)
95		(set-pvoc-output pv 0)
96		(if (not (pvoc-in-data pv))
97		    (begin
98		      (set-pvoc-in-data pv (make-float-vector N))
99		      (do ((i 0 (+ i 1)))
100			  ((= i N))
101			(set! ((pvoc-in-data pv) i) (input))))
102		    (let ((indat (pvoc-in-data pv)))
103		      ;; extra loop here since I find the optimized case confusing (we could dispense with the data move)
104		      (float-vector-move! indat 0 D)
105		      (do ((i (- N D) (+ i 1)))
106			  ((= i N))
107			(set! (indat i) (input)))))
108		(let ((buf (modulo filptr N)))
109		  (if (= buf 0)
110		      (begin
111			(fill! amps 0.0)
112			(float-vector-add! amps (pvoc-in-data pv))
113			(float-vector-multiply! amps (pvoc-window pv)))
114		      (do ((k 0 (+ k 1)))
115			  ((= k N))
116			(set! (amps buf) (* ((pvoc-window pv) k) ((pvoc-in-data pv) k)))
117			(set! buf (+ 1 buf))
118			(if (= buf N) (set! buf 0)))))
119		(set-pvoc-filptr pv (+ filptr D))
120		(mus-fft amps freqs N 1)
121		(rectangular->polar amps freqs)))
122
123	  (if (pvoc-edit pv)
124	      ((pvoc-edit pv) pv)
125	      (let ((lp (pvoc-lastphase pv))
126		    (pscl (/ 1.0 D))
127		    (kscl (/ pi2 N))
128		    (lim (floor (/ N 2))))
129		;; if no editing func:
130		(do ((k 0 (+ k 1)))
131		    ((= k lim))
132		  (let ((phasediff (remainder (- (freqs k) (lp k)) pi2)))
133		    (float-vector-set! lp k (freqs k))
134		    (if (> phasediff pi) (set! phasediff (- phasediff pi2))
135			(if (< phasediff (- pi)) (set! phasediff (+ phasediff pi2))))
136		    (set! (freqs k) (+ (* pscl phasediff) (* k kscl)))))))
137
138	  (let ((scl (/ 1.0 (pvoc-interp pv))))
139	    (float-vector-subtract! amps (pvoc-amps pv))
140	    (float-vector-subtract! freqs (pvoc-phaseinc pv))
141	    (float-vector-scale! amps scl)
142	    (float-vector-scale! freqs scl)
143	    )))
144
145      (set-pvoc-output pv (+ 1 (pvoc-output pv)))
146
147      (if (pvoc-synthesize pv)
148	  ((pvoc-synthesize pv) pv)
149	  ;; if no synthesis func:
150	  ;; synthesize next sample
151	  (begin
152	    (float-vector-add! (pvoc-amps pv) (pvoc-ampinc pv))
153	    (float-vector-add! (pvoc-phaseinc pv) (pvoc-freqs pv))
154	    (float-vector-add! (pvoc-phases pv) (pvoc-phaseinc pv))
155	    (sine-bank (pvoc-amps pv) (pvoc-phases pv))))
156      )))
157
158#|
159(let* ((ind (open-sound "oboe.snd"))
160       (pv (make-pvocoder 256 4 64))
161       (rd (make-sampler 0)))
162  (map-channel (lambda (y) (pvocoder pv (lambda () (rd))))))
163|#
164
165#|
166;;; ---------------- same thing using phase-vocoder gen
167
168(define test-pv-1
169  (lambda (freq)
170    (let* ((reader (make-sampler 0))
171	   (pv (make-phase-vocoder (lambda (dir) (next-sample reader))
172				   512 4 128 1.0
173				   #f ;no change to analysis
174				   #f ;no change to edits
175				   #f ;no change to synthesis
176				   )))
177      (map-channel (lambda (val) (phase-vocoder pv))))))
178
179(define test-pv-2
180  (lambda (freq)
181    (let* ((reader (make-sampler 0))
182	   (pv (make-phase-vocoder (lambda (dir) (next-sample reader))
183				   512 4 128 freq
184				   #f ;no change to analysis
185				   #f
186				   #f ; no change to synthesis
187				   )))
188      (map-channel (lambda (val) (phase-vocoder pv))))))
189
190(define test-pv-3
191  (lambda (time)
192    (let* ((reader (make-sampler 0))
193	   (pv (make-phase-vocoder (lambda (dir) (next-sample reader))
194				   512 4 (floor (* 128 time)) 1.0
195				   #f ;no change to analysis
196				   #f ;no change to edits
197				   #f ;no change to synthesis
198				   ))
199	   (len (floor (* time (framples))))
200	   (data (make-float-vector len)))
201      (do ((i 0 (+ i 1)))
202	  ((= i len))
203	(set! (data i) (phase-vocoder pv)))
204      (free-sampler reader)
205      (float-vector->channel data 0 len))))
206
207(define test-pv-4
208  (lambda (gate)
209    (let* ((reader (make-sampler 0))
210	   (pv (make-phase-vocoder (lambda (dir) (next-sample reader))
211				   512 4 128 1.0
212				   #f ;no change to analysis
213				   (lambda (v)
214				     (let ((N (mus-length v)))
215				       (do ((i 0 (+ i 1)))
216					   ((= i N))
217					 (if (< ((phase-vocoder-amp-increments v) i) gate)
218					     (float-vector-set! (phase-vocoder-amp-increments v) i 0.0)))
219				       #t))
220				   #f ;no change to synthesis
221				   )))
222      (map-channel (lambda (val) (phase-vocoder pv))))))
223|#
224
225
226;;; -------- another version of the phase vocoder --------
227
228(define pvoc
229  (let ((+documentation+ "(pvoc fftsize overlap time pitch gate hoffset snd chn) applies the phase vocoder
230  algorithm to the current sound (i.e. fft analysis, oscil bank resynthesis). 'pitch'
231  specifies the pitch transposition ratio, 'time' - specifies the time dilation ratio,
232  'gate' specifies a resynthesis gate in dB (partials with amplitudes lower than
233  the gate value will not be synthesized), 'hoffset is a pitch offset in Hz.")
234	(pi2 (* 2.0 pi)))
235
236    (lambda* ((fftsize 512) (overlap 4) (time 1.0)
237	      (pitch 1.0) (gate 0.0) (hoffset 0.0)
238	      (snd 0) (chn 0))
239      (let* ((len (framples))
240	     (N2 (floor (/ fftsize 2)))
241	     (window (make-fft-window hamming-window fftsize))
242	     (in-data (channel->float-vector 0 (* fftsize 2) snd chn))
243	     (lastamp (make-float-vector N2))
244	     (lastfreq (make-float-vector N2))
245	     (outlen (floor (* time len)))
246	     (interp (* (floor (/ fftsize overlap)) time))
247	     (obank (make-oscil-bank lastfreq (make-float-vector N2) lastamp))
248	     (filptr 0)
249	     (D (floor (/ fftsize overlap)))
250	     (syngate (if (= 0.0 gate)    ; take a resynthesis gate specificed in dB, convert to linear amplitude
251			  0.0000
252			  (expt 10 (/ (- (abs gate)) 20))))
253	     (poffset (hz->radians hoffset))
254	     (fdr (make-float-vector fftsize))
255	     (fdi (make-float-vector fftsize))
256	     (lastphase (make-float-vector N2))
257	     (ampinc (make-float-vector N2))
258	     (freqinc (make-float-vector N2))
259	     (fundamental (/ pi2 fftsize))
260	     (output interp)
261	     (out-data (make-float-vector (max len outlen)))
262	     (in-data-beg 0))
263
264	(set! window (float-vector-scale! window (/ 2.0 (* 0.54 fftsize)))) ;den = hamming window integrated
265
266	(do ((i 0 (+ i 1)))
267	    ((>= i outlen))
268	  (when (>= output interp) ;; if all the samples have been output then do the next frame
269	    (let ((buffix (modulo filptr fftsize)))
270					; buffix is the index into the input buffer
271					; it wraps around circularly as time increases in the input
272	      (set! output 0)       ; reset the output sample counter
273	      ;; save the old amplitudes and frequencies
274	      (fill! lastamp 0.0)
275	      (fill! lastfreq 0.0)
276	      (float-vector-add! lastamp fdr)
277	      (float-vector-add! lastfreq fdi)
278	      (do ((k 0 (+ k 1)))
279		  ((= k fftsize))
280		;; apply the window and then stuff into the input array
281		(set! (fdr buffix) (* (window k) (in-data (- filptr in-data-beg))))
282		(set! filptr (+ 1 filptr))
283		;; increment the buffer index with wrap around
284		(set! buffix (+ 1 buffix))
285		(if (>= buffix fftsize) (set! buffix 0)))
286	      ;; rewind the file for the next hop
287	      (set! filptr (- (+ filptr D) fftsize))
288	      (if (> filptr (+ in-data-beg fftsize))
289		  (begin
290		    (set! in-data-beg filptr)
291		    (set! in-data (channel->float-vector in-data-beg (* fftsize 2) snd chn))))
292	      ;; no imaginary component input so zero out fdi
293	      (fill! fdi 0.0)
294	      ;; compute the fft
295	      (mus-fft fdr fdi fftsize 1)
296	      ;; now convert into magnitude and interpolated frequency
297	      (do ((k 0 (+ k 1)))
298		  ((= k N2))
299		(let ((a (fdr k))
300		      (b (fdi k))
301		      (phasediff 0))
302		  (let ((mag (sqrt (+ (* a a) (* b b)))))
303		    (set! (fdr k) mag)    ;; current amp stored in fdr
304		    ;; mag is always positive
305		    ;; if it is zero then the phase difference is zero
306		    (if (> mag 0)
307			(let ((phase (- (atan b a))))
308			  (set! phasediff (- phase (lastphase k)))
309			  (set! (lastphase k) phase)
310			  ;; frequency wrapping from Moore p. 254
311			  (if (> phasediff pi) (do () ((<= phasediff pi)) (set! phasediff (- phasediff pi2))))
312			  (if (< phasediff (- pi)) (do () ((>= phasediff (- pi))) (set! phasediff (+ phasediff pi2)))))))
313		  ;; current frequency stored in fdi
314		  ;; scale by the pitch transposition
315		  (set! (fdi k) (* pitch (+ (/ phasediff D) (* k fundamental) poffset))))
316		;; resynthesis gating
317		(if (< (fdr k) syngate) (set! (fdr k) 0.0))
318		;; take (lastamp k) and count up to (fdr k)
319		;; interpolating by ampinc
320		(set! (ampinc k) (/ (- (fdr k) (lastamp k)) interp))
321		;; take (lastfreq k) and count up to (fdi k)
322		;; interpolating by freqinc
323		(set! (freqinc k) (/ (- (fdi k) (lastfreq k)) interp)))))
324	  ;; loop over the partials interpolate frequency and amplitude
325	  (float-vector-add! lastamp ampinc)
326	  (float-vector-add! lastfreq freqinc)
327	  (float-vector-set! out-data i (oscil-bank obank))
328	  (set! output (+ 1 output)))
329	(float-vector->channel out-data 0 (max len outlen))))))
330