comparison dsp/chromagram/ConstantQ.cpp @ 483:fdaa63607c15

Untabify, indent, tidy
author Chris Cannam <cannam@all-day-breakfast.com>
date Fri, 31 May 2019 11:54:32 +0100
parents 73fc1de3254a
children 5998ee1042d3
comparison
equal deleted inserted replaced
482:cbe668c7d724 483:fdaa63607c15
54 double* hammingWindowRe = new double [ m_FFTLength ]; 54 double* hammingWindowRe = new double [ m_FFTLength ];
55 double* hammingWindowIm = new double [ m_FFTLength ]; 55 double* hammingWindowIm = new double [ m_FFTLength ];
56 double* transfHammingWindowRe = new double [ m_FFTLength ]; 56 double* transfHammingWindowRe = new double [ m_FFTLength ];
57 double* transfHammingWindowIm = new double [ m_FFTLength ]; 57 double* transfHammingWindowIm = new double [ m_FFTLength ];
58 58
59 for (unsigned u=0; u < m_FFTLength; u++) 59 for (unsigned u=0; u < m_FFTLength; u++) {
60 { 60 hammingWindowRe[u] = 0;
61 hammingWindowRe[u] = 0; 61 hammingWindowIm[u] = 0;
62 hammingWindowIm[u] = 0;
63 } 62 }
64 63
65 // Here, fftleng*2 is a guess of the number of sparse cells in the matrix 64 // Here, fftleng*2 is a guess of the number of sparse cells in the matrix
66 // The matrix K x fftlength but the non-zero cells are an antialiased 65 // The matrix K x fftlength but the non-zero cells are an antialiased
67 // square root function. So mostly is a line, with some grey point. 66 // square root function. So mostly is a line, with some grey point.
68 sk->is.reserve( m_FFTLength*2 ); 67 sk->is.reserve( m_FFTLength*2 );
69 sk->js.reserve( m_FFTLength*2 ); 68 sk->js.reserve( m_FFTLength*2 );
70 sk->real.reserve( m_FFTLength*2 ); 69 sk->real.reserve( m_FFTLength*2 );
71 sk->imag.reserve( m_FFTLength*2 ); 70 sk->imag.reserve( m_FFTLength*2 );
72 71
73 // for each bin value K, calculate temporal kernel, take its fft to 72 // for each bin value K, calculate temporal kernel, take its fft to
74 //calculate the spectral kernel then threshold it to make it sparse and 73 //calculate the spectral kernel then threshold it to make it sparse and
75 //add it to the sparse kernels matrix 74 //add it to the sparse kernels matrix
76 double squareThreshold = m_CQThresh * m_CQThresh; 75 double squareThreshold = m_CQThresh * m_CQThresh;
77 76
78 FFT m_FFT(m_FFTLength); 77 FFT m_FFT(m_FFTLength);
79 78
80 for (unsigned k = m_uK; k--; ) 79 for (unsigned k = m_uK; k--; ) {
81 { 80 for (unsigned u=0; u < m_FFTLength; u++) {
82 for (unsigned u=0; u < m_FFTLength; u++)
83 {
84 hammingWindowRe[u] = 0; 81 hammingWindowRe[u] = 0;
85 hammingWindowIm[u] = 0; 82 hammingWindowIm[u] = 0;
86 } 83 }
87 84
88 // Computing a hamming window 85 // Computing a hamming window
89 const unsigned hammingLength = (int) ceil( m_dQ * m_FS / ( m_FMin * pow(2,((double)(k))/(double)m_BPO))); 86 const unsigned hammingLength = (int) ceil( m_dQ * m_FS / ( m_FMin * pow(2,((double)(k))/(double)m_BPO)));
90 87
91 // cerr << "k = " << k << ", q = " << m_dQ << ", m_FMin = " << m_FMin << ", hammingLength = " << hammingLength << " (rounded up from " << (m_dQ * m_FS / ( m_FMin * pow(2,((double)(k))/(double)m_BPO))) << ")" << endl; 88 // cerr << "k = " << k << ", q = " << m_dQ << ", m_FMin = " << m_FMin << ", hammingLength = " << hammingLength << " (rounded up from " << (m_dQ * m_FS / ( m_FMin * pow(2,((double)(k))/(double)m_BPO))) << ")" << endl;
92 89
93 90
94 unsigned origin = m_FFTLength/2 - hammingLength/2; 91 unsigned origin = m_FFTLength/2 - hammingLength/2;
95 92
96 for (unsigned i=0; i<hammingLength; i++) 93 for (unsigned i=0; i<hammingLength; i++) {
97 { 94 const double angle = 2*PI*m_dQ*i/hammingLength;
98 const double angle = 2*PI*m_dQ*i/hammingLength; 95 const double real = cos(angle);
99 const double real = cos(angle); 96 const double imag = sin(angle);
100 const double imag = sin(angle); 97 const double absol = hamming(hammingLength, i)/hammingLength;
101 const double absol = hamming(hammingLength, i)/hammingLength; 98 hammingWindowRe[ origin + i ] = absol*real;
102 hammingWindowRe[ origin + i ] = absol*real; 99 hammingWindowIm[ origin + i ] = absol*imag;
103 hammingWindowIm[ origin + i ] = absol*imag; 100 }
104 }
105 101
106 for (unsigned i = 0; i < m_FFTLength/2; ++i) { 102 for (unsigned i = 0; i < m_FFTLength/2; ++i) {
107 double temp = hammingWindowRe[i]; 103 double temp = hammingWindowRe[i];
108 hammingWindowRe[i] = hammingWindowRe[i + m_FFTLength/2]; 104 hammingWindowRe[i] = hammingWindowRe[i + m_FFTLength/2];
109 hammingWindowRe[i + m_FFTLength/2] = temp; 105 hammingWindowRe[i + m_FFTLength/2] = temp;
110 temp = hammingWindowIm[i]; 106 temp = hammingWindowIm[i];
111 hammingWindowIm[i] = hammingWindowIm[i + m_FFTLength/2]; 107 hammingWindowIm[i] = hammingWindowIm[i + m_FFTLength/2];
112 hammingWindowIm[i + m_FFTLength/2] = temp; 108 hammingWindowIm[i + m_FFTLength/2] = temp;
113 } 109 }
114 110
115 //do fft of hammingWindow 111 //do fft of hammingWindow
116 m_FFT.process( 0, hammingWindowRe, hammingWindowIm, transfHammingWindowRe, transfHammingWindowIm ); 112 m_FFT.process( 0, hammingWindowRe, hammingWindowIm, transfHammingWindowRe, transfHammingWindowIm );
117 113
118 114 for (unsigned j=0; j<( m_FFTLength ); j++) {
119 for (unsigned j=0; j<( m_FFTLength ); j++) 115 // perform thresholding
120 { 116 const double squaredBin = squaredModule( transfHammingWindowRe[ j ], transfHammingWindowIm[ j ]);
121 // perform thresholding 117 if (squaredBin <= squareThreshold) continue;
122 const double squaredBin = squaredModule( transfHammingWindowRe[ j ], transfHammingWindowIm[ j ]); 118
123 if (squaredBin <= squareThreshold) continue; 119 // Insert non-zero position indexes
124 120 sk->is.push_back(j);
125 // Insert non-zero position indexes 121 sk->js.push_back(k);
126 sk->is.push_back(j); 122
127 sk->js.push_back(k); 123 // take conjugate, normalise and add to array sparkernel
128 124 sk->real.push_back( transfHammingWindowRe[ j ]/m_FFTLength);
129 // take conjugate, normalise and add to array sparkernel 125 sk->imag.push_back(-transfHammingWindowIm[ j ]/m_FFTLength);
130 sk->real.push_back( transfHammingWindowRe[ j ]/m_FFTLength); 126 }
131 sk->imag.push_back(-transfHammingWindowIm[ j ]/m_FFTLength);
132 }
133
134 } 127 }
135 128
136 delete [] hammingWindowRe; 129 delete [] hammingWindowRe;
137 delete [] hammingWindowIm; 130 delete [] hammingWindowIm;
138 delete [] transfHammingWindowRe; 131 delete [] transfHammingWindowRe;
152 return m_CQdata; 145 return m_CQdata;
153 } 146 }
154 147
155 SparseKernel *sk = m_sparseKernel; 148 SparseKernel *sk = m_sparseKernel;
156 149
157 for (unsigned row=0; row<2*m_uK; row++) 150 for (unsigned row=0; row<2*m_uK; row++) {
158 { 151 m_CQdata[ row ] = 0;
159 m_CQdata[ row ] = 0; 152 m_CQdata[ row+1 ] = 0;
160 m_CQdata[ row+1 ] = 0;
161 } 153 }
162 const unsigned *fftbin = &(sk->is[0]); 154 const unsigned *fftbin = &(sk->is[0]);
163 const unsigned *cqbin = &(sk->js[0]); 155 const unsigned *cqbin = &(sk->js[0]);
164 const double *real = &(sk->real[0]); 156 const double *real = &(sk->real[0]);
165 const double *imag = &(sk->imag[0]); 157 const double *imag = &(sk->imag[0]);
166 const unsigned int sparseCells = sk->real.size(); 158 const unsigned int sparseCells = sk->real.size();
167 159
168 for (unsigned i = 0; i<sparseCells; i++) 160 for (unsigned i = 0; i<sparseCells; i++) {
169 { 161 const unsigned row = cqbin[i];
170 const unsigned row = cqbin[i]; 162 const unsigned col = fftbin[i];
171 const unsigned col = fftbin[i];
172 if (col == 0) continue; 163 if (col == 0) continue;
173 const double & r1 = real[i]; 164 const double & r1 = real[i];
174 const double & i1 = imag[i]; 165 const double & i1 = imag[i];
175 const double & r2 = fftdata[ (2*m_FFTLength) - 2*col - 2 ]; 166 const double & r2 = fftdata[ (2*m_FFTLength) - 2*col - 2 ];
176 const double & i2 = fftdata[ (2*m_FFTLength) - 2*col - 2 + 1 ]; 167 const double & i2 = fftdata[ (2*m_FFTLength) - 2*col - 2 + 1 ];
177 // add the multiplication 168 // add the multiplication
178 m_CQdata[ 2*row ] += (r1*r2 - i1*i2); 169 m_CQdata[ 2*row ] += (r1*r2 - i1*i2);
179 m_CQdata[ 2*row+1] += (r1*i2 + i1*r2); 170 m_CQdata[ 2*row+1] += (r1*i2 + i1*r2);
180 } 171 }
181 172
182 return m_CQdata; 173 return m_CQdata;
183 } 174 }
184 175
185 176
186 void ConstantQ::initialise( CQConfig Config ) 177 void ConstantQ::initialise( CQConfig Config )
187 { 178 {
188 m_FS = Config.FS; 179 m_FS = Config.FS;
189 m_FMin = Config.min; // min freq 180 m_FMin = Config.min; // min freq
190 m_FMax = Config.max; // max freq 181 m_FMax = Config.max; // max freq
191 m_BPO = Config.BPO; // bins per octave 182 m_BPO = Config.BPO; // bins per octave
192 m_CQThresh = Config.CQThresh;// ConstantQ threshold for kernel generation 183 m_CQThresh = Config.CQThresh;// ConstantQ threshold for kernel generation
193 184
194 m_dQ = 1/(pow(2,(1/(double)m_BPO))-1); // Work out Q value for Filter bank 185 m_dQ = 1/(pow(2,(1/(double)m_BPO))-1); // Work out Q value for Filter bank
195 m_uK = (unsigned int) ceil(m_BPO * log(m_FMax/m_FMin)/log(2.0)); // No. of constant Q bins 186 m_uK = (unsigned int) ceil(m_BPO * log(m_FMax/m_FMin)/log(2.0)); // No. of constant Q bins
196 187
197 // std::cerr << "ConstantQ::initialise: rate = " << m_FS << ", fmin = " << m_FMin << ", fmax = " << m_FMax << ", bpo = " << m_BPO << ", K = " << m_uK << ", Q = " << m_dQ << std::endl; 188 // std::cerr << "ConstantQ::initialise: rate = " << m_FS << ", fmin = " << m_FMin << ", fmax = " << m_FMax << ", bpo = " << m_BPO << ", K = " << m_uK << ", Q = " << m_dQ << std::endl;
198 189
199 // work out length of fft required for this constant Q Filter bank 190 // work out length of fft required for this constant Q Filter bank
200 m_FFTLength = (int) pow(2, nextpow2(ceil( m_dQ*m_FS/m_FMin ))); 191 m_FFTLength = (int) pow(2, nextpow2(ceil( m_dQ*m_FS/m_FMin )));
221 return; 212 return;
222 } 213 }
223 214
224 SparseKernel *sk = m_sparseKernel; 215 SparseKernel *sk = m_sparseKernel;
225 216
226 for (unsigned row=0; row<m_uK; row++) 217 for (unsigned row=0; row<m_uK; row++) {
227 { 218 CQRe[ row ] = 0;
228 CQRe[ row ] = 0; 219 CQIm[ row ] = 0;
229 CQIm[ row ] = 0;
230 } 220 }
231 221
232 const unsigned *fftbin = &(sk->is[0]); 222 const unsigned *fftbin = &(sk->is[0]);
233 const unsigned *cqbin = &(sk->js[0]); 223 const unsigned *cqbin = &(sk->js[0]);
234 const double *real = &(sk->real[0]); 224 const double *real = &(sk->real[0]);
235 const double *imag = &(sk->imag[0]); 225 const double *imag = &(sk->imag[0]);
236 const unsigned int sparseCells = sk->real.size(); 226 const unsigned int sparseCells = sk->real.size();
237 227
238 for (unsigned i = 0; i<sparseCells; i++) 228 for (unsigned i = 0; i<sparseCells; i++) {
239 { 229 const unsigned row = cqbin[i];
240 const unsigned row = cqbin[i]; 230 const unsigned col = fftbin[i];
241 const unsigned col = fftbin[i];
242 if (col == 0) continue; 231 if (col == 0) continue;
243 const double & r1 = real[i]; 232 const double & r1 = real[i];
244 const double & i1 = imag[i]; 233 const double & i1 = imag[i];
245 const double & r2 = FFTRe[ m_FFTLength - col ]; 234 const double & r2 = FFTRe[ m_FFTLength - col ];
246 const double & i2 = FFTIm[ m_FFTLength - col ]; 235 const double & i2 = FFTIm[ m_FFTLength - col ];
247 // add the multiplication 236 // add the multiplication
248 CQRe[ row ] += (r1*r2 - i1*i2); 237 CQRe[ row ] += (r1*r2 - i1*i2);
249 CQIm[ row ] += (r1*i2 + i1*r2); 238 CQIm[ row ] += (r1*i2 + i1*r2);
250 } 239 }
251 } 240 }