comparison audioio/PhaseVocoderTimeStretcher.cpp @ 14:085f34c73939

* IntegerTimeStretcher -> PhaseVocoderTimeStretcher (no longer confined to integer multiples)
author Chris Cannam
date Wed, 13 Sep 2006 11:06:28 +0000
parents audioio/IntegerTimeStretcher.cpp@00ed645f4175
children cc566264c935
comparison
equal deleted inserted replaced
13:00ed645f4175 14:085f34c73939
1 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */
2
3 /*
4 Sonic Visualiser
5 An audio file viewer and annotation editor.
6 Centre for Digital Music, Queen Mary, University of London.
7 This file copyright 2006 Chris Cannam.
8
9 This program is free software; you can redistribute it and/or
10 modify it under the terms of the GNU General Public License as
11 published by the Free Software Foundation; either version 2 of the
12 License, or (at your option) any later version. See the file
13 COPYING included with this distribution for more information.
14 */
15
16 #include "PhaseVocoderTimeStretcher.h"
17
18 #include <iostream>
19 #include <cassert>
20
21 //#define DEBUG_PHASE_VOCODER_TIME_STRETCHER 1
22
23 PhaseVocoderTimeStretcher::PhaseVocoderTimeStretcher(float ratio,
24 size_t maxProcessInputBlockSize,
25 size_t inputIncrement,
26 size_t windowSize,
27 WindowType windowType) :
28 m_ratio(ratio),
29 m_n1(inputIncrement),
30 m_n2(lrintf(m_n1 * ratio)),
31 m_wlen(std::max(windowSize, m_n2 * 2)),
32 m_inbuf(m_wlen),
33 m_outbuf(maxProcessInputBlockSize * ratio + 1024) //!!!
34 {
35 m_window = new Window<float>(windowType, m_wlen),
36
37 m_time = (fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex) * m_wlen);
38 m_freq = (fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex) * m_wlen);
39 m_dbuf = (float *)fftwf_malloc(sizeof(float) * m_wlen);
40 m_mashbuf = (float *)fftwf_malloc(sizeof(float) * m_wlen);
41 m_modulationbuf = (float *)fftwf_malloc(sizeof(float) * m_wlen);
42 m_prevPhase = (float *)fftwf_malloc(sizeof(float) * m_wlen);
43 m_prevAdjustedPhase = (float *)fftwf_malloc(sizeof(float) * m_wlen);
44
45 m_plan = fftwf_plan_dft_1d(m_wlen, m_time, m_freq, FFTW_FORWARD, FFTW_ESTIMATE);
46 m_iplan = fftwf_plan_dft_c2r_1d(m_wlen, m_freq, m_dbuf, FFTW_ESTIMATE);
47
48 for (int i = 0; i < m_wlen; ++i) {
49 m_mashbuf[i] = 0.0;
50 m_modulationbuf[i] = 0.0;
51 m_prevPhase[i] = 0.0;
52 m_prevAdjustedPhase[i] = 0.0;
53 }
54 }
55
56 PhaseVocoderTimeStretcher::~PhaseVocoderTimeStretcher()
57 {
58 std::cerr << "PhaseVocoderTimeStretcher::~PhaseVocoderTimeStretcher" << std::endl;
59
60 fftwf_destroy_plan(m_plan);
61 fftwf_destroy_plan(m_iplan);
62
63 fftwf_free(m_time);
64 fftwf_free(m_freq);
65 fftwf_free(m_dbuf);
66 fftwf_free(m_mashbuf);
67 fftwf_free(m_modulationbuf);
68 fftwf_free(m_prevPhase);
69 fftwf_free(m_prevAdjustedPhase);
70
71 delete m_window;
72 }
73
74 size_t
75 PhaseVocoderTimeStretcher::getProcessingLatency() const
76 {
77 return getWindowSize() - getInputIncrement();
78 }
79
80 void
81 PhaseVocoderTimeStretcher::process(float *input, float *output, size_t samples)
82 {
83 // We need to add samples from input to our internal buffer. When
84 // we have m_windowSize samples in the buffer, we can process it,
85 // move the samples back by m_n1 and write the output onto our
86 // internal output buffer. If we have (samples * ratio) samples
87 // in that, we can write m_n2 of them back to output and return
88 // (otherwise we have to write zeroes).
89
90 // When we process, we write m_wlen to our fixed output buffer
91 // (m_mashbuf). We then pull out the first m_n2 samples from that
92 // buffer, push them into the output ring buffer, and shift
93 // m_mashbuf left by that amount.
94
95 // The processing latency is then m_wlen - m_n2.
96
97 size_t consumed = 0;
98
99 #ifdef DEBUG_PHASE_VOCODER_TIME_STRETCHER
100 std::cerr << "PhaseVocoderTimeStretcher::process(" << samples << ", consumed = " << consumed << "), writable " << m_inbuf.getWriteSpace() <<", readable "<< m_outbuf.getReadSpace() << std::endl;
101 #endif
102
103 while (consumed < samples) {
104
105 size_t writable = m_inbuf.getWriteSpace();
106 writable = std::min(writable, samples - consumed);
107
108 if (writable == 0) {
109 //!!! then what? I don't think this should happen, but
110 std::cerr << "WARNING: PhaseVocoderTimeStretcher::process: writable == 0" << std::endl;
111 break;
112 }
113
114 #ifdef DEBUG_PHASE_VOCODER_TIME_STRETCHER
115 std::cerr << "writing " << writable << " from index " << consumed << " to inbuf, consumed will be " << consumed + writable << std::endl;
116 #endif
117 m_inbuf.write(input + consumed, writable);
118 consumed += writable;
119
120 while (m_inbuf.getReadSpace() >= m_wlen &&
121 m_outbuf.getWriteSpace() >= m_n2) {
122
123 // We know we have at least m_wlen samples available
124 // in m_inbuf. We need to peek m_wlen of them for
125 // processing, and then read m_n1 to advance the read
126 // pointer.
127
128 size_t got = m_inbuf.peek(m_dbuf, m_wlen);
129 assert(got == m_wlen);
130
131 processBlock(m_dbuf, m_mashbuf, m_modulationbuf);
132
133 #ifdef DEBUG_PHASE_VOCODER_TIME_STRETCHER
134 std::cerr << "writing first " << m_n2 << " from mashbuf, skipping " << m_n1 << " on inbuf " << std::endl;
135 #endif
136 m_inbuf.skip(m_n1);
137
138 for (size_t i = 0; i < m_n2; ++i) {
139 if (m_modulationbuf[i] > 0.f) {
140 m_mashbuf[i] /= m_modulationbuf[i];
141 }
142 }
143
144 m_outbuf.write(m_mashbuf, m_n2);
145
146 for (size_t i = 0; i < m_wlen - m_n2; ++i) {
147 m_mashbuf[i] = m_mashbuf[i + m_n2];
148 m_modulationbuf[i] = m_modulationbuf[i + m_n2];
149 }
150
151 for (size_t i = m_wlen - m_n2; i < m_wlen; ++i) {
152 m_mashbuf[i] = 0.0f;
153 m_modulationbuf[i] = 0.0f;
154 }
155 }
156
157 // std::cerr << "WARNING: PhaseVocoderTimeStretcher::process: writespace not enough for output increment (" << m_outbuf.getWriteSpace() << " < " << m_n2 << ")" << std::endl;
158 // }
159
160 #ifdef DEBUG_PHASE_VOCODER_TIME_STRETCHER
161 std::cerr << "loop ended: inbuf read space " << m_inbuf.getReadSpace() << ", outbuf write space " << m_outbuf.getWriteSpace() << std::endl;
162 #endif
163 }
164
165 size_t toRead = lrintf(samples * m_ratio);
166
167 if (m_outbuf.getReadSpace() < toRead) {
168 std::cerr << "WARNING: PhaseVocoderTimeStretcher::process: not enough data (yet?) (" << m_outbuf.getReadSpace() << " < " << toRead << ")" << std::endl;
169 size_t fill = toRead - m_outbuf.getReadSpace();
170 for (size_t i = 0; i < fill; ++i) {
171 output[i] = 0.0;
172 }
173 m_outbuf.read(output + fill, m_outbuf.getReadSpace());
174 } else {
175 #ifdef DEBUG_PHASE_VOCODER_TIME_STRETCHER
176 std::cerr << "enough data - writing " << toRead << " from outbuf" << std::endl;
177 #endif
178 m_outbuf.read(output, toRead);
179 }
180
181 #ifdef DEBUG_PHASE_VOCODER_TIME_STRETCHER
182 std::cerr << "PhaseVocoderTimeStretcher::process returning" << std::endl;
183 #endif
184 }
185
186 void
187 PhaseVocoderTimeStretcher::processBlock(float *buf, float *out, float *modulation)
188 {
189 size_t i;
190
191 // buf contains m_wlen samples; out contains enough space for
192 // m_wlen * ratio samples (we mix into out, rather than replacing)
193
194 #ifdef DEBUG_PHASE_VOCODER_TIME_STRETCHER
195 std::cerr << "PhaseVocoderTimeStretcher::processBlock" << std::endl;
196 #endif
197
198 m_window->cut(buf);
199
200 for (i = 0; i < m_wlen/2; ++i) {
201 float temp = buf[i];
202 buf[i] = buf[i + m_wlen/2];
203 buf[i + m_wlen/2] = temp;
204 }
205
206 for (i = 0; i < m_wlen; ++i) {
207 m_time[i][0] = buf[i];
208 m_time[i][1] = 0.0;
209 }
210
211 fftwf_execute(m_plan); // m_time -> m_freq
212
213 for (i = 0; i < m_wlen; ++i) {
214
215 float mag = sqrtf(m_freq[i][0] * m_freq[i][0] +
216 m_freq[i][1] * m_freq[i][1]);
217
218 float phase = princargf(atan2f(m_freq[i][1], m_freq[i][0]));
219
220 float omega = (2 * M_PI * m_n1 * i) / m_wlen;
221
222 float expectedPhase = m_prevPhase[i] + omega;
223
224 float phaseError = princargf(phase - expectedPhase);
225
226 float phaseIncrement = (omega + phaseError) / m_n1;
227
228 float adjustedPhase = m_prevAdjustedPhase[i] + m_n2 * phaseIncrement;
229
230 float real = mag * cosf(adjustedPhase);
231 float imag = mag * sinf(adjustedPhase);
232 m_freq[i][0] = real;
233 m_freq[i][1] = imag;
234
235 m_prevPhase[i] = phase;
236 m_prevAdjustedPhase[i] = adjustedPhase;
237 }
238
239 fftwf_execute(m_iplan); // m_freq -> in, inverse fft
240
241 for (i = 0; i < m_wlen/2; ++i) {
242 float temp = buf[i] / m_wlen;
243 buf[i] = buf[i + m_wlen/2] / m_wlen;
244 buf[i + m_wlen/2] = temp;
245 }
246
247 m_window->cut(buf);
248 /*
249 int div = m_wlen / m_n2;
250 if (div > 1) div /= 2;
251 for (i = 0; i < m_wlen; ++i) {
252 buf[i] /= div;
253 }
254 */
255 for (i = 0; i < m_wlen; ++i) {
256 out[i] += buf[i];
257 modulation[i] += m_window->getValue(i);
258 }
259 }