← Back to branch summary

~ubuntu-branches/ubuntu/utopic/ardour3/utopic

« back to all changes in this revision

Viewing changes to libs/qm-dsp/dsp/chromagram/ConstantQ.cpp

  • Committer: Package Import Robot
  • Author(s): Felipe Sateler
  • Date: 2013-09-21 19:05:02 UTC
  • Revision ID: package-import@ubuntu.com-20130921190502-8gsftrku6jnzhd7v
Tags: upstream-3.4~dfsg
Import upstream version 3.4~dfsg

Show diffs side-by-side

added added

removed removed

Lines of Context:
libs/qm-dsp/dsp/chromagram/ConstantQ.cpp
 
1
/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
 
2
/*
 
3
    QM DSP Library
 
4
 
 
5
    Centre for Digital Music, Queen Mary, University of London.
 
6
    This file 2005-2006 Christian Landone.
 
7
 
 
8
    This program is free software; you can redistribute it and/or
 
9
    modify it under the terms of the GNU General Public License as
 
10
    published by the Free Software Foundation; either version 2 of the
 
11
    License, or (at your option) any later version.  See the file
 
12
    COPYING included with this distribution for more information.
 
13
*/
 
14
 
 
15
#include "ConstantQ.h"
 
16
#include "dsp/transforms/FFT.h"
 
17
 
 
18
#include <iostream>
 
19
 
 
20
#ifdef NOT_DEFINED
 
21
// see note in CQprecalc
 
22
 
 
23
#include "CQprecalc.cpp"
 
24
 
 
25
static bool push_precalculated(int uk, int fftlength,
 
26
                               std::vector<unsigned> &is,
 
27
                               std::vector<unsigned> &js,
 
28
                               std::vector<double> &real,
 
29
                               std::vector<double> &imag)
 
30
{
 
31
    if (uk == 76 && fftlength == 16384) {
 
32
        push_76_16384(is, js, real, imag);
 
33
        return true;
 
34
    }
 
35
    if (uk == 144 && fftlength == 4096) {
 
36
        push_144_4096(is, js, real, imag);
 
37
        return true;
 
38
    }
 
39
    if (uk == 65 && fftlength == 2048) {
 
40
        push_65_2048(is, js, real, imag);
 
41
        return true;
 
42
    }
 
43
    if (uk == 84 && fftlength == 65536) {
 
44
        push_84_65536(is, js, real, imag);
 
45
        return true;
 
46
    }
 
47
    return false;
 
48
}
 
49
#endif
 
50
 
 
51
//---------------------------------------------------------------------------
 
52
// nextpow2 returns the smallest integer n such that 2^n >= x.
 
53
static double nextpow2(double x) {
 
54
    double y = ceil(log(x)/log(2.0));
 
55
    return(y);
 
56
}
 
57
 
 
58
static double squaredModule(const double & xx, const double & yy) {
 
59
    return xx*xx + yy*yy;
 
60
}
 
61
 
 
62
//----------------------------------------------------------------------------
 
63
 
 
64
ConstantQ::ConstantQ( CQConfig Config ) :
 
65
    m_sparseKernel(0)
 
66
{
 
67
    initialise( Config );
 
68
}
 
69
 
 
70
ConstantQ::~ConstantQ()
 
71
{
 
72
    deInitialise();
 
73
}
 
74
 
 
75
//----------------------------------------------------------------------------
 
76
void ConstantQ::sparsekernel()
 
77
{
 
78
//    std::cerr << "ConstantQ: initialising sparse kernel, uK = " << m_uK << ", FFTLength = " << m_FFTLength << "...";
 
79
 
 
80
    SparseKernel *sk = new SparseKernel();
 
81
 
 
82
#ifdef NOT_DEFINED
 
83
    if (push_precalculated(m_uK, m_FFTLength,
 
84
                           sk->is, sk->js, sk->real, sk->imag)) {
 
85
//        std::cerr << "using precalculated kernel" << std::endl;
 
86
        m_sparseKernel = sk;
 
87
        return;
 
88
    }
 
89
#endif
 
90
 
 
91
    //generates spectral kernel matrix (upside down?)
 
92
    // initialise temporal kernel with zeros, twice length to deal w. complex numbers
 
93
 
 
94
    double* hammingWindowRe = new double [ m_FFTLength ];
 
95
    double* hammingWindowIm = new double [ m_FFTLength ];
 
96
    double* transfHammingWindowRe = new double [ m_FFTLength ];
 
97
    double* transfHammingWindowIm = new double [ m_FFTLength ];
 
98
 
 
99
    for (unsigned u=0; u < m_FFTLength; u++) 
 
100
    {
 
101
        hammingWindowRe[u] = 0;
 
102
        hammingWindowIm[u] = 0;
 
103
    }
 
104
 
 
105
    // Here, fftleng*2 is a guess of the number of sparse cells in the matrix
 
106
    // The matrix K x fftlength but the non-zero cells are an antialiased
 
107
    // square root function. So mostly is a line, with some grey point.
 
108
    sk->is.reserve( m_FFTLength*2 );
 
109
    sk->js.reserve( m_FFTLength*2 );
 
110
    sk->real.reserve( m_FFTLength*2 );
 
111
    sk->imag.reserve( m_FFTLength*2 );
 
112
        
 
113
    // for each bin value K, calculate temporal kernel, take its fft to
 
114
    //calculate the spectral kernel then threshold it to make it sparse and 
 
115
    //add it to the sparse kernels matrix
 
116
    double squareThreshold = m_CQThresh * m_CQThresh;
 
117
 
 
118
    FFT m_FFT(m_FFTLength);
 
119
        
 
120
    for (unsigned k = m_uK; k--; ) 
 
121
    {
 
122
        for (unsigned u=0; u < m_FFTLength; u++) 
 
123
        {
 
124
            hammingWindowRe[u] = 0;
 
125
            hammingWindowIm[u] = 0;
 
126
        }
 
127
        
 
128
        // Computing a hamming window
 
129
        const unsigned hammingLength = (int) ceil( m_dQ * m_FS / ( m_FMin * pow(2,((double)(k))/(double)m_BPO)));
 
130
 
 
131
        unsigned origin = m_FFTLength/2 - hammingLength/2;
 
132
 
 
133
        for (unsigned i=0; i<hammingLength; i++) 
 
134
        {
 
135
            const double angle = 2*PI*m_dQ*i/hammingLength;
 
136
            const double real = cos(angle);
 
137
            const double imag = sin(angle);
 
138
            const double absol = hamming(hammingLength, i)/hammingLength;
 
139
            hammingWindowRe[ origin + i ] = absol*real;
 
140
            hammingWindowIm[ origin + i ] = absol*imag;
 
141
        }
 
142
 
 
143
        for (unsigned i = 0; i < m_FFTLength/2; ++i) {
 
144
            double temp = hammingWindowRe[i];
 
145
            hammingWindowRe[i] = hammingWindowRe[i + m_FFTLength/2];
 
146
            hammingWindowRe[i + m_FFTLength/2] = temp;
 
147
            temp = hammingWindowIm[i];
 
148
            hammingWindowIm[i] = hammingWindowIm[i + m_FFTLength/2];
 
149
            hammingWindowIm[i + m_FFTLength/2] = temp;
 
150
        }
 
151
    
 
152
        //do fft of hammingWindow
 
153
        m_FFT.process( 0, hammingWindowRe, hammingWindowIm, transfHammingWindowRe, transfHammingWindowIm );
 
154
 
 
155
                
 
156
        for (unsigned j=0; j<( m_FFTLength ); j++) 
 
157
        {
 
158
            // perform thresholding
 
159
            const double squaredBin = squaredModule( transfHammingWindowRe[ j ], transfHammingWindowIm[ j ]);
 
160
            if (squaredBin <= squareThreshold) continue;
 
161
                
 
162
            // Insert non-zero position indexes, doubled because they are floats
 
163
            sk->is.push_back(j);
 
164
            sk->js.push_back(k);
 
165
 
 
166
            // take conjugate, normalise and add to array sparkernel
 
167
            sk->real.push_back( transfHammingWindowRe[ j ]/m_FFTLength);
 
168
            sk->imag.push_back(-transfHammingWindowIm[ j ]/m_FFTLength);
 
169
        }
 
170
 
 
171
    }
 
172
 
 
173
    delete [] hammingWindowRe;
 
174
    delete [] hammingWindowIm;
 
175
    delete [] transfHammingWindowRe;
 
176
    delete [] transfHammingWindowIm;
 
177
 
 
178
/*
 
179
    using std::cout;
 
180
    using std::endl;
 
181
 
 
182
    cout.precision(28);
 
183
 
 
184
    int n = sk->is.size();
 
185
    int w = 8;
 
186
    cout << "static unsigned int sk_i_" << m_uK << "_" << m_FFTLength << "[" << n << "] = {" << endl;
 
187
    for (int i = 0; i < n; ++i) {
 
188
        if (i % w == 0) cout << "    ";
 
189
        cout << sk->is[i];
 
190
        if (i + 1 < n) cout << ", ";
 
191
        if (i % w == w-1) cout << endl;
 
192
    };
 
193
    if (n % w != 0) cout << endl;
 
194
    cout << "};" << endl;
 
195
 
 
196
    n = sk->js.size();
 
197
    cout << "static unsigned int sk_j_" << m_uK << "_" << m_FFTLength << "[" << n << "] = {" << endl;
 
198
    for (int i = 0; i < n; ++i) {
 
199
        if (i % w == 0) cout << "    ";
 
200
        cout << sk->js[i];
 
201
        if (i + 1 < n) cout << ", ";
 
202
        if (i % w == w-1) cout << endl;
 
203
    };
 
204
    if (n % w != 0) cout << endl;
 
205
    cout << "};" << endl;
 
206
 
 
207
    w = 2;
 
208
    n = sk->real.size();
 
209
    cout << "static double sk_real_" << m_uK << "_" << m_FFTLength << "[" << n << "] = {" << endl;
 
210
    for (int i = 0; i < n; ++i) {
 
211
        if (i % w == 0) cout << "    ";
 
212
        cout << sk->real[i];
 
213
        if (i + 1 < n) cout << ", ";
 
214
        if (i % w == w-1) cout << endl;
 
215
    };
 
216
    if (n % w != 0) cout << endl;
 
217
    cout << "};" << endl;
 
218
 
 
219
    n = sk->imag.size();
 
220
    cout << "static double sk_imag_" << m_uK << "_" << m_FFTLength << "[" << n << "] = {" << endl;
 
221
    for (int i = 0; i < n; ++i) {
 
222
        if (i % w == 0) cout << "    ";
 
223
        cout << sk->imag[i];
 
224
        if (i + 1 < n) cout << ", ";
 
225
        if (i % w == w-1) cout << endl;
 
226
    };
 
227
    if (n % w != 0) cout << endl;
 
228
    cout << "};" << endl;
 
229
 
 
230
    cout << "static void push_" << m_uK << "_" << m_FFTLength << "(vector<unsigned int> &is, vector<unsigned int> &js, vector<double> &real, vector<double> &imag)" << endl;
 
231
    cout << "{\n    is.reserve(" << n << ");\n";
 
232
    cout << "    js.reserve(" << n << ");\n";
 
233
    cout << "    real.reserve(" << n << ");\n";
 
234
    cout << "    imag.reserve(" << n << ");\n";
 
235
    cout << "    for (int i = 0; i < " << n << "; ++i) {" << endl;
 
236
    cout << "        is.push_back(sk_i_" << m_uK << "_" << m_FFTLength << "[i]);" << endl;
 
237
    cout << "        js.push_back(sk_j_" << m_uK << "_" << m_FFTLength << "[i]);" << endl;
 
238
    cout << "        real.push_back(sk_real_" << m_uK << "_" << m_FFTLength << "[i]);" << endl;
 
239
    cout << "        imag.push_back(sk_imag_" << m_uK << "_" << m_FFTLength << "[i]);" << endl;
 
240
    cout << "    }" << endl;
 
241
    cout << "}" << endl;
 
242
*/
 
243
//    std::cerr << "done\n -> is: " << sk->is.size() << ", js: " << sk->js.size() << ", reals: " << sk->real.size() << ", imags: " << sk->imag.size() << std::endl;
 
244
    
 
245
    m_sparseKernel = sk;
 
246
    return;
 
247
}
 
248
 
 
249
//-----------------------------------------------------------------------------
 
250
double* ConstantQ::process( const double* fftdata )
 
251
{
 
252
    if (!m_sparseKernel) {
 
253
        std::cerr << "ERROR: ConstantQ::process: Sparse kernel has not been initialised" << std::endl;
 
254
        return m_CQdata;
 
255
    }
 
256
 
 
257
    SparseKernel *sk = m_sparseKernel;
 
258
 
 
259
    for (unsigned row=0; row<2*m_uK; row++) 
 
260
    {
 
261
        m_CQdata[ row ] = 0;
 
262
        m_CQdata[ row+1 ] = 0;
 
263
    }
 
264
    const unsigned *fftbin = &(sk->is[0]);
 
265
    const unsigned *cqbin  = &(sk->js[0]);
 
266
    const double   *real   = &(sk->real[0]);
 
267
    const double   *imag   = &(sk->imag[0]);
 
268
    const unsigned int sparseCells = sk->real.size();
 
269
        
 
270
    for (unsigned i = 0; i<sparseCells; i++)
 
271
    {
 
272
        const unsigned row = cqbin[i];
 
273
        const unsigned col = fftbin[i];
 
274
        const double & r1  = real[i];
 
275
        const double & i1  = imag[i];
 
276
        const double & r2  = fftdata[ (2*m_FFTLength) - 2*col - 2 ];
 
277
        const double & i2  = fftdata[ (2*m_FFTLength) - 2*col - 2 + 1 ];
 
278
        // add the multiplication
 
279
        m_CQdata[ 2*row  ] += (r1*r2 - i1*i2);
 
280
        m_CQdata[ 2*row+1] += (r1*i2 + i1*r2);
 
281
    }
 
282
 
 
283
    return m_CQdata;
 
284
}
 
285
 
 
286
 
 
287
void ConstantQ::initialise( CQConfig Config )
 
288
{
 
289
    m_FS = Config.FS;
 
290
    m_FMin = Config.min;                // min freq
 
291
    m_FMax = Config.max;                // max freq
 
292
    m_BPO = Config.BPO;         // bins per octave
 
293
    m_CQThresh = Config.CQThresh;// ConstantQ threshold for kernel generation
 
294
 
 
295
    m_dQ = 1/(pow(2,(1/(double)m_BPO))-1);      // Work out Q value for Filter bank
 
296
    m_uK = (unsigned int) ceil(m_BPO * log(m_FMax/m_FMin)/log(2.0));    // No. of constant Q bins
 
297
 
 
298
//    std::cerr << "ConstantQ::initialise: rate = " << m_FS << ", fmin = " << m_FMin << ", fmax = " << m_FMax << ", bpo = " << m_BPO << ", K = " << m_uK << ", Q = " << m_dQ << std::endl;
 
299
 
 
300
    // work out length of fft required for this constant Q Filter bank
 
301
    m_FFTLength = (int) pow(2, nextpow2(ceil( m_dQ*m_FS/m_FMin )));
 
302
 
 
303
    m_hop = m_FFTLength/8; // <------ hop size is window length divided by 32
 
304
 
 
305
//    std::cerr << "ConstantQ::initialise: -> fft length = " << m_FFTLength << ", hop = " << m_hop << std::endl;
 
306
 
 
307
    // allocate memory for cqdata
 
308
    m_CQdata = new double [2*m_uK];
 
309
}
 
310
 
 
311
void ConstantQ::deInitialise()
 
312
{
 
313
    delete [] m_CQdata;
 
314
    delete m_sparseKernel;
 
315
}
 
316
 
 
317
void ConstantQ::process(const double *FFTRe, const double* FFTIm,
 
318
                        double *CQRe, double *CQIm)
 
319
{
 
320
    if (!m_sparseKernel) {
 
321
        std::cerr << "ERROR: ConstantQ::process: Sparse kernel has not been initialised" << std::endl;
 
322
        return;
 
323
    }
 
324
 
 
325
    SparseKernel *sk = m_sparseKernel;
 
326
 
 
327
    for (unsigned row=0; row<m_uK; row++) 
 
328
    {
 
329
        CQRe[ row ] = 0;
 
330
        CQIm[ row ] = 0;
 
331
    }
 
332
 
 
333
    const unsigned *fftbin = &(sk->is[0]);
 
334
    const unsigned *cqbin  = &(sk->js[0]);
 
335
    const double   *real   = &(sk->real[0]);
 
336
    const double   *imag   = &(sk->imag[0]);
 
337
    const unsigned int sparseCells = sk->real.size();
 
338
        
 
339
    for (unsigned i = 0; i<sparseCells; i++)
 
340
    {
 
341
        const unsigned row = cqbin[i];
 
342
        const unsigned col = fftbin[i];
 
343
        const double & r1  = real[i];
 
344
        const double & i1  = imag[i];
 
345
        const double & r2  = FFTRe[ m_FFTLength - col - 1 ];
 
346
        const double & i2  = FFTIm[ m_FFTLength - col - 1 ];
 
347
        // add the multiplication
 
348
        CQRe[ row ] += (r1*r2 - i1*i2);
 
349
        CQIm[ row ] += (r1*i2 + i1*r2);
 
350
    }
 
351
}