46 #include "fixdictmp.h"
73 void FixDictMp::matching_pursuit(MatrixXd signal, qint32 max_iterations, qreal epsilon, qint32 boost, QString path, qreal delta)
75 Eigen::initParallel();
77 std::cout <<
"\nFixDict Matching Pursuit Algorithm started...\n";
79 QList<Dictionary> parsed_dicts;
81 qint32 sample_count = signal.rows();
82 qint32 channel_count = signal.cols();
85 qreal residuum_energy = 0;
86 qreal energy_threshold = 0;
87 qreal last_energy = 0;
88 bool sample_count_mismatch =
false;
90 this->residuum = signal;
91 parsed_dicts = parse_xml_dict(path);
94 for(qint32 channel = 0; channel < channel_count; channel++)
96 for(qint32 sample = 0; sample < sample_count; sample++)
97 signal_energy += (signal(sample, channel) * signal(sample, channel));
99 energy_threshold = 0.01 * epsilon * signal_energy;
100 residuum_energy = signal_energy;
103 std::cout <<
"absolute energy of signal: " << residuum_energy <<
"\n";
105 while(it < max_iterations && energy_threshold < residuum_energy)
109 QList<find_best_matching> list_of_best;
110 for(qint32 i = 0; i < parsed_dicts.length(); i++)
112 find_best_matching current_best_matching;
113 current_best_matching.pdict = parsed_dicts.at(i);
114 current_best_matching.current_resid = this->residuum;
115 current_best_matching.boost = boost;
116 list_of_best.append(current_best_matching);
119 QFuture<FixDictAtom> mapped_best_matchings = QtConcurrent::mapped(list_of_best, &find_best_matching::parallel_correlation);
120 mapped_best_matchings.waitForFinished();
122 QFuture<FixDictAtom>::const_iterator i;
123 for (i = mapped_best_matchings.constBegin(); i != mapped_best_matchings.constEnd(); i++)
126 if(i == mapped_best_matchings.constBegin())
127 global_best_matching = *i;
129 else if(abs(i->max_scalar_product) > abs(global_best_matching.max_scalar_product))
130 global_best_matching = *i;
133 global_best_matching.display_text = create_display_text(global_best_matching);
135 VectorXd fitted_atom = VectorXd::Zero(this->residuum.rows());
136 VectorXd resized_atom = VectorXd::Zero(this->residuum.rows());
138 if(global_best_matching.atom_samples.rows() > this->residuum.rows())
139 for(qint32 k = 0; k < this->residuum.rows(); k++)
140 resized_atom[k] = global_best_matching.atom_samples[k + floor(global_best_matching.atom_samples.rows() / 2) - floor(this->residuum.rows() / 2)];
141 else resized_atom = global_best_matching.atom_samples;
144 for(qint32 k = 0; k < resized_atom.rows(); k++)
145 if((k + global_best_matching.translation - floor(resized_atom.rows() / 2)) >= 0
146 && (k + global_best_matching.translation - floor(resized_atom.rows() / 2)) < fitted_atom.rows())
147 fitted_atom[(k + global_best_matching.translation - floor(resized_atom.rows() / 2))] += resized_atom[k];
151 norm = fitted_atom.norm();
152 if(norm != 0) fitted_atom /= norm;
154 for(qint32 chn = 0; chn < this->residuum.cols(); chn++)
156 qreal scalarproduct = 0;
157 for(qint32 sample = 0; sample < this->residuum.rows(); sample++)
158 scalarproduct += (this->residuum(sample, chn) * fitted_atom[sample]);
160 global_best_matching.max_scalar_list.append(scalarproduct);
162 for(qint32 k = 0; k < fitted_atom.rows(); k++)
164 this->residuum(k,chn) -= global_best_matching.max_scalar_list.at(chn) * fitted_atom[k];
165 global_best_matching.energy += pow(global_best_matching.max_scalar_list.at(chn) * fitted_atom[k], 2);
169 global_best_matching.atom_samples = fitted_atom;
172 last_energy = residuum_energy;
180 residuum_energy -= global_best_matching.energy;
181 current_energy += global_best_matching.energy;
183 fix_dict_list.append(global_best_matching);
185 if(global_best_matching.sample_count != signal.rows() && !sample_count_mismatch)
187 std::cout <<
"\n=============================================\n"
188 <<
" INFO\n\n(part-)dictionary does not fit the signal!\n\nResults might be better by creating dictionaries\ncontaining atoms of the same length as the signal\n"
189 <<
"\nsignal samples: " << signal.rows() <<
"\n"
190 <<
"atom samples: " << global_best_matching.sample_count <<
"\n"
191 <<
"=============================================\n";
192 sample_count_mismatch =
true;
195 std::cout <<
"\n" <<
"===============" << it + 1 <<
"th atom found" <<
"===============" <<
":\n\n"<<
196 qPrintable(global_best_matching.display_text) <<
"\n\n" <<
"sample_count: " << global_best_matching.sample_count <<
197 " source_dict: " << qPrintable(global_best_matching.dict_source) <<
" Atom ID: " << global_best_matching.id <<
198 "\nsclr_prdct: " << global_best_matching.max_scalar_product <<
"\n" <<
199 "\nabsolute energy of residue: " << residuum_energy <<
"\n";
203 emit current_result(it, max_iterations, current_energy, signal_energy, residuum, adaptive_list, fix_dict_list);
206 if(((last_energy * 100 / signal_energy) - (residuum_energy * 100 / signal_energy)) < delta)
208 std::cout <<
"\n=============================================\n"
209 <<
" ALGORITHM ABORTED\n\ndictionary excludes atoms to reduce further residual energy\n"
210 <<
"=============================================\n";
211 emit send_warning(1);
214 if( QThread::currentThread()->isInterruptionRequested())
222 std::cout <<
"\nFixDict Matching Pursuit Algorithm finished.\n";
223 emit finished_calc();
232 qint32 channel_count = current_resid.cols() * (boost / 100.0);
233 if(boost == 0 || channel_count == 0)
236 Eigen::FFT<double> fft;
237 std::ptrdiff_t max_index;
238 VectorXcd fft_atom = VectorXcd::Zero(current_resid.rows());
241 qreal max_scalar_product = 0;
243 for(qint32 i = 0; i < current_pdict.atoms.length(); i++)
245 VectorXd fitted_atom = VectorXd::Zero(current_resid.rows());
246 qint32 p = floor(current_resid.rows() / 2);
248 VectorXd resized_atom = VectorXd::Zero(current_resid.rows());
250 if(current_pdict.atoms.at(i).atom_samples.rows() > current_resid.rows())
251 for(qint32 k = 0; k < current_resid.rows(); k++)
252 resized_atom[k] = current_pdict.atoms.at(i).atom_samples[k + floor(current_pdict.atoms.at(i).atom_samples.rows() / 2) - floor(current_resid.rows() / 2)];
253 else resized_atom = current_pdict.atoms.at(i).atom_samples;
255 if(resized_atom.rows() < current_resid.rows())
256 for(qint32 k = 0; k < resized_atom.rows(); k++)
257 fitted_atom[(k + p - floor(resized_atom.rows() / 2))] = resized_atom[k];
258 else fitted_atom = resized_atom;
262 norm = fitted_atom.norm();
263 if(norm != 0) fitted_atom /= norm;
265 fft.fwd(fft_atom, fitted_atom);
267 for(qint32 chn = 0; chn < channel_count; chn++)
269 p = floor(current_resid.rows() / 2);
270 VectorXd corr_coeffs = VectorXd::Zero(current_resid.rows());
271 VectorXcd fft_signal = VectorXcd::Zero(current_resid.rows());
272 VectorXcd fft_sig_atom = VectorXcd::Zero(current_resid.rows());
274 fft.fwd(fft_signal, current_resid.col(chn));
276 for( qint32 m = 0; m < current_resid.rows(); m++)
277 fft_sig_atom[m] = fft_signal[m] * conj(fft_atom[m]);
279 fft.inv(corr_coeffs, fft_sig_atom);
282 max_scalar_product = corr_coeffs.maxCoeff(&max_index);
286 best_matching = current_pdict.atoms.at(i);
287 best_matching.max_scalar_product = max_scalar_product;
290 if(max_index >= p && current_resid.rows() % (2) == 0) p = max_index - p;
291 else if(max_index >= p && current_resid.rows() % (2) != 0) p = max_index - p - 1;
292 else p = max_index + p;
294 best_matching.translation = p;
297 else if(abs(max_scalar_product) > abs(best_matching.max_scalar_product))
299 best_matching = current_pdict.atoms.at(i);
300 best_matching.max_scalar_product = max_scalar_product;
303 if(max_index >= p && current_resid.rows() % (2) == 0) p = max_index - p;
304 else if(max_index >= p && current_resid.rows() % (2) != 0) p = max_index - p - 1;
305 else p = max_index + p;
307 best_matching.translation = p;
311 best_matching.atom_formula = current_pdict.atom_formula;
312 best_matching.dict_source = current_pdict.source;
313 best_matching.type = current_pdict.type;
314 best_matching.sample_count = current_pdict.sample_count;
316 return best_matching;
320 QList<Dictionary> FixDictMp::parse_xml_dict(QString path)
322 QFile current_dict(path);
323 std::cout <<
"\nparsing dictionary, please be patient...";
325 QDomDocument dictionary;
326 dictionary.setContent(¤t_dict);
327 QList<Dictionary> parsed_dict;
328 QDomElement current_element = dictionary.documentElement();
329 QDomNodeList node_list = current_element.childNodes();
330 QList<parse_node> pdict_nodes;
331 QList<QDomNode> nodes_listed;
332 bool is_emitted =
false;
334 for(qint32 i = 0; i < node_list.length(); i++)
337 temp_parse.node = node_list.at(i);
338 pdict_nodes.append(temp_parse);
339 nodes_listed.append(node_list.at(i));
343 QFuture<Dictionary> mapped_dicts = QtConcurrent::mapped(pdict_nodes, &parse_node::fill_dict_in_map);
344 mapped_dicts.waitForFinished();
346 QFuture<Dictionary>::const_iterator i;
348 for (i = mapped_dicts.constBegin(); i != mapped_dicts.constEnd(); i++)
353 if(i->sample_count != this->residuum.rows() && !is_emitted)
356 emit send_warning(2);
359 parsed_dict.append(*i);
362 std::cout <<
" done.\n\n";
368 Dictionary FixDictMp::fill_dict(
const QDomNode &pdict)
371 QDomElement atom = pdict.firstChildElement(
"ATOM");
373 current_dict.source = pdict.toElement().attribute(
"source_dict");
374 current_dict.atom_count();
375 current_dict.atom_formula = pdict.toElement().attribute(
"formula");
376 current_dict.sample_count = pdict.toElement().attribute(
"sample_count").toInt();
378 if(current_dict.atom_formula == QString(
"Gaboratom"))
380 current_dict.type = GABORATOM;
382 while(!atom.isNull())
386 current_atom.id = atom.attribute(
"ID").toInt();
387 current_atom.gabor_atom.scale = atom.attribute(
"scale").toDouble();
388 current_atom.gabor_atom.modulation = atom.attribute(
"modu").toDouble();
389 current_atom.gabor_atom.phase = atom.attribute(
"phase").toDouble();
391 if(atom.hasChildNodes())
393 QString sample_string = atom.firstChild().toElement().attribute(
"samples");
394 QStringList sample_list = sample_string.split(
":", QString::SkipEmptyParts);
395 current_atom.atom_samples = VectorXd::Zero(sample_list.length());
396 for(qint32 i = 0; i < sample_list.length(); i++)
397 current_atom.atom_samples[i] = sample_list.at(i).toDouble();
399 current_dict.atoms.append(current_atom);
400 atom = atom.nextSiblingElement(
"ATOM");
403 else if(current_dict.atom_formula == QString(
"Chirpatom"))
405 current_dict.type = CHIRPATOM;
407 while(!atom.isNull())
411 current_atom.id = atom.attribute(
"id").toInt();
412 current_atom.chirp_atom.scale = atom.attribute(
"scale").toDouble();
413 current_atom.chirp_atom.modulation = atom.attribute(
"modu").toDouble();
414 current_atom.chirp_atom.phase = atom.attribute(
"phase").toDouble();
415 current_atom.chirp_atom.chirp = atom.attribute(
"chirp").toDouble();
417 if(atom.hasChildNodes())
419 QString sample_string = atom.firstChild().toElement().attribute(
"samples");
420 QStringList sample_list = sample_string.split(
":", QString::SkipEmptyParts);
421 current_atom.atom_samples = VectorXd::Zero(sample_list.length());
422 for(qint32 i = 0; i < sample_list.length(); i++)
423 current_atom.atom_samples[i] = sample_list.at(i).toDouble();
425 current_dict.atoms.append(current_atom);
426 atom = atom.nextSiblingElement(
"ATOM");
431 current_dict.type = FORMULAATOM;
433 while(!atom.isNull())
437 current_atom.id = atom.attribute(
"id").toInt();
438 current_atom.formula_atom.a = atom.attribute(
"a").toDouble();
439 current_atom.formula_atom.b = atom.attribute(
"b").toDouble();
440 current_atom.formula_atom.c = atom.attribute(
"c").toDouble();
441 current_atom.formula_atom.d = atom.attribute(
"d").toDouble();
442 current_atom.formula_atom.e = atom.attribute(
"e").toDouble();
443 current_atom.formula_atom.f = atom.attribute(
"f").toDouble();
444 current_atom.formula_atom.g = atom.attribute(
"g").toDouble();
445 current_atom.formula_atom.h = atom.attribute(
"h").toDouble();
447 if(atom.hasChildNodes())
449 QString sample_string = atom.firstChild().toElement().attribute(
"samples");
450 QStringList sample_list = sample_string.split(
":", QString::SkipEmptyParts);
451 current_atom.atom_samples = VectorXd::Zero(sample_list.length());
452 for(qint32 i = 0; i < sample_list.length(); i++)
453 current_atom.atom_samples[i] = sample_list.at(i).toDouble();
455 current_dict.atoms.append(current_atom);
456 atom = atom.nextSiblingElement(
"ATOM");
465 QString FixDictMp::create_display_text(
FixDictAtom global_best_matching)
467 QString display_text =
"";
468 if(global_best_matching.type == GABORATOM)
470 display_text = QString(
"Gaboratom: scale: %0, translation: %1, modulation: %2, phase: %3")
471 .arg(QString::number(global_best_matching.gabor_atom.scale,
'f', 2))
472 .arg(QString::number(global_best_matching.translation,
'f', 2))
473 .arg(QString::number(global_best_matching.gabor_atom.modulation,
'f', 2))
474 .arg(QString::number(global_best_matching.gabor_atom.phase,
'f', 2));
476 else if(global_best_matching.type == CHIRPATOM)
478 display_text = QString(
"Chripatom: scale: %0, translation: %1, modulation: %2, phase: %3, chirp: %4")
479 .arg(QString::number(global_best_matching.chirp_atom.scale,
'f', 2))
480 .arg(QString::number(global_best_matching.translation,
'f', 2))
481 .arg(QString::number(global_best_matching.chirp_atom.modulation,
'f', 2))
482 .arg(QString::number(global_best_matching.chirp_atom.phase,
'f', 2))
483 .arg(QString::number(global_best_matching.chirp_atom.chirp,
'f', 2));
485 else if(global_best_matching.type == FORMULAATOM)
487 display_text = QString(
"%0: transl: %1 a: %2, b: %3 c: %4, d: %5, e: %6, f: %7, g: %8, h: %9")
488 .arg(global_best_matching.atom_formula)
489 .arg(QString::number(global_best_matching.translation,
'f', 2))
490 .arg(QString::number(global_best_matching.formula_atom.a,
'f', 2))
491 .arg(QString::number(global_best_matching.formula_atom.b,
'f', 2))
492 .arg(QString::number(global_best_matching.formula_atom.c,
'f', 2))
493 .arg(QString::number(global_best_matching.formula_atom.d,
'f', 2))
494 .arg(QString::number(global_best_matching.formula_atom.e,
'f', 2))
495 .arg(QString::number(global_best_matching.formula_atom.f,
'f', 2))
496 .arg(QString::number(global_best_matching.formula_atom.g,
'f', 2))
497 .arg(QString::number(global_best_matching.formula_atom.h,
'f', 2));
505 void FixDictMp::recieve_input(MatrixXd signal, qint32 max_iterations, qreal epsilon, qint32 boost, QString path, qreal delta)
507 matching_pursuit(signal, max_iterations, epsilon, boost, path, delta);
512 qint32 Dictionary::atom_count()
514 return atoms.length();
519 void Dictionary::clear()
522 this->atom_formula =
"";
523 this->sample_count = 0;
FixDictAtom used in fix dict MP Algorithm.