Commit 9aebaf9d authored by thomas.forbriger's avatar thomas.forbriger
Browse files

croposp [FEATURE]: support selection of nspd mode

parent 55e55d1d
...@@ -73,7 +73,8 @@ int main(int iargc, char* argv[]) ...@@ -73,7 +73,8 @@ int main(int iargc, char* argv[])
" [-pattern p] [-overwrite] [-DEBUG]" "\n" " [-pattern p] [-overwrite] [-DEBUG]" "\n"
" [-log n] [-avgout]\n" " [-log n] [-avgout]\n"
" [-nsegments n] [-divisor n] [-padfactor n] [-overlap f]\n" " [-nsegments n] [-divisor n] [-padfactor n] [-overlap f]\n"
" [-psd f] [-npsd f] [-transfer f] [-phase f] [-coherency f]" "\n" " [-psd f] [-transfer f] [-phase f] [-coherency f]" "\n"
" [-npsd f] [-npsdmode f]\n"
" [-prefixpsd s] [-prefixnpsd s] [-prefixtransfer s]\n" " [-prefixpsd s] [-prefixnpsd s] [-prefixtransfer s]\n"
" [-prefixphase s] [-prefixcoherency s]" "\n" " [-prefixphase s] [-prefixcoherency s]" "\n"
" file [t:sel] [f:format] [n:label] [file [t:s] [f:f] [n:l]] [...]\n" " file [t:sel] [f:format] [n:label] [file [t:s] [f:f] [n:l]] [...]\n"
...@@ -90,7 +91,7 @@ int main(int iargc, char* argv[]) ...@@ -90,7 +91,7 @@ int main(int iargc, char* argv[])
"-trim trim time series to identical begin and length\n" "-trim trim time series to identical begin and length\n"
"-datetolerance t set tolerance as a fraction of sampling interval\n" "-datetolerance t set tolerance as a fraction of sampling interval\n"
" when comparing time series header date values\n" " when comparing time series header date values\n"
"-pattern p set pattern for time series label\n" "-pattern p set pattern for signal label\n"
" %C: channel (WID2 header field)\n" " %C: channel (WID2 header field)\n"
" %S: station (WID2 header field)\n" " %S: station (WID2 header field)\n"
" %I: instype (WID2 header field)\n" " %I: instype (WID2 header field)\n"
...@@ -98,11 +99,24 @@ int main(int iargc, char* argv[]) ...@@ -98,11 +99,24 @@ int main(int iargc, char* argv[])
" %N: label set for file name\n" " %N: label set for file name\n"
" %F: file name\n" " %F: file name\n"
" %NT: number of trace in file\n" " %NT: number of trace in file\n"
"-prefixpsd s set label prefix for psd output to \"s\"\n" "-patpsd s set label pattern for psd output to \"s\"\n"
"-prefixnpsd s set label prefix for incoherent psd output to \"s\"\n" " %1: signal label\n"
"-prefixtransfer s set label prefix for amplitude ratio output to \"s\"\n" "-patnpsd s set label pattern for incoherent psd output to \"s\"\n"
"-prefixphase s set label prefix for phase delay output to \"s\"\n" " %1: label of signal\n"
"-prefixcoherency s set label prefix for coherency output to \"s\"\n" " %2: label of first reference signal\n"
" %3: label of second reference signal\n"
" %M: label of npsd mode\n"
"-pattransfer s set label pattern for amplitude ratio output to \"s\"\n"
" %1: label of first signal\n"
" %2: label of second signal\n"
" %3: label of reference signal\n"
"-patphase s set label pattern for phase delay output to \"s\"\n"
" %1: label of first signal\n"
" %2: label of second signal\n"
" %3: label of reference signal\n"
"-patcoherency s set label pattern for coherency output to \"s\"\n"
" %1: label of first signal\n"
" %2: label of second signal\n"
"\n" "\n"
"computation options:\n" "computation options:\n"
"-log n map averages to logarithmic scale with \"n\"\n" "-log n map averages to logarithmic scale with \"n\"\n"
...@@ -119,6 +133,7 @@ int main(int iargc, char* argv[]) ...@@ -119,6 +133,7 @@ int main(int iargc, char* argv[])
"-padfactor n pad with zeroes to increase length of series by a factor\n" "-padfactor n pad with zeroes to increase length of series by a factor\n"
" of \"n\"\n" " of \"n\"\n"
"-overlap f let segments overlap by a fraction of \"f\"\n" "-overlap f let segments overlap by a fraction of \"f\"\n"
"-npsdmode m select npsd mode (see below)\n"
"\n" "\n"
"output options:\n" "output options:\n"
"-overwrite overwrite existing output files\n" "-overwrite overwrite existing output files\n"
...@@ -181,18 +196,20 @@ int main(int iargc, char* argv[]) ...@@ -181,18 +196,20 @@ int main(int iargc, char* argv[])
{"overlap",arg_yes,"0.5"}, {"overlap",arg_yes,"0.5"},
// 17: compute phase delay // 17: compute phase delay
{"phase",arg_yes,"-"}, {"phase",arg_yes,"-"},
// 18: set psd prefix // 18: set psd label pattern
{"prefixpsd",arg_yes,"PSD"}, {"patpsd",arg_yes,"PSD %1"},
// 19: set npsd prefix // 19: set npsd label pattern
{"prefixnpsd",arg_yes,"NPSD"}, {"patnpsd",arg_yes,"NPSD %1 (ref: %2 %3) %M"},
// 20: set transfer prefix // 20: set transfer label pattern
{"prefixtransfer",arg_yes,"AMP"}, {"pattransfer",arg_yes,"AMP %1/%2 (ref: %3)"},
// 21: set phase prefix // 21: set phase label pattern
{"prefixphase",arg_yes,"PHA"}, {"patphase",arg_yes,"PHA %1-%2 (ref: %3)"},
// 22: set coherency prefix // 22: set coherency label pattern
{"prefixcoherency",arg_yes,"COH"}, {"patcoherency",arg_yes,"COH %1 %2"},
// 23: average upon output // 23: average upon output
{"avgout",arg_no,"-"}, {"avgout",arg_no,"-"},
// 24: average upon output
{"npsdmode",arg_yes,"0"},
{NULL} {NULL}
}; };
...@@ -204,6 +221,14 @@ int main(int iargc, char* argv[]) ...@@ -204,6 +221,14 @@ int main(int iargc, char* argv[])
// define commandline argument modifier keys // define commandline argument modifier keys
static const char* cmdlinekeys[]={tracekey, formatkey, labelkey, 0}; static const char* cmdlinekeys[]={tracekey, formatkey, labelkey, 0};
// nspd mode strings
static const char npsd_label_abs[]="|Pkk-(Plk*Pkj/Plj)|";
static const char npsd_label_abs2[]="|Pkk-|Plk*Pkj/Plj||";
static const char npsd_label_real[]="|Pkk-real(Plk*Pkj/Plj)|";
static const char npsd_label_imag[]="|imag(Plk*Pkj/Plj)|";
static const char* npsd_mode_label[]
={npsd_label_abs, npsd_label_abs2, npsd_label_real, npsd_label_imag, 0};
/* ====================================================================== */ /* ====================================================================== */
/* analyze command line /* analyze command line
* -------------------- * --------------------
...@@ -225,6 +250,12 @@ int main(int iargc, char* argv[]) ...@@ -225,6 +250,12 @@ int main(int iargc, char* argv[])
cerr << usage_text << endl; cerr << usage_text << endl;
cerr << help_text << endl; cerr << help_text << endl;
cerr << endl; cerr << endl;
cerr << "npsd modes:" << endl;
for (unsigned int l=0; l<4; ++l)
{
cerr << l << ": " << npsd_mode_label[l] << endl;
}
cerr << endl << endl;
cerr << "References" << endl cerr << "References" << endl
<< "----------" << endl; << "----------" << endl;
cerr << reference_sleeman_et_al_2006 << endl; cerr << reference_sleeman_et_al_2006 << endl;
...@@ -259,14 +290,17 @@ int main(int iargc, char* argv[]) ...@@ -259,14 +290,17 @@ int main(int iargc, char* argv[])
opt.compute_phase=cmdline.optset(17); opt.compute_phase=cmdline.optset(17);
opt.outfile_phase=cmdline.string_arg(17); opt.outfile_phase=cmdline.string_arg(17);
opt.prefix_psd=cmdline.string_arg(18); opt.pattern_psd=cmdline.string_arg(18);
opt.prefix_npsd=cmdline.string_arg(19); opt.pattern_npsd=cmdline.string_arg(19);
opt.prefix_transfer=cmdline.string_arg(20); opt.pattern_transfer=cmdline.string_arg(20);
opt.prefix_phase=cmdline.string_arg(21); opt.pattern_phase=cmdline.string_arg(21);
opt.prefix_coherency=cmdline.string_arg(22); opt.pattern_coherency=cmdline.string_arg(22);
opt.avg_on_output=cmdline.optset(23); opt.avg_on_output=cmdline.optset(23);
opt.npsd_mode=cmdline.int_arg(24);
TFXX_assert((opt.npsd_mode>=0) && (opt.npsd_mode<4),
"npsd mode out of range");
TFXX_assert(opt.n_per_decade>0, TFXX_assert(opt.n_per_decade>0,
"number of samples per decade must be finite and positive"); "number of samples per decade must be finite and positive");
TFXX_assert(opt.nsegments>0, TFXX_assert(opt.nsegments>0,
...@@ -466,7 +500,9 @@ int main(int iargc, char* argv[]) ...@@ -466,7 +500,9 @@ int main(int iargc, char* argv[])
named_psd.series=psd.data; named_psd.series=psd.data;
} // if (opt.logscale && !opt.avg_on_output) ... else } // if (opt.logscale && !opt.avg_on_output) ... else
named_psd.label=opt.prefix_psd+" "+i_meta->label; named_psd.label=opt.pattern_psd;
named_psd.label=croposp::patsubst(named_psd.label, "%1",
i_meta->label);
PSD_vector.push_back(named_psd); PSD_vector.push_back(named_psd);
if (opt.verbose) if (opt.verbose)
...@@ -619,8 +655,13 @@ int main(int iargc, char* argv[]) ...@@ -619,8 +655,13 @@ int main(int iargc, char* argv[])
psd::TDseries::Tcoc acpsd=psd::abs(CPSD_vector[pairs(k,l)].series); psd::TDseries::Tcoc acpsd=psd::abs(CPSD_vector[pairs(k,l)].series);
coherency_vector[pairs(k,l)].series coherency_vector[pairs(k,l)].series
=acpsd/(sqrtpsd_vector[k].series*sqrtpsd_vector[l].series); =acpsd/(sqrtpsd_vector[k].series*sqrtpsd_vector[l].series);
coherency_vector[pairs(k,l)].label= coherency_vector[pairs(k,l)].label=opt.pattern_coherency;
opt.prefix_coherency+" "+CPSD_vector[pairs(k,l)].label; coherency_vector[pairs(k,l)].label
=croposp::patsubst(coherency_vector[pairs(k,l)].label, "%1",
vector_of_metadata[k].label);
coherency_vector[pairs(k,l)].label
=croposp::patsubst(coherency_vector[pairs(k,l)].label, "%2",
vector_of_metadata[l].label);
if (opt.verbose) if (opt.verbose)
{ {
...@@ -674,10 +715,16 @@ int main(int iargc, char* argv[]) ...@@ -674,10 +715,16 @@ int main(int iargc, char* argv[])
cpsd2=CPSD_vector[pairs(l,i)].series; cpsd2=CPSD_vector[pairs(l,i)].series;
if (pairs.swap(l,i)) { cpsd2 = psd::conj(cpsd2); } if (pairs.swap(l,i)) { cpsd2 = psd::conj(cpsd2); }
phase_vector[triples(i,k,l)].series=psd::arg(cpsd1/cpsd2); phase_vector[triples(i,k,l)].series=psd::arg(cpsd1/cpsd2);
phase_vector[triples(i,k,l)].label=opt.pattern_phase;
phase_vector[triples(i,k,l)].label
=croposp::patsubst(phase_vector[triples(i,k,l)].label, "%1",
vector_of_metadata[k].label);
phase_vector[triples(i,k,l)].label
=croposp::patsubst(phase_vector[triples(i,k,l)].label, "%2",
vector_of_metadata[l].label);
phase_vector[triples(i,k,l)].label phase_vector[triples(i,k,l)].label
=opt.prefix_phase+" "+vector_of_metadata[k].label =croposp::patsubst(phase_vector[triples(i,k,l)].label, "%3",
+"-"+vector_of_metadata[l].label+" (ref: " vector_of_metadata[i].label);
+vector_of_metadata[i].label+")";
} }
} }
} }
...@@ -725,10 +772,16 @@ int main(int iargc, char* argv[]) ...@@ -725,10 +772,16 @@ int main(int iargc, char* argv[])
cpsd1=psd::abs(CPSD_vector[pairs(k,i)].series); cpsd1=psd::abs(CPSD_vector[pairs(k,i)].series);
cpsd2=psd::abs(CPSD_vector[pairs(l,i)].series); cpsd2=psd::abs(CPSD_vector[pairs(l,i)].series);
transfer_vector[triples(i,k,l)].series=cpsd1/cpsd2; transfer_vector[triples(i,k,l)].series=cpsd1/cpsd2;
transfer_vector[triples(i,k,l)].label=opt.pattern_transfer;
transfer_vector[triples(i,k,l)].label transfer_vector[triples(i,k,l)].label
=opt.prefix_transfer+" "+vector_of_metadata[k].label =croposp::patsubst(transfer_vector[triples(i,k,l)].label, "%1",
+"/"+vector_of_metadata[l].label+" (ref: " vector_of_metadata[k].label);
+vector_of_metadata[i].label+")"; transfer_vector[triples(i,k,l)].label
=croposp::patsubst(transfer_vector[triples(i,k,l)].label, "%2",
vector_of_metadata[l].label);
transfer_vector[triples(i,k,l)].label
=croposp::patsubst(transfer_vector[triples(i,k,l)].label, "%3",
vector_of_metadata[i].label);
} }
} }
} }
...@@ -752,6 +805,8 @@ int main(int iargc, char* argv[]) ...@@ -752,6 +805,8 @@ int main(int iargc, char* argv[])
cout << endl cout << endl
<< "Compute incoherent component." << "Compute incoherent component."
<< endl; << endl;
cout << "selected mode: "
<< npsd_mode_label[opt.npsd_mode] << endl;
} }
TFXX_assert(pairs.size()>2, "requires at least three input time series"); TFXX_assert(pairs.size()>2, "requires at least three input time series");
...@@ -800,18 +855,41 @@ int main(int iargc, char* argv[]) ...@@ -800,18 +855,41 @@ int main(int iargc, char* argv[])
while(iNPSD.valid() && iCPSD.valid() && iPSD.valid()) while(iNPSD.valid() && iCPSD.valid() && iPSD.valid())
{ {
*iNPSD = std::abs(*iPSD - (*iCPSD)); switch (opt.npsd_mode) {
case 0:
*iNPSD = std::abs(*iPSD - (*iCPSD));
break;
case 1:
*iNPSD = std::abs(*iPSD-std::abs(*iCPSD));
break;
case 2:
*iNPSD = std::abs(*iPSD-(iCPSD->real()));
break;
case 3:
*iNPSD = std::abs(iCPSD->imag());
break;
default:
TFXX_abort("illegal npsd mode");
}
++iNPSD; ++iNPSD;
++iCPSD; ++iCPSD;
++iPSD; ++iPSD;
} }
npsd_vector[triples(i,k,l)].series=npsd; npsd_vector[triples(i,k,l)].series=npsd;
npsd_vector[triples(i,k,l)].label=opt.pattern_npsd;
npsd_vector[triples(i,k,l)].label
=croposp::patsubst(npsd_vector[triples(i,k,l)].label, "%1",
vector_of_metadata[i].label);
npsd_vector[triples(i,k,l)].label
=croposp::patsubst(npsd_vector[triples(i,k,l)].label, "%2",
vector_of_metadata[k].label);
npsd_vector[triples(i,k,l)].label
=croposp::patsubst(npsd_vector[triples(i,k,l)].label, "%3",
vector_of_metadata[l].label);
npsd_vector[triples(i,k,l)].label npsd_vector[triples(i,k,l)].label
=opt.prefix_npsd+" " =croposp::patsubst(npsd_vector[triples(i,k,l)].label, "%M",
+vector_of_metadata[i].label+" (ref: " npsd_mode_label[opt.npsd_mode]);
+vector_of_metadata[k].label
+" "+vector_of_metadata[l].label+")";
} }
} }
} }
......
...@@ -32,7 +32,7 @@ ...@@ -32,7 +32,7 @@
* ============================================================================ * ============================================================================
*/ */
#define CROPOSP_VERSION \ #define CROPOSP_VERSION \
"CROPOSP V1.1 Cross power spectral density" "CROPOSP V1.2 Cross power spectral density"
// include guard // include guard
#ifndef CROPOSP_H_VERSION #ifndef CROPOSP_H_VERSION
...@@ -63,8 +63,9 @@ namespace croposp { ...@@ -63,8 +63,9 @@ namespace croposp {
bool compute_phase; bool compute_phase;
std::string outfile_psd, outfile_npsd, outfile_coherency, outfile_transfer; std::string outfile_psd, outfile_npsd, outfile_coherency, outfile_transfer;
std::string outfile_phase; std::string outfile_phase;
std::string prefix_psd, prefix_npsd, prefix_transfer; std::string pattern_psd, pattern_npsd, pattern_transfer;
std::string prefix_phase, prefix_coherency; std::string pattern_phase, pattern_coherency;
unsigned int npsd_mode;
}; // struct Options }; // struct Options
// a struct to hold meta data // a struct to hold meta data
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment