MNE-CPP  beta 1.0
adaptivemp.cpp
Go to the documentation of this file.
1 //=============================================================================================================
41 //*************************************************************************************************************
42 //=============================================================================================================
43 // INCLUDES
44 //=============================================================================================================
45 
46 #include "adaptivemp.h"
47 #include <vector>
48 #include <QFuture>
49 #include <QtConcurrent/QtConcurrent>
50 
51 //*************************************************************************************************************
52 //=============================================================================================================
53 // USED NAMESPACES
54 //=============================================================================================================
55 
56 using namespace UTILSLIB;
57 
58 //*************************************************************************************************************
59 //=============================================================================================================
60 // DEFINE MEMBER METHODS
61 //=============================================================================================================
62 
64  : it(0)
65  , max_it(0)
66  , signal_energy(0)
67  , current_energy(0)
68 {
69 
70 }
71 
72 //*************************************************************************************************************
73 
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)
77 {
78  //stop_running = false;
79  std::cout << "\nAdaptive Matching Pursuit Algorithm started...\n";
80 
81  max_it = max_iterations;
82  Eigen::FFT<double> fft;
83  MatrixXd residuum = signal; //residuum initialised with signal
84  qint32 sample_count = signal.rows();
85  qint32 channel_count = signal.cols();
86 
87  signal_energy = 0;
88  qreal residuum_energy = 0;
89  qreal energy_threshold = 0;
90 
91  //calculate signal_energy
92  for(qint32 channel = 0; channel < channel_count; channel++)
93  {
94  for(qint32 sample = 0; sample < sample_count; sample++)
95  signal_energy += (signal(sample, channel) * signal(sample, channel));
96 
97  energy_threshold = 0.01 * epsilon * signal_energy;
98  residuum_energy = signal_energy;
99  }
100  std::cout << "absolute energy of signal: " << residuum_energy << "\n";
101 
102  while(it < max_iterations && (energy_threshold < residuum_energy) && sample_count > 1)
103  {
104  channel_count = channel_count * (boost / 100.0); //reducing the number of observed channels in the algorithm to increase speed performance
105  if(boost == 0 || channel_count == 0)
106  channel_count = 1;
107 
108  //variables for dyadic sampling
109  qreal s = 1; //scale
110  qint32 j = 1;
111  VectorXd max_scalar_product = VectorXd::Zero(channel_count); //inner product for choosing the best matching atom
112  qreal k = 0; //for modulation 2*pi*k/N
113  qint32 p = floor(sample_count / 2); //translation
114  GaborAtom *gabor_Atom = new GaborAtom();
115  gabor_Atom->sample_count = sample_count;
116  gabor_Atom->energy = 0;
117  qreal phase = 0;
118 
119  while(s < sample_count)
120  {
121  k = 0; //for modulation 2*pi*k/N
122  p = floor(sample_count / 2); //translation
123  VectorXd envelope = GaborAtom::gauss_function(sample_count, s, p);
124  VectorXcd fft_envelope = RowVectorXcd::Zero(sample_count);
125  fft.fwd(fft_envelope, envelope);
126 
127  while(k < sample_count/2)
128  {
129  p = floor(sample_count/2);
130  VectorXcd modulation = modulation_function(sample_count, k);
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);
135 
136  //iteration for multichannel, depending on boost setting
137  for(qint32 chn = 0; chn < channel_count; chn++)
138  {
139  qint32 max_index = 0;
140  qreal maximum = 0;
141  phase = 0;
142  p = floor(sample_count/2);//here is difference to dr. gratkowski´s code (he didn´t reset parameter p)
143 
144  //complex correlation of signal and sinus-modulated gaussfunction
145  for(qint32 l = 0; l< sample_count; l++)
146  modulated_resid[l] = residuum(l, chn) * modulation[l];
147 
148  fft.fwd(fft_modulated_resid, modulated_resid);
149 
150  for( qint32 m = 0; m < sample_count; m++)
151  fft_m_e_resid[m] = fft_modulated_resid[m] * conj(fft_envelope[m]);
152 
153  fft.inv(corr_coeffs, fft_m_e_resid);
154  maximum = corr_coeffs[0];
155 
156  //find index of maximum correlation-coefficient to use in translation
157  for(qint32 i = 1; i < corr_coeffs.rows(); i++)
158  if(maximum < corr_coeffs[i])
159  {
160  maximum = corr_coeffs[i];
161  max_index = i;
162  }
163 
164  //adapting translation p to create atomtranslation correctly
165  if(max_index >= p) p = max_index - p + 1;
166  else p = max_index + p;
167 
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];
172 
173  if(abs(atom_parameters[4]) >= abs(temp_scalar_product))
174  {
175  //set highest scalarproduct, in comparison to best matching atom
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;
182 
183  if(trial_separation)
184  {
185  max_scalar_product[chn] = atom_parameters[4];
186 
187  if(atoms_in_chns.length() < channel_count)
188  atoms_in_chns.append(*gabor_Atom);
189  else
190  atoms_in_chns.replace(chn, *gabor_Atom);
191  }
192  else
193  max_scalar_product[0] = atom_parameters[4];
194 
195  }
196 
197  }
198  k += pow(2.0,(-j))*sample_count/2;
199 
200  }
201  j++;
202  s = pow(2.0,j);
203  }
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";
207 
208  //replace atoms with s==N and p = floor(N/2) by such atoms that do not have an envelope
209  k = 0;
210  s = sample_count;
211  p = floor(sample_count / 2);
212  j = floor(log2(sample_count));//log(sample_count) / log(2));
213  phase = 0;
214 
215  //iteration for multichannel, depending on boost setting
216  for(qint32 chn = 0; chn < channel_count; chn++)
217  {
218  k = 0;
219 
220  while(k < sample_count / 2)
221  {
222  VectorXd parameters_no_envelope = calculate_atom(sample_count, s, p, k, chn, residuum, RETURNPARAMETERS, fix_phase);
223 
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))
228  {
229  //set highest scalarproduct, in comparison to best matching atom
230 
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;
237 
238  if(trial_separation)
239  {
240  max_scalar_product[chn] = parameters_no_envelope[4];
241  atoms_in_chns.replace(chn, *gabor_Atom);
242  }
243  else
244  max_scalar_product[0] = parameters_no_envelope[4];
245 
246  }
247  k += pow(2.0,(-j))*sample_count/2;
248 
249  }
250 
251  }
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";
254 
255  //simplexfunction to find minimum of target among parameters s, p, k
256 
257  if(trial_separation && simplex_it != 0)
258  {
259  for(qint32 chn = 0; chn < atoms_in_chns.length(); chn++)
260  {
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);
264 
265  }
266  }
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);
270 
271  //calc multichannel parameters phase and max_scalar_product
272  channel_count = signal.cols();
273  VectorXd best_match = VectorXd::Zero(sample_count);
274 
275 
276  for(qint32 chn = 0; chn < channel_count; chn++)
277  {
278 
279  if(trial_separation)
280  {
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);
285  }
286  else
287  {
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]);
291 
292  gabor_Atom->max_scalar_list.append(channel_params[4]);
293 
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));
296  }
297 
298  //substract best matching Atom from Residuum in each channel
299  for(qint32 j = 0; j < gabor_Atom->sample_count; j++)
300  {
301  if(!trial_separation)
302  {
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);
305  }
306  else
307  {
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);
310  }
311  }
312  if(trial_separation)
313  {
314  atoms_in_chns.replace(chn, *gabor_Atom); // change energy
315  residuum_energy -= atoms_in_chns.at(chn).energy;// / channel_count;;
316  current_energy += atoms_in_chns.at(chn).energy;// / channel_count;;
317  }
318  }
319 
320  if(!trial_separation)
321  {
322  residuum_energy -= gabor_Atom->energy;
323  current_energy += gabor_Atom->energy;
324  }
325  else
326  /*for(qint32 i = 0; i < atoms_in_chns.length(); i++)
327  {
328  residuum_energy -= atoms_in_chns.at(i).energy / channel_count;
329  current_energy += atoms_in_chns.at(i).energy / channel_count;
330  }*/
331 
332  std::cout << "absolute energy of residue: " << residuum_energy << "\n";
333 
334  if(!trial_separation)
335  atoms_in_chns.append(*gabor_Atom);
336 
337  atom_list.append(atoms_in_chns);
338 
339  atoms_in_chns.clear();
340  delete gabor_Atom;
341  it++;
342 
343  emit current_result(it, max_it, current_energy, signal_energy, residuum, atom_list, fix_dict_list);
344 
345  if( QThread::currentThread()->isInterruptionRequested())
346  {
347  send_warning(10);
348  break;
349  }
350 
351  }//end iterations
352  std::cout << "\nAdaptive Matching Pursuit Algorithm finished.\n";
353  emit finished_calc();
354  return atom_list;
355 }
356 
357 //*************************************************************************************************************
358 
359 VectorXcd AdaptiveMp::modulation_function(qint32 N, qreal k)
360 {
361  VectorXcd modulation = VectorXcd::Zero(N);
362 
363  for(qint32 n = 0; n < N; n++)
364  {
365  modulation[n] = std::polar(1 / sqrt(qreal(N)), 2 * PI * k / qreal(N) * qreal(n));
366  }
367  return modulation;
368 }
369 
370 //*************************************************************************************************************
371 
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)
373 {
374  GaborAtom *gabor_Atom = new GaborAtom();
375  qreal phase = 0;
376  //create complex Gaboratom
377  VectorXcd complex_gabor_atom = gabor_Atom->create_complex(sample_count, scale, translation, modulation);
378 
379  //calculate Inner Product: preparation to find the parameter phase
380  std::complex<double> inner_product(0, 0);
381 
382  if(fix_phase == false)
383  {
384  for(qint32 i = 0; i < sample_count; i++)
385  inner_product += residuum(i, channel) * conj(complex_gabor_atom[i]);
386  }
387  else
388  {
389  for(qint32 chn = 0; chn < residuum.cols(); chn++)
390  {
391  for(qint32 i = 0; i < sample_count; i++)
392  inner_product += residuum(i, chn) * conj(complex_gabor_atom[i]);
393  }
394  if(residuum.cols() != 0)
395  inner_product /= residuum.cols();
396  }
397 
398  //calculate phase to create realGaborAtoms
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);
402 
403  delete gabor_Atom;
404 
405  switch(return_value)
406  {
407  case RETURNPARAMETERS:
408  {
409 
410  qreal scalar_product = 0;
411 
412  for(qint32 i = 0; i < sample_count; i++)
413  scalar_product += real_gabor_atom[i] * residuum(i, channel);
414 
415  VectorXd atom_parameters = VectorXd::Zero(5);
416 
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;
422 
423 
424  return atom_parameters;
425  }
426 
427  case RETURNATOM: {return real_gabor_atom;} //returns normalized realGaborAtom
428 
429  default: return real_gabor_atom;// only to avoid compiler warnings
430  }
431 }
432 
433 //*************************************************************************************************************
434 
435 void AdaptiveMp::simplex_maximisation(qint32 simplex_it, qreal simplex_reflection, qreal simplex_expansion, qreal simplex_contraction, qreal simplex_full_contraction,
436  GaborAtom *gabor_Atom, VectorXd max_scalar_product, qint32 sample_count, bool fix_phase, MatrixXd residuum, bool trial_separation, qint32 chn)
437 {
438  //Maximisation Simplex Algorithm implemented by Botao Jia, adapted to the MP Algorithm by Martin Henfling. Copyright (C) 2010 Botao Jia
439  //ToDo: change to clean use of EIGEN, @present its mixed with Namespace std and <vector>
440  std::vector<double> init;
441 
442  init.push_back(gabor_Atom->scale);
443  init.push_back(gabor_Atom->translation);
444  init.push_back(gabor_Atom->modulation);
445 
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(); //space dimension
450 
451  VectorXd atom_fxc_params = VectorXd::Zero(5); //initialisation for contraction coefficients
452 
453  qreal a = simplex_reflection, b = simplex_expansion, g = simplex_contraction, h = simplex_full_contraction;
454  //coefficients default a = 1, b = 0.2, g = 0.5, h = 0.5
455  //a: reflection -> xr step away from worst siplex found
456  //b: expansion -> xe if better with a so go in this direction with b
457  //g: contraction -> xc calc new worst point an bring closer to middle of simplex
458  //h: full contraction to x1
459  std::vector<double> xcentroid_old(N,0); //simplex center * (N+1)
460  std::vector<double> xcentroid_new(N,0); //simplex center * (N+1)
461  std::vector<double> vf(N+1,0); //f evaluated at simplex vertices
462  qint32 x1 = 0, xn = 0, xnp1 = 0; //x1: f(x1) = min { f(x1), f(x2)...f(x_{n+1} }
463  //xnp1: f(xnp1) = max { f(x1), f(x2)...f(x_{n+1} }
464  //xn: f(xn)<f(xnp1) && f(xn)> all other f(x_i)
465  qint32 cnt = 0; //iteration step number
466 
467  if(x.size()== 0) //if no initial simplex is specified
468  {
469  //construct the trial simplex
470  //based upon the initial guess parameters
471  std::vector<double> del( init );
472  std::transform(del.begin(), del.end(), del.begin(),
473  std::bind2nd( std::divides<double>() , 20) );//'20' is picked
474  //assuming initial trail close to true
475 
476  for(qint32 i = 0; i < N; ++i)
477  {
478  std::vector<double> tmp( init );
479  tmp[i] += del[i];
480  x.push_back( tmp );
481  }
482 
483  x.push_back(init);//x.size()=N+1, x[i].size()=N
484 
485  //xcentriod
486  std::transform(init.begin(), init.end(), xcentroid_old.begin(), std::bind2nd(std::multiplies<double>(), N+1) );
487  }//constructing the simplex finished
488 
489  //optimization begins
490  for(cnt=0; cnt<iterations; ++cnt)
491  {
492  for(qint32 i=0; i < N+1; ++i)
493  {
494  VectorXd atom_fx = VectorXd::Zero(sample_count);
495 
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);
498 
499  else
500  atom_fx = calculate_atom(sample_count, x[i][0], x[i][1], x[i][2], chn, residuum, RETURNATOM, fix_phase);
501 
502  //create targetfunction of realGaborAtom and Residuum
503  double target = 0;
504  for(qint32 k = 0; k < atom_fx.rows(); k++)
505  {
506  target -=atom_fx[k]*residuum(k, chn); //ToDo: old residuum(k,0)
507  }
508 
509  vf[i] = target;
510  }
511 
512  x1=0; xn=0; xnp1=0;//find index of max, second max, min of vf.
513 
514  for(quint32 i=0; i < vf.size(); ++i)
515  {
516  if(vf[i]<vf[x1]) x1 = i;
517  if(vf[i]>vf[xnp1]) xnp1 = i;
518  }
519 
520  xn = x1;
521 
522  for(quint32 i=0; i<vf.size();++i) if(vf[i]<vf[xnp1] && vf[i]>vf[xn]) xn=i;
523 
524  //x1, xn, xnp1 are found
525 
526  std::vector<double> xg(N, 0);//xg: centroid of the N best vertexes
527 
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>() );
529 
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) );
532  //xg found, xcentroid_new updated
533 
534  //termination condition
535  double diff=0; //calculate the difference of the simplex centers
536 
537  //see if the difference is less than the termination criteria
538  for(qint32 i=0; i<N; ++i) diff += fabs(xcentroid_old[i]-xcentroid_new[i]);
539 
540  if (diff/N < tol) break; //terminate the optimizer
541  else xcentroid_old.swap(xcentroid_new); //update simplex center
542 
543  //reflection:
544  std::vector<double> xr(N,0);
545 
546  for( qint32 i=0; i<N; ++i) xr[i]=xg[i]+a*(xg[i]-x[xnp1][i]);
547  //reflection, xr found
548 
549  VectorXd atom_fxr = VectorXd::Zero(sample_count);
550 
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);
553 
554  else
555  atom_fxr = calculate_atom(sample_count, xr[0], xr[1], xr[2], chn, residuum, RETURNATOM, fix_phase);
556 
557  //create targetfunction of realGaborAtom and Residuum
558  double fxr = 0;
559  for(qint32 k = 0; k < atom_fxr.rows(); k++) fxr -=atom_fxr[k]*residuum(k,chn);//ToDo: old residuum(k,0)
560 
561  //double fxr = target;//record function at xr
562 
563  if(vf[x1]<=fxr && fxr<=vf[xn]) std::copy(xr.begin(), xr.end(), x[xnp1].begin());
564 
565  //expansion:
566  else if(fxr<vf[x1])
567  {
568  std::vector<double> xe(N,0);
569 
570  for( qint32 i=0; i<N; ++i) xe[i]=xr[i]+b*(xr[i]-xg[i]);
571 
572  VectorXd atom_fxe = VectorXd::Zero(sample_count);
573 
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);
576 
577  else
578  atom_fxe = calculate_atom(sample_count, xe[0], xe[1], xe[2], chn, residuum, RETURNATOM, fix_phase);
579 
580  //create targetfunction of realGaborAtom and Residuum
581  double fxe = 0;
582  for(qint32 k = 0; k < atom_fxe.rows(); k++) fxe -=atom_fxe[k]*residuum(k,chn);//ToDo: old residuum(k,0)
583 
584  if( fxe < fxr ) std::copy(xe.begin(), xe.end(), x[xnp1].begin() );
585  else std::copy(xr.begin(), xr.end(), x[xnp1].begin() );
586  }//expansion finished, xe is not used outside the scope
587 
588  //contraction:
589  else if( fxr > vf[xn] )
590  {
591  std::vector<double> xc(N,0);
592 
593  for( qint32 i=0; i<N; ++i)
594  xc[i]=xg[i]+g*(x[xnp1][i]-xg[i]);
595 
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);
598 
599  else
600  atom_fxc_params = AdaptiveMp::calculate_atom(sample_count, xc[0], xc[1], xc[2], chn, residuum, RETURNPARAMETERS, fix_phase);
601 
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]);
603 
604  atom_fxc_params[4] = 0;
605 
606  for(qint32 i = 0; i < sample_count; i++)
607  atom_fxc_params[4] += atom_fxc[i] * residuum(i, chn);
608 
609  //create targetfunction of realGaborAtom and Residuum
610  double fxc = 0;
611 
612  for(qint32 k = 0; k < atom_fxc.rows(); k++)
613  fxc -=atom_fxc[k]*residuum(k,chn);//ToDo: old residuum(k,0)
614 
615  if( fxc < vf[xnp1] )
616  std::copy(xc.begin(), xc.end(), x[xnp1].begin() );
617 
618  else
619  for( quint32 i=0; i<x.size(); ++i )
620  if( i!=x1 )
621  for(qint32 j=0; j<N; ++j)
622  x[i][j] = x[x1][j] + h * ( x[i][j]-x[x1][j] );
623  }//contraction finished, xc is not used outside the scope
624  }//optimization is finished
625  //end Maximisation Copyright (C) 2010 Botao Jia
626 
627  if(iterations != 0)
628  {
629 
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);
632 
633  else
634  atom_fxc_params = AdaptiveMp::calculate_atom(sample_count, x[x1][0], x[x1][1], x[x1][2], chn, residuum, RETURNPARAMETERS, fix_phase);
635 
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];
639 
640  if(abs(atom_fxc_params[4]) > abs(temp_scalar_product) && atom_fxc_params[1] < sample_count && atom_fxc_params[1] > 0)//ToDo: find a way to make the simplex not running out of bounds
641  {
642  //set highest scalarproduct, in comparison to best matching atom
643 
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];
649 
650  if(trial_separation)
651  atoms_in_chns.replace(chn, *gabor_Atom);
652  }
653 
654  if(cnt==iterations)//max number of iteration achieves before tol is satisfied
655  {
656  send_warning(11);
657  std::cout<<"Simplex Iteration limit of "<<iterations<<" achieved, result may not be optimal\n";
658  }
659 
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";
662  }
663 }
664 
665 //*************************************************************************************************************
666 
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)
669 {
670  matching_pursuit(signal, max_iterations, epsilon, fix_phase, boost, simplex_it, simplex_reflection, simplex_expansion, simplex_contraction, simplex_full_contraction, trial_separation);
671 }
672 
673 //*************************************************************************************************************
674 
676 {
677 }
GaborAtom used in adaptive MP Algorithm.
Definition: atom.h:192
static VectorXd gauss_function(qint32 sample_count, qreal scale, quint32 translation)
Definition: atom.cpp:135
ADAPIVEMP class declaration, providing the implemetation of the Matching Pursuit Algorithm introduced...
VectorXcd create_complex(qint32 sample_count, qreal scale, quint32 translation, qreal modulation)
Definition: atom.cpp:150
VectorXd create_real(qint32 sample_count, qreal scale, quint32 translation, qreal modulation, qreal phase)
Definition: atom.cpp:185
static VectorXd calculate_atom(qint32 sample_count, qreal scale, qint32 translation, qreal modulation, qint32 channel, MatrixXd residuum, ReturnValue return_value, bool fix_phase)
Definition: adaptivemp.cpp:372
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)
Definition: adaptivemp.cpp:435
VectorXcd modulation_function(qint32 N, qreal k)
Definition: adaptivemp.cpp:359