Mercurial > hg > sonic-visualiser
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 } |