Mercurial > hg > qm-dsp
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 } |