49 #include <QtConcurrent/QtConcurrent>
74 QList<QList<GaborAtom> > AdaptiveMp::matching_pursuit(MatrixXd signal, qint32 max_iterations, qreal epsilon,
bool fix_phase =
false, qint32 boost = 0,
75 qint32 simplex_it = 1E3, qreal simplex_reflection = 1.0, qreal simplex_expansion = 0.2 ,
76 qreal simplex_contraction = 0.5, qreal simplex_full_contraction = 0.5,
bool trial_separation =
false)
79 std::cout <<
"\nAdaptive Matching Pursuit Algorithm started...\n";
81 max_it = max_iterations;
82 Eigen::FFT<double> fft;
83 MatrixXd residuum = signal;
84 qint32 sample_count = signal.rows();
85 qint32 channel_count = signal.cols();
88 qreal residuum_energy = 0;
89 qreal energy_threshold = 0;
92 for(qint32 channel = 0; channel < channel_count; channel++)
94 for(qint32 sample = 0; sample < sample_count; sample++)
95 signal_energy += (signal(sample, channel) * signal(sample, channel));
97 energy_threshold = 0.01 * epsilon * signal_energy;
98 residuum_energy = signal_energy;
100 std::cout <<
"absolute energy of signal: " << residuum_energy <<
"\n";
102 while(it < max_iterations && (energy_threshold < residuum_energy) && sample_count > 1)
104 channel_count = channel_count * (boost / 100.0);
105 if(boost == 0 || channel_count == 0)
111 VectorXd max_scalar_product = VectorXd::Zero(channel_count);
113 qint32 p = floor(sample_count / 2);
115 gabor_Atom->sample_count = sample_count;
116 gabor_Atom->energy = 0;
119 while(s < sample_count)
122 p = floor(sample_count / 2);
124 VectorXcd fft_envelope = RowVectorXcd::Zero(sample_count);
125 fft.fwd(fft_envelope, envelope);
127 while(k < sample_count/2)
129 p = floor(sample_count/2);
131 VectorXcd modulated_resid = VectorXcd::Zero(sample_count);
132 VectorXcd fft_modulated_resid = VectorXcd::Zero(sample_count);
133 VectorXcd fft_m_e_resid = VectorXcd::Zero(sample_count);
134 VectorXd corr_coeffs = VectorXd::Zero(sample_count);
137 for(qint32 chn = 0; chn < channel_count; chn++)
139 qint32 max_index = 0;
142 p = floor(sample_count/2);
145 for(qint32 l = 0; l< sample_count; l++)
146 modulated_resid[l] = residuum(l, chn) * modulation[l];
148 fft.fwd(fft_modulated_resid, modulated_resid);
150 for( qint32 m = 0; m < sample_count; m++)
151 fft_m_e_resid[m] = fft_modulated_resid[m] * conj(fft_envelope[m]);
153 fft.inv(corr_coeffs, fft_m_e_resid);
154 maximum = corr_coeffs[0];
157 for(qint32 i = 1; i < corr_coeffs.rows(); i++)
158 if(maximum < corr_coeffs[i])
160 maximum = corr_coeffs[i];
165 if(max_index >= p) p = max_index - p + 1;
166 else p = max_index + p;
168 VectorXd atom_parameters =
calculate_atom(sample_count, s, p, k, chn, residuum, RETURNPARAMETERS, fix_phase);
169 qreal temp_scalar_product = 0;
170 if(trial_separation) temp_scalar_product = max_scalar_product[chn];
171 else temp_scalar_product = max_scalar_product[0];
173 if(abs(atom_parameters[4]) >= abs(temp_scalar_product))
176 gabor_Atom->scale = atom_parameters[0];
177 gabor_Atom->translation = atom_parameters[1];
178 gabor_Atom->modulation = atom_parameters[2];
179 gabor_Atom->phase = atom_parameters[3];
180 gabor_Atom->max_scalar_product = atom_parameters[4];
181 gabor_Atom->bm_channel = chn;
185 max_scalar_product[chn] = atom_parameters[4];
187 if(atoms_in_chns.length() < channel_count)
188 atoms_in_chns.append(*gabor_Atom);
190 atoms_in_chns.replace(chn, *gabor_Atom);
193 max_scalar_product[0] = atom_parameters[4];
198 k += pow(2.0,(-j))*sample_count/2;
204 std::cout <<
"\n" <<
"===============" <<
" found parameters " << it + 1 <<
"===============" <<
":\n\n"<<
205 "scale: " << gabor_Atom->scale <<
" trans: " << gabor_Atom->translation <<
206 " modu: " << gabor_Atom->modulation <<
" phase: " << gabor_Atom->phase <<
" sclr_prdct: " << gabor_Atom->max_scalar_product <<
"\n";
211 p = floor(sample_count / 2);
212 j = floor(log2(sample_count));
216 for(qint32 chn = 0; chn < channel_count; chn++)
220 while(k < sample_count / 2)
222 VectorXd parameters_no_envelope =
calculate_atom(sample_count, s, p, k, chn, residuum, RETURNPARAMETERS, fix_phase);
224 qreal temp_scalar_product = 0;
225 if(trial_separation) temp_scalar_product = max_scalar_product[chn];
226 else temp_scalar_product = max_scalar_product[0];
227 if(abs(parameters_no_envelope[4]) > abs(temp_scalar_product))
231 gabor_Atom->scale = parameters_no_envelope[0];
232 gabor_Atom->translation = parameters_no_envelope[1];
233 gabor_Atom->modulation = parameters_no_envelope[2];
234 gabor_Atom->phase = parameters_no_envelope[3];
235 gabor_Atom->max_scalar_product = parameters_no_envelope[4];
236 gabor_Atom->bm_channel = chn;
240 max_scalar_product[chn] = parameters_no_envelope[4];
241 atoms_in_chns.replace(chn, *gabor_Atom);
244 max_scalar_product[0] = parameters_no_envelope[4];
247 k += pow(2.0,(-j))*sample_count/2;
252 std::cout <<
" after comparison to NoEnvelope " <<
":\n"<<
"scale: " << gabor_Atom->scale <<
" trans: " << gabor_Atom->translation <<
253 " modu: " << gabor_Atom->modulation <<
" phase: " << gabor_Atom->phase <<
" sclr_prdct: " << gabor_Atom->max_scalar_product <<
"\n\n";
257 if(trial_separation && simplex_it != 0)
259 for(qint32 chn = 0; chn < atoms_in_chns.length(); chn++)
261 *gabor_Atom = atoms_in_chns.at(chn);
262 simplex_maximisation(simplex_it, simplex_reflection, simplex_expansion, simplex_contraction, simplex_full_contraction,
263 gabor_Atom, max_scalar_product, sample_count, fix_phase, residuum, trial_separation, chn);
267 else if(simplex_it != 0)
268 simplex_maximisation(simplex_it, simplex_reflection, simplex_expansion, simplex_contraction, simplex_full_contraction,
269 gabor_Atom, max_scalar_product, sample_count, fix_phase, residuum, trial_separation, gabor_Atom->bm_channel);
272 channel_count = signal.cols();
273 VectorXd best_match = VectorXd::Zero(sample_count);
276 for(qint32 chn = 0; chn < channel_count; chn++)
281 *gabor_Atom = atoms_in_chns.at(chn);
282 gabor_Atom->energy = 0;
283 best_match = gabor_Atom->
create_real(gabor_Atom->sample_count, gabor_Atom->scale, gabor_Atom->translation,
284 gabor_Atom->modulation, gabor_Atom->phase);
288 VectorXd channel_params =
calculate_atom(sample_count, gabor_Atom->scale, gabor_Atom->translation,
289 gabor_Atom->modulation, chn, residuum, RETURNPARAMETERS, fix_phase);
290 gabor_Atom->phase_list.append(channel_params[3]);
292 gabor_Atom->max_scalar_list.append(channel_params[4]);
294 best_match = gabor_Atom->
create_real(gabor_Atom->sample_count, gabor_Atom->scale, gabor_Atom->translation,
295 gabor_Atom->modulation, gabor_Atom->phase_list.at(chn));
299 for(qint32 j = 0; j < gabor_Atom->sample_count; j++)
301 if(!trial_separation)
303 residuum(j,chn) -= gabor_Atom->max_scalar_list.at(chn) * best_match[j];
304 gabor_Atom->energy += pow(gabor_Atom->max_scalar_list.at(chn) * best_match[j], 2);
308 residuum(j,chn) -= gabor_Atom->max_scalar_product * best_match[j];
309 gabor_Atom->energy += pow(gabor_Atom->max_scalar_product * best_match[j], 2);
314 atoms_in_chns.replace(chn, *gabor_Atom);
315 residuum_energy -= atoms_in_chns.at(chn).energy;
316 current_energy += atoms_in_chns.at(chn).energy;
320 if(!trial_separation)
322 residuum_energy -= gabor_Atom->energy;
323 current_energy += gabor_Atom->energy;
332 std::cout <<
"absolute energy of residue: " << residuum_energy <<
"\n";
334 if(!trial_separation)
335 atoms_in_chns.append(*gabor_Atom);
337 atom_list.append(atoms_in_chns);
339 atoms_in_chns.clear();
343 emit current_result(it, max_it, current_energy, signal_energy, residuum, atom_list, fix_dict_list);
345 if( QThread::currentThread()->isInterruptionRequested())
352 std::cout <<
"\nAdaptive Matching Pursuit Algorithm finished.\n";
353 emit finished_calc();
361 VectorXcd modulation = VectorXcd::Zero(N);
363 for(qint32 n = 0; n < N; n++)
365 modulation[n] = std::polar(1 / sqrt(qreal(N)), 2 * PI * k / qreal(N) * qreal(n));
372 VectorXd
AdaptiveMp::calculate_atom(qint32 sample_count, qreal scale, qint32 translation, qreal modulation, qint32 channel, MatrixXd residuum, ReturnValue return_value = RETURNATOM,
bool fix_phase =
false)
377 VectorXcd complex_gabor_atom = gabor_Atom->
create_complex(sample_count, scale, translation, modulation);
380 std::complex<double> inner_product(0, 0);
382 if(fix_phase ==
false)
384 for(qint32 i = 0; i < sample_count; i++)
385 inner_product += residuum(i, channel) * conj(complex_gabor_atom[i]);
389 for(qint32 chn = 0; chn < residuum.cols(); chn++)
391 for(qint32 i = 0; i < sample_count; i++)
392 inner_product += residuum(i, chn) * conj(complex_gabor_atom[i]);
394 if(residuum.cols() != 0)
395 inner_product /= residuum.cols();
399 phase = std::arg(inner_product);
400 if (phase < 0) phase = 2 * PI - phase;
401 VectorXd real_gabor_atom = gabor_Atom->
create_real(sample_count, scale, translation, modulation, phase);
407 case RETURNPARAMETERS:
410 qreal scalar_product = 0;
412 for(qint32 i = 0; i < sample_count; i++)
413 scalar_product += real_gabor_atom[i] * residuum(i, channel);
415 VectorXd atom_parameters = VectorXd::Zero(5);
417 atom_parameters[0] = scale;
418 atom_parameters[1] = translation;
419 atom_parameters[2] = modulation;
420 atom_parameters[3] = phase;
421 atom_parameters[4] = scalar_product;
424 return atom_parameters;
427 case RETURNATOM: {
return real_gabor_atom;}
429 default:
return real_gabor_atom;
436 GaborAtom *gabor_Atom, VectorXd max_scalar_product, qint32 sample_count,
bool fix_phase, MatrixXd residuum,
bool trial_separation, qint32 chn)
440 std::vector<double> init;
442 init.push_back(gabor_Atom->scale);
443 init.push_back(gabor_Atom->translation);
444 init.push_back(gabor_Atom->modulation);
446 double tol = 1E8 * std::numeric_limits<double>::epsilon();
447 std::vector<std::vector<double> > x = std::vector<std::vector<double> >();
448 qint32 iterations = simplex_it;
449 qint32 N = init.size();
451 VectorXd atom_fxc_params = VectorXd::Zero(5);
453 qreal a = simplex_reflection, b = simplex_expansion, g = simplex_contraction, h = simplex_full_contraction;
459 std::vector<double> xcentroid_old(N,0);
460 std::vector<double> xcentroid_new(N,0);
461 std::vector<double> vf(N+1,0);
462 qint32 x1 = 0, xn = 0, xnp1 = 0;
471 std::vector<double> del( init );
472 std::transform(del.begin(), del.end(), del.begin(),
473 std::bind2nd( std::divides<double>() , 20) );
476 for(qint32 i = 0; i < N; ++i)
478 std::vector<double> tmp( init );
486 std::transform(init.begin(), init.end(), xcentroid_old.begin(), std::bind2nd(std::multiplies<double>(), N+1) );
490 for(cnt=0; cnt<iterations; ++cnt)
492 for(qint32 i=0; i < N+1; ++i)
494 VectorXd atom_fx = VectorXd::Zero(sample_count);
496 if(gabor_Atom->scale == sample_count && gabor_Atom->translation == floor(sample_count / 2))
497 atom_fx =
calculate_atom(sample_count, sample_count, floor(sample_count / 2), x[i][2], chn, residuum, RETURNATOM, fix_phase);
500 atom_fx =
calculate_atom(sample_count, x[i][0], x[i][1], x[i][2], chn, residuum, RETURNATOM, fix_phase);
504 for(qint32 k = 0; k < atom_fx.rows(); k++)
506 target -=atom_fx[k]*residuum(k, chn);
514 for(quint32 i=0; i < vf.size(); ++i)
516 if(vf[i]<vf[x1]) x1 = i;
517 if(vf[i]>vf[xnp1]) xnp1 = i;
522 for(quint32 i=0; i<vf.size();++i)
if(vf[i]<vf[xnp1] && vf[i]>vf[xn]) xn=i;
526 std::vector<double> xg(N, 0);
528 for(quint32 i=0; i<x.size(); ++i)
if(i!=xnp1) std::transform(xg.begin(), xg.end(), x[i].begin(), xg.begin(), std::plus<double>() );
530 std::transform(xg.begin(), xg.end(), x[xnp1].begin(), xcentroid_new.begin(), std::plus<double>());
531 std::transform(xg.begin(), xg.end(), xg.begin(), std::bind2nd(std::divides<double>(), N) );
538 for(qint32 i=0; i<N; ++i) diff += fabs(xcentroid_old[i]-xcentroid_new[i]);
540 if (diff/N < tol)
break;
541 else xcentroid_old.swap(xcentroid_new);
544 std::vector<double> xr(N,0);
546 for( qint32 i=0; i<N; ++i) xr[i]=xg[i]+a*(xg[i]-x[xnp1][i]);
549 VectorXd atom_fxr = VectorXd::Zero(sample_count);
551 if(gabor_Atom->scale == sample_count && gabor_Atom->translation == floor(sample_count / 2))
552 atom_fxr =
calculate_atom(sample_count, sample_count, floor(sample_count / 2), xr[2], chn, residuum, RETURNATOM, fix_phase);
555 atom_fxr =
calculate_atom(sample_count, xr[0], xr[1], xr[2], chn, residuum, RETURNATOM, fix_phase);
559 for(qint32 k = 0; k < atom_fxr.rows(); k++) fxr -=atom_fxr[k]*residuum(k,chn);
563 if(vf[x1]<=fxr && fxr<=vf[xn]) std::copy(xr.begin(), xr.end(), x[xnp1].begin());
568 std::vector<double> xe(N,0);
570 for( qint32 i=0; i<N; ++i) xe[i]=xr[i]+b*(xr[i]-xg[i]);
572 VectorXd atom_fxe = VectorXd::Zero(sample_count);
574 if(gabor_Atom->scale == sample_count && gabor_Atom->translation == floor(sample_count / 2))
575 atom_fxe =
calculate_atom(sample_count, sample_count, floor(sample_count / 2), xe[2], chn, residuum, RETURNATOM, fix_phase);
578 atom_fxe =
calculate_atom(sample_count, xe[0], xe[1], xe[2], chn, residuum, RETURNATOM, fix_phase);
582 for(qint32 k = 0; k < atom_fxe.rows(); k++) fxe -=atom_fxe[k]*residuum(k,chn);
584 if( fxe < fxr ) std::copy(xe.begin(), xe.end(), x[xnp1].begin() );
585 else std::copy(xr.begin(), xr.end(), x[xnp1].begin() );
589 else if( fxr > vf[xn] )
591 std::vector<double> xc(N,0);
593 for( qint32 i=0; i<N; ++i)
594 xc[i]=xg[i]+g*(x[xnp1][i]-xg[i]);
596 if(gabor_Atom->scale == sample_count && gabor_Atom->translation == floor(sample_count / 2))
597 atom_fxc_params =
AdaptiveMp::calculate_atom(sample_count, sample_count, floor(sample_count / 2), xc[2], chn, residuum, RETURNPARAMETERS, fix_phase);
602 VectorXd atom_fxc = gabor_Atom->
create_real(gabor_Atom->sample_count, atom_fxc_params[0], atom_fxc_params[1], atom_fxc_params[2], atom_fxc_params[3]);
604 atom_fxc_params[4] = 0;
606 for(qint32 i = 0; i < sample_count; i++)
607 atom_fxc_params[4] += atom_fxc[i] * residuum(i, chn);
612 for(qint32 k = 0; k < atom_fxc.rows(); k++)
613 fxc -=atom_fxc[k]*residuum(k,chn);
616 std::copy(xc.begin(), xc.end(), x[xnp1].begin() );
619 for( quint32 i=0; i<x.size(); ++i )
621 for(qint32 j=0; j<N; ++j)
622 x[i][j] = x[x1][j] + h * ( x[i][j]-x[x1][j] );
630 if(gabor_Atom->scale == sample_count && gabor_Atom->translation == floor(sample_count / 2))
631 atom_fxc_params =
AdaptiveMp::calculate_atom(sample_count, sample_count, floor(sample_count / 2), x[x1][2], chn, residuum, RETURNPARAMETERS, fix_phase);
634 atom_fxc_params =
AdaptiveMp::calculate_atom(sample_count, x[x1][0], x[x1][1], x[x1][2], chn, residuum, RETURNPARAMETERS, fix_phase);
636 qreal temp_scalar_product = 0;
637 if(trial_separation) temp_scalar_product = max_scalar_product[chn];
638 else temp_scalar_product = max_scalar_product[0];
640 if(abs(atom_fxc_params[4]) > abs(temp_scalar_product) && atom_fxc_params[1] < sample_count && atom_fxc_params[1] > 0)
644 gabor_Atom->scale = atom_fxc_params[0];
645 gabor_Atom->translation = atom_fxc_params[1];
646 gabor_Atom->modulation = atom_fxc_params[2];
647 gabor_Atom->phase = atom_fxc_params[3];
648 gabor_Atom->max_scalar_product = atom_fxc_params[4];
651 atoms_in_chns.replace(chn, *gabor_Atom);
657 std::cout<<
"Simplex Iteration limit of "<<iterations<<
" achieved, result may not be optimal\n";
660 std::cout <<
" after simplex optimization " <<
":\n"<<
"scale: " << gabor_Atom->scale <<
" trans: " << gabor_Atom->translation <<
661 " modu: " << gabor_Atom->modulation <<
" phase: " << gabor_Atom->phase <<
" sclr_prdct: " << gabor_Atom->max_scalar_product <<
"\n\n";
667 void AdaptiveMp::recieve_input(Eigen::MatrixXd signal, qint32 max_iterations, qreal epsilon,
bool fix_phase =
false, qint32 boost = 0, qint32 simplex_it = 1E3,
668 qreal simplex_reflection = 1.0, qreal simplex_expansion = 0.2, qreal simplex_contraction = 0.5, qreal simplex_full_contraction = 0.5,
bool trial_separation =
false)
670 matching_pursuit(signal, max_iterations, epsilon, fix_phase, boost, simplex_it, simplex_reflection, simplex_expansion, simplex_contraction, simplex_full_contraction, trial_separation);
GaborAtom used in adaptive MP Algorithm.
static VectorXd gauss_function(qint32 sample_count, qreal scale, quint32 translation)
ADAPIVEMP class declaration, providing the implemetation of the Matching Pursuit Algorithm introduced...
VectorXcd create_complex(qint32 sample_count, qreal scale, quint32 translation, qreal modulation)
VectorXd create_real(qint32 sample_count, qreal scale, quint32 translation, qreal modulation, qreal phase)
static VectorXd calculate_atom(qint32 sample_count, qreal scale, qint32 translation, qreal modulation, qint32 channel, MatrixXd residuum, ReturnValue return_value, bool fix_phase)
void simplex_maximisation(qint32 simplex_it, qreal simplex_reflection, qreal simplex_expansion, qreal simplex_contraction, qreal simplex_full_contraction, GaborAtom *gabor_Atom, VectorXd max_scalar_product, qint32 sample_count, bool fix_phase, MatrixXd residuum, bool trial_separation, qint32 chn)
VectorXcd modulation_function(qint32 N, qreal k)