cannam@64
|
1 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */
|
cannam@64
|
2
|
cannam@64
|
3 /*
|
cannam@64
|
4 Vamp
|
cannam@64
|
5
|
cannam@64
|
6 An API for audio analysis and feature extraction plugins.
|
cannam@64
|
7
|
cannam@64
|
8 Centre for Digital Music, Queen Mary, University of London.
|
cannam@64
|
9 Copyright 2006 Chris Cannam.
|
cannam@64
|
10
|
cannam@64
|
11 Permission is hereby granted, free of charge, to any person
|
cannam@64
|
12 obtaining a copy of this software and associated documentation
|
cannam@64
|
13 files (the "Software"), to deal in the Software without
|
cannam@64
|
14 restriction, including without limitation the rights to use, copy,
|
cannam@64
|
15 modify, merge, publish, distribute, sublicense, and/or sell copies
|
cannam@64
|
16 of the Software, and to permit persons to whom the Software is
|
cannam@64
|
17 furnished to do so, subject to the following conditions:
|
cannam@64
|
18
|
cannam@64
|
19 The above copyright notice and this permission notice shall be
|
cannam@64
|
20 included in all copies or substantial portions of the Software.
|
cannam@64
|
21
|
cannam@64
|
22 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
cannam@64
|
23 EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
|
cannam@64
|
24 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
cannam@64
|
25 NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR
|
cannam@64
|
26 ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
|
cannam@64
|
27 CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
|
cannam@64
|
28 WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
|
cannam@64
|
29
|
cannam@64
|
30 Except as contained in this notice, the names of the Centre for
|
cannam@64
|
31 Digital Music; Queen Mary, University of London; and Chris Cannam
|
cannam@64
|
32 shall not be used in advertising or otherwise to promote the sale,
|
cannam@64
|
33 use or other dealings in this Software without prior written
|
cannam@64
|
34 authorization.
|
cannam@64
|
35 */
|
cannam@64
|
36
|
cannam@64
|
37 #include "PluginInputDomainAdapter.h"
|
cannam@64
|
38
|
cannam@64
|
39 #include <cmath>
|
cannam@64
|
40
|
cannam@64
|
41 namespace Vamp {
|
cannam@64
|
42
|
cannam@64
|
43 namespace HostExt {
|
cannam@64
|
44
|
cannam@70
|
45 class PluginInputDomainAdapter::Impl
|
cannam@70
|
46 {
|
cannam@70
|
47 public:
|
cannam@70
|
48 Impl(Plugin *plugin, float inputSampleRate);
|
cannam@70
|
49 ~Impl();
|
cannam@70
|
50
|
cannam@70
|
51 bool initialise(size_t channels, size_t stepSize, size_t blockSize);
|
cannam@70
|
52
|
cannam@70
|
53 size_t getPreferredStepSize() const;
|
cannam@70
|
54 size_t getPreferredBlockSize() const;
|
cannam@70
|
55
|
cannam@70
|
56 FeatureSet process(const float *const *inputBuffers, RealTime timestamp);
|
cannam@70
|
57
|
cannam@70
|
58 protected:
|
cannam@70
|
59 Plugin *m_plugin;
|
cannam@70
|
60 float m_inputSampleRate;
|
cannam@70
|
61 size_t m_channels;
|
cannam@70
|
62 size_t m_blockSize;
|
cannam@70
|
63 float **m_freqbuf;
|
cannam@70
|
64 double *m_ri;
|
cannam@70
|
65 double *m_ro;
|
cannam@70
|
66 double *m_io;
|
cannam@70
|
67
|
cannam@70
|
68 void fft(unsigned int n, bool inverse,
|
cannam@70
|
69 double *ri, double *ii, double *ro, double *io);
|
cannam@70
|
70
|
cannam@70
|
71 size_t makeBlockSizeAcceptable(size_t) const;
|
cannam@70
|
72 };
|
cannam@70
|
73
|
cannam@64
|
74 PluginInputDomainAdapter::PluginInputDomainAdapter(Plugin *plugin) :
|
cannam@70
|
75 PluginWrapper(plugin)
|
cannam@70
|
76 {
|
cannam@70
|
77 m_impl = new Impl(plugin, m_inputSampleRate);
|
cannam@70
|
78 }
|
cannam@70
|
79
|
cannam@70
|
80 PluginInputDomainAdapter::~PluginInputDomainAdapter()
|
cannam@70
|
81 {
|
cannam@70
|
82 delete m_impl;
|
cannam@70
|
83 }
|
cannam@70
|
84
|
cannam@70
|
85 bool
|
cannam@70
|
86 PluginInputDomainAdapter::initialise(size_t channels, size_t stepSize, size_t blockSize)
|
cannam@70
|
87 {
|
cannam@70
|
88 return m_impl->initialise(channels, stepSize, blockSize);
|
cannam@70
|
89 }
|
cannam@70
|
90
|
cannam@70
|
91 Plugin::InputDomain
|
cannam@70
|
92 PluginInputDomainAdapter::getInputDomain() const
|
cannam@70
|
93 {
|
cannam@70
|
94 return TimeDomain;
|
cannam@70
|
95 }
|
cannam@70
|
96
|
cannam@70
|
97 size_t
|
cannam@70
|
98 PluginInputDomainAdapter::getPreferredStepSize() const
|
cannam@70
|
99 {
|
cannam@70
|
100 return m_impl->getPreferredStepSize();
|
cannam@70
|
101 }
|
cannam@70
|
102
|
cannam@70
|
103 size_t
|
cannam@70
|
104 PluginInputDomainAdapter::getPreferredBlockSize() const
|
cannam@70
|
105 {
|
cannam@70
|
106 return m_impl->getPreferredBlockSize();
|
cannam@70
|
107 }
|
cannam@70
|
108
|
cannam@70
|
109 Plugin::FeatureSet
|
cannam@70
|
110 PluginInputDomainAdapter::process(const float *const *inputBuffers, RealTime timestamp)
|
cannam@70
|
111 {
|
cannam@70
|
112 return m_impl->process(inputBuffers, timestamp);
|
cannam@70
|
113 }
|
cannam@70
|
114
|
cannam@70
|
115 PluginInputDomainAdapter::Impl::Impl(Plugin *plugin, float inputSampleRate) :
|
cannam@70
|
116 m_plugin(plugin),
|
cannam@70
|
117 m_inputSampleRate(inputSampleRate),
|
cannam@64
|
118 m_channels(0),
|
cannam@64
|
119 m_blockSize(0),
|
cannam@64
|
120 m_freqbuf(0)
|
cannam@64
|
121 {
|
cannam@64
|
122 }
|
cannam@64
|
123
|
cannam@70
|
124 PluginInputDomainAdapter::Impl::~Impl()
|
cannam@64
|
125 {
|
cannam@70
|
126 // the adapter will delete the plugin
|
cannam@70
|
127
|
cannam@70
|
128 if (m_channels > 0) {
|
cannam@70
|
129 for (size_t c = 0; c < m_channels; ++c) {
|
cannam@70
|
130 delete[] m_freqbuf[c];
|
cannam@70
|
131 }
|
cannam@70
|
132 delete[] m_freqbuf;
|
cannam@70
|
133 delete[] m_ri;
|
cannam@70
|
134 delete[] m_ro;
|
cannam@70
|
135 delete[] m_io;
|
cannam@70
|
136 }
|
cannam@64
|
137 }
|
cannam@64
|
138
|
cannam@64
|
139 bool
|
cannam@70
|
140 PluginInputDomainAdapter::Impl::initialise(size_t channels, size_t stepSize, size_t blockSize)
|
cannam@64
|
141 {
|
cannam@64
|
142 if (m_plugin->getInputDomain() == TimeDomain) {
|
cannam@64
|
143
|
cannam@64
|
144 m_blockSize = blockSize;
|
cannam@64
|
145 m_channels = channels;
|
cannam@64
|
146
|
cannam@64
|
147 return m_plugin->initialise(channels, stepSize, blockSize);
|
cannam@64
|
148 }
|
cannam@64
|
149
|
cannam@64
|
150 if (blockSize < 2) {
|
cannam@70
|
151 std::cerr << "ERROR: Vamp::HostExt::PluginInputDomainAdapter::Impl::initialise: blocksize < 2 not supported" << std::endl;
|
cannam@64
|
152 return false;
|
cannam@64
|
153 }
|
cannam@64
|
154
|
cannam@64
|
155 if (blockSize & (blockSize-1)) {
|
cannam@70
|
156 std::cerr << "ERROR: Vamp::HostExt::PluginInputDomainAdapter::Impl::initialise: non-power-of-two\nblocksize " << blockSize << " not supported" << std::endl;
|
cannam@64
|
157 return false;
|
cannam@64
|
158 }
|
cannam@64
|
159
|
cannam@64
|
160 if (m_channels > 0) {
|
cannam@64
|
161 for (size_t c = 0; c < m_channels; ++c) {
|
cannam@64
|
162 delete[] m_freqbuf[c];
|
cannam@64
|
163 }
|
cannam@64
|
164 delete[] m_freqbuf;
|
cannam@64
|
165 delete[] m_ri;
|
cannam@64
|
166 delete[] m_ro;
|
cannam@64
|
167 delete[] m_io;
|
cannam@64
|
168 }
|
cannam@64
|
169
|
cannam@64
|
170 m_blockSize = blockSize;
|
cannam@64
|
171 m_channels = channels;
|
cannam@64
|
172
|
cannam@64
|
173 m_freqbuf = new float *[m_channels];
|
cannam@64
|
174 for (size_t c = 0; c < m_channels; ++c) {
|
cannam@64
|
175 m_freqbuf[c] = new float[m_blockSize + 2];
|
cannam@64
|
176 }
|
cannam@64
|
177 m_ri = new double[m_blockSize];
|
cannam@64
|
178 m_ro = new double[m_blockSize];
|
cannam@64
|
179 m_io = new double[m_blockSize];
|
cannam@64
|
180
|
cannam@64
|
181 return m_plugin->initialise(channels, stepSize, blockSize);
|
cannam@64
|
182 }
|
cannam@64
|
183
|
cannam@64
|
184 size_t
|
cannam@70
|
185 PluginInputDomainAdapter::Impl::getPreferredStepSize() const
|
cannam@64
|
186 {
|
cannam@64
|
187 size_t step = m_plugin->getPreferredStepSize();
|
cannam@64
|
188
|
cannam@64
|
189 if (step == 0 && (m_plugin->getInputDomain() == FrequencyDomain)) {
|
cannam@64
|
190 step = getPreferredBlockSize() / 2;
|
cannam@64
|
191 }
|
cannam@64
|
192
|
cannam@64
|
193 return step;
|
cannam@64
|
194 }
|
cannam@64
|
195
|
cannam@64
|
196 size_t
|
cannam@70
|
197 PluginInputDomainAdapter::Impl::getPreferredBlockSize() const
|
cannam@64
|
198 {
|
cannam@64
|
199 size_t block = m_plugin->getPreferredBlockSize();
|
cannam@64
|
200
|
cannam@64
|
201 if (m_plugin->getInputDomain() == FrequencyDomain) {
|
cannam@64
|
202 if (block == 0) {
|
cannam@64
|
203 block = 1024;
|
cannam@64
|
204 } else {
|
cannam@64
|
205 block = makeBlockSizeAcceptable(block);
|
cannam@64
|
206 }
|
cannam@64
|
207 }
|
cannam@64
|
208
|
cannam@64
|
209 return block;
|
cannam@64
|
210 }
|
cannam@64
|
211
|
cannam@64
|
212 size_t
|
cannam@70
|
213 PluginInputDomainAdapter::Impl::makeBlockSizeAcceptable(size_t blockSize) const
|
cannam@64
|
214 {
|
cannam@64
|
215 if (blockSize < 2) {
|
cannam@64
|
216
|
cannam@70
|
217 std::cerr << "WARNING: Vamp::HostExt::PluginInputDomainAdapter::Impl::initialise: blocksize < 2 not" << std::endl
|
cannam@64
|
218 << "supported, increasing from " << blockSize << " to 2" << std::endl;
|
cannam@64
|
219 blockSize = 2;
|
cannam@64
|
220
|
cannam@64
|
221 } else if (blockSize & (blockSize-1)) {
|
cannam@64
|
222
|
cannam@64
|
223 // not a power of two, can't handle that with our current fft
|
cannam@64
|
224 // implementation
|
cannam@64
|
225
|
cannam@64
|
226 size_t nearest = blockSize;
|
cannam@64
|
227 size_t power = 0;
|
cannam@64
|
228 while (nearest > 1) {
|
cannam@64
|
229 nearest >>= 1;
|
cannam@64
|
230 ++power;
|
cannam@64
|
231 }
|
cannam@64
|
232 nearest = 1;
|
cannam@64
|
233 while (power) {
|
cannam@64
|
234 nearest <<= 1;
|
cannam@64
|
235 --power;
|
cannam@64
|
236 }
|
cannam@64
|
237
|
cannam@64
|
238 if (blockSize - nearest > (nearest*2) - blockSize) {
|
cannam@64
|
239 nearest = nearest*2;
|
cannam@64
|
240 }
|
cannam@64
|
241
|
cannam@70
|
242 std::cerr << "WARNING: Vamp::HostExt::PluginInputDomainAdapter::Impl::initialise: non-power-of-two\nblocksize " << blockSize << " not supported, using blocksize " << nearest << " instead" << std::endl;
|
cannam@64
|
243 blockSize = nearest;
|
cannam@64
|
244 }
|
cannam@64
|
245
|
cannam@64
|
246 return blockSize;
|
cannam@64
|
247 }
|
cannam@64
|
248
|
cannam@68
|
249 // for some visual studii apparently
|
cannam@68
|
250 #ifndef M_PI
|
cannam@68
|
251 #define M_PI 3.14159265358979232846
|
cannam@68
|
252 #endif
|
cannam@68
|
253
|
cannam@64
|
254 Plugin::FeatureSet
|
cannam@70
|
255 PluginInputDomainAdapter::Impl::process(const float *const *inputBuffers,
|
cannam@70
|
256 RealTime timestamp)
|
cannam@64
|
257 {
|
cannam@64
|
258 if (m_plugin->getInputDomain() == TimeDomain) {
|
cannam@64
|
259 return m_plugin->process(inputBuffers, timestamp);
|
cannam@64
|
260 }
|
cannam@64
|
261
|
cannam@64
|
262 // The timestamp supplied should be (according to the Vamp::Plugin
|
cannam@64
|
263 // spec) the time of the start of the time-domain input block.
|
cannam@64
|
264 // However, we want to pass to the plugin an FFT output calculated
|
cannam@64
|
265 // from the block of samples _centred_ on that timestamp.
|
cannam@64
|
266 //
|
cannam@64
|
267 // We have two options:
|
cannam@64
|
268 //
|
cannam@64
|
269 // 1. Buffer the input, calculating the fft of the values at the
|
cannam@64
|
270 // passed-in block minus blockSize/2 rather than starting at the
|
cannam@64
|
271 // passed-in block. So each time we call process on the plugin,
|
cannam@64
|
272 // we are passing in the same timestamp as was passed to our own
|
cannam@64
|
273 // process plugin, but not (the frequency domain representation
|
cannam@64
|
274 // of) the same set of samples. Advantages: avoids confusion in
|
cannam@64
|
275 // the host by ensuring the returned values have timestamps
|
cannam@64
|
276 // comparable with that passed in to this function (in fact this
|
cannam@64
|
277 // is pretty much essential for one-value-per-block outputs);
|
cannam@64
|
278 // consistent with hosts such as SV that deal with the
|
cannam@64
|
279 // frequency-domain transform themselves. Disadvantages: means
|
cannam@64
|
280 // making the not necessarily correct assumption that the samples
|
cannam@64
|
281 // preceding the first official block are all zero (or some other
|
cannam@64
|
282 // known value).
|
cannam@64
|
283 //
|
cannam@64
|
284 // 2. Increase the passed-in timestamps by half the blocksize. So
|
cannam@64
|
285 // when we call process, we are passing in the frequency domain
|
cannam@64
|
286 // representation of the same set of samples as passed to us, but
|
cannam@64
|
287 // with a different timestamp. Advantages: simplicity; avoids
|
cannam@64
|
288 // iffy assumption mentioned above. Disadvantages: inconsistency
|
cannam@64
|
289 // with SV in cases where stepSize != blockSize/2; potential
|
cannam@64
|
290 // confusion arising from returned timestamps being calculated
|
cannam@64
|
291 // from the adjusted input timestamps rather than the original
|
cannam@64
|
292 // ones (and inaccuracy where the returned timestamp is implied,
|
cannam@64
|
293 // as in one-value-per-block).
|
cannam@64
|
294 //
|
cannam@64
|
295 // Neither way is ideal, but I don't think either is strictly
|
cannam@64
|
296 // incorrect either. I think this is just a case where the same
|
cannam@64
|
297 // plugin can legitimately produce differing results from the same
|
cannam@64
|
298 // input data, depending on how that data is packaged.
|
cannam@64
|
299 //
|
cannam@64
|
300 // We'll go for option 2, adjusting the timestamps. Note in
|
cannam@64
|
301 // particular that this means some results can differ from those
|
cannam@64
|
302 // produced by SV.
|
cannam@64
|
303
|
cannam@65
|
304 // std::cerr << "PluginInputDomainAdapter: sampleRate " << m_inputSampleRate << ", blocksize " << m_blockSize << ", adjusting time from " << timestamp;
|
cannam@64
|
305
|
cannam@68
|
306 timestamp = timestamp + RealTime::frame2RealTime
|
cannam@68
|
307 (m_blockSize/2, int(m_inputSampleRate + 0.5));
|
cannam@64
|
308
|
cannam@65
|
309 // std::cerr << " to " << timestamp << std::endl;
|
cannam@64
|
310
|
cannam@64
|
311 for (size_t c = 0; c < m_channels; ++c) {
|
cannam@64
|
312
|
cannam@64
|
313 for (size_t i = 0; i < m_blockSize; ++i) {
|
cannam@64
|
314 // Hanning window
|
cannam@64
|
315 m_ri[i] = double(inputBuffers[c][i])
|
cannam@64
|
316 * (0.50 - 0.50 * cos((2 * M_PI * i)
|
cannam@64
|
317 / m_blockSize));
|
cannam@64
|
318 }
|
cannam@64
|
319
|
cannam@64
|
320 for (size_t i = 0; i < m_blockSize/2; ++i) {
|
cannam@64
|
321 // FFT shift
|
cannam@64
|
322 double value = m_ri[i];
|
cannam@64
|
323 m_ri[i] = m_ri[i + m_blockSize/2];
|
cannam@64
|
324 m_ri[i + m_blockSize/2] = value;
|
cannam@64
|
325 }
|
cannam@64
|
326
|
cannam@64
|
327 fft(m_blockSize, false, m_ri, 0, m_ro, m_io);
|
cannam@64
|
328
|
cannam@64
|
329 for (size_t i = 0; i <= m_blockSize/2; ++i) {
|
cannam@64
|
330 m_freqbuf[c][i * 2] = m_ro[i];
|
cannam@64
|
331 m_freqbuf[c][i * 2 + 1] = m_io[i];
|
cannam@64
|
332 }
|
cannam@64
|
333 }
|
cannam@64
|
334
|
cannam@64
|
335 return m_plugin->process(m_freqbuf, timestamp);
|
cannam@64
|
336 }
|
cannam@64
|
337
|
cannam@64
|
338 void
|
cannam@70
|
339 PluginInputDomainAdapter::Impl::fft(unsigned int n, bool inverse,
|
cannam@70
|
340 double *ri, double *ii, double *ro, double *io)
|
cannam@64
|
341 {
|
cannam@64
|
342 if (!ri || !ro || !io) return;
|
cannam@64
|
343
|
cannam@64
|
344 unsigned int bits;
|
cannam@64
|
345 unsigned int i, j, k, m;
|
cannam@64
|
346 unsigned int blockSize, blockEnd;
|
cannam@64
|
347
|
cannam@64
|
348 double tr, ti;
|
cannam@64
|
349
|
cannam@64
|
350 if (n < 2) return;
|
cannam@64
|
351 if (n & (n-1)) return;
|
cannam@64
|
352
|
cannam@64
|
353 double angle = 2.0 * M_PI;
|
cannam@64
|
354 if (inverse) angle = -angle;
|
cannam@64
|
355
|
cannam@64
|
356 for (i = 0; ; ++i) {
|
cannam@64
|
357 if (n & (1 << i)) {
|
cannam@64
|
358 bits = i;
|
cannam@64
|
359 break;
|
cannam@64
|
360 }
|
cannam@64
|
361 }
|
cannam@64
|
362
|
cannam@64
|
363 static unsigned int tableSize = 0;
|
cannam@64
|
364 static int *table = 0;
|
cannam@64
|
365
|
cannam@64
|
366 if (tableSize != n) {
|
cannam@64
|
367
|
cannam@64
|
368 delete[] table;
|
cannam@64
|
369
|
cannam@64
|
370 table = new int[n];
|
cannam@64
|
371
|
cannam@64
|
372 for (i = 0; i < n; ++i) {
|
cannam@64
|
373
|
cannam@64
|
374 m = i;
|
cannam@64
|
375
|
cannam@64
|
376 for (j = k = 0; j < bits; ++j) {
|
cannam@64
|
377 k = (k << 1) | (m & 1);
|
cannam@64
|
378 m >>= 1;
|
cannam@64
|
379 }
|
cannam@64
|
380
|
cannam@64
|
381 table[i] = k;
|
cannam@64
|
382 }
|
cannam@64
|
383
|
cannam@64
|
384 tableSize = n;
|
cannam@64
|
385 }
|
cannam@64
|
386
|
cannam@64
|
387 if (ii) {
|
cannam@64
|
388 for (i = 0; i < n; ++i) {
|
cannam@64
|
389 ro[table[i]] = ri[i];
|
cannam@64
|
390 io[table[i]] = ii[i];
|
cannam@64
|
391 }
|
cannam@64
|
392 } else {
|
cannam@64
|
393 for (i = 0; i < n; ++i) {
|
cannam@64
|
394 ro[table[i]] = ri[i];
|
cannam@64
|
395 io[table[i]] = 0.0;
|
cannam@64
|
396 }
|
cannam@64
|
397 }
|
cannam@64
|
398
|
cannam@64
|
399 blockEnd = 1;
|
cannam@64
|
400
|
cannam@64
|
401 for (blockSize = 2; blockSize <= n; blockSize <<= 1) {
|
cannam@64
|
402
|
cannam@64
|
403 double delta = angle / (double)blockSize;
|
cannam@64
|
404 double sm2 = -sin(-2 * delta);
|
cannam@64
|
405 double sm1 = -sin(-delta);
|
cannam@64
|
406 double cm2 = cos(-2 * delta);
|
cannam@64
|
407 double cm1 = cos(-delta);
|
cannam@64
|
408 double w = 2 * cm1;
|
cannam@64
|
409 double ar[3], ai[3];
|
cannam@64
|
410
|
cannam@64
|
411 for (i = 0; i < n; i += blockSize) {
|
cannam@64
|
412
|
cannam@64
|
413 ar[2] = cm2;
|
cannam@64
|
414 ar[1] = cm1;
|
cannam@64
|
415
|
cannam@64
|
416 ai[2] = sm2;
|
cannam@64
|
417 ai[1] = sm1;
|
cannam@64
|
418
|
cannam@64
|
419 for (j = i, m = 0; m < blockEnd; j++, m++) {
|
cannam@64
|
420
|
cannam@64
|
421 ar[0] = w * ar[1] - ar[2];
|
cannam@64
|
422 ar[2] = ar[1];
|
cannam@64
|
423 ar[1] = ar[0];
|
cannam@64
|
424
|
cannam@64
|
425 ai[0] = w * ai[1] - ai[2];
|
cannam@64
|
426 ai[2] = ai[1];
|
cannam@64
|
427 ai[1] = ai[0];
|
cannam@64
|
428
|
cannam@64
|
429 k = j + blockEnd;
|
cannam@64
|
430 tr = ar[0] * ro[k] - ai[0] * io[k];
|
cannam@64
|
431 ti = ar[0] * io[k] + ai[0] * ro[k];
|
cannam@64
|
432
|
cannam@64
|
433 ro[k] = ro[j] - tr;
|
cannam@64
|
434 io[k] = io[j] - ti;
|
cannam@64
|
435
|
cannam@64
|
436 ro[j] += tr;
|
cannam@64
|
437 io[j] += ti;
|
cannam@64
|
438 }
|
cannam@64
|
439 }
|
cannam@64
|
440
|
cannam@64
|
441 blockEnd = blockSize;
|
cannam@64
|
442 }
|
cannam@64
|
443
|
cannam@64
|
444 if (inverse) {
|
cannam@64
|
445
|
cannam@64
|
446 double denom = (double)n;
|
cannam@64
|
447
|
cannam@64
|
448 for (i = 0; i < n; i++) {
|
cannam@64
|
449 ro[i] /= denom;
|
cannam@64
|
450 io[i] /= denom;
|
cannam@64
|
451 }
|
cannam@64
|
452 }
|
cannam@64
|
453 }
|
cannam@64
|
454
|
cannam@64
|
455 }
|
cannam@64
|
456
|
cannam@64
|
457 }
|
cannam@64
|
458
|