1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
|
/***********************************************/
/**
* @file graceSstResidualAnalysis.cpp
*
* @brief Multiresolution analysis of GRACE SST residuals.
*
* @author Saniya Behzadpour
* @date 2018-07-15
*/
/***********************************************/
// Latex documentation
#define DOCSTRING docstring
static const char *docstring = R"(
This program applies the Multi-Resolution Analysis (MRA) using
Discrete Wavelet Transform (DWT) to the monthly GRACE SST post-fit residuals.
First, the residuals are transferred into wavelet domain by applying an 8 level
Daubechies wavelet transform (default).
In the next step, detail coefficients are merged into three major groups
due to their approximate frequency subbands:
\begin{itemize}
\item Low scale details, corresponding to the frequency band above 10 mHz;
\item Intermediate scale details, corresponding to the approximate frequency
range above 3 mHz up to 10 mHz;
\item High scale details, corresponding to the approximate frequency range
above 0.5 mHz up to 10 mHz.
\end{itemize}
In the last step, each group is reconstructed back into time domain.
)";
/***********************************************/
#include "base/wavelets.h"
#include "files/fileMatrix.h"
#include "files/fileInstrument.h"
#include "parallel/parallel.h"
#include "programs/program.h"
/***** CLASS ***********************************/
/** @brief Multiresolution analysis of GRACE SST residuals.
* @ingroup programsGroup */
class GraceSstResidualAnalysis
{
public:
void run(Config &config, Parallel::CommunicatorPtr comm);
};
GROOPS_REGISTER_PROGRAM(GraceSstResidualAnalysis, SINGLEPROCESS, "Multiresolution analysis of GRACE SST residuals.", Grace)
/***********************************************/
void GraceSstResidualAnalysis::run(Config &config, Parallel::CommunicatorPtr /*comm*/)
{
try
{
FileName fileNameOutHigh, fileNameOutMid, fileNameOutLow;
FileName fileNameInResiduals, fileNameInWavelet;
UInt level = 8;
readConfig(config, "outputfileInstrumentHighScale", fileNameOutHigh, Config::MUSTSET, "", "High scale details");
readConfig(config, "outputfileInstrumentMidScale", fileNameOutMid, Config::MUSTSET, "", "Intermediate scale details");
readConfig(config, "outputfileInstrumentLowScale", fileNameOutLow, Config::MUSTSET, "", "Low scale details");
readConfig(config, "inputfileInstrument", fileNameInResiduals, Config::MUSTSET, "", "GRACE SST Residuals");
readConfig(config, "inputfileWavelet", fileNameInWavelet, Config::MUSTSET, "{groopsDataDir}/wavelets/db20.txt", "wavelet coefficients");
if(isCreateSchema(config)) return;
// =======================
logStatus<<"read GRACE SST residuals "<<"<"<<fileNameInResiduals<<">"<<Log::endl;
SatelliteTrackingArc sstArc = InstrumentFile::read(fileNameInResiduals);
UInt posCount = sstArc.size();
std::vector<Double> signal(posCount);
for(UInt i=0; i<posCount; i++)
signal.at(i) = sstArc.at(i).rangeRate;
logStatus<<"read wavelet filters "<<"<"<<fileNameInWavelet<<">"<<Log::endl;
Vector wl;
readFileMatrix(fileNameInWavelet, wl);
std::vector<UInt> length;
std::vector<UInt> flag, flag2;
std::vector<Double> dwtCoeff,dwt_a08,idwt_a08,idwt_out;
std::vector<Double> dwt_d03, dwt_d02, dwt_d01;
std::vector<Double> idwt_d03,idwt_d02,idwt_d01;
//perform 8-Level DWT
Wavelets::discreteWaveletTransform(signal, wl, level, dwtCoeff,flag, length);
//compute cumulative flag vector
flag2.resize(length.size());
for (UInt i=0; i<length.size();i++)
flag2.at(i) = length.at(i);
for (UInt i=1; i<flag2.size();i++)
flag2.at(i) += flag2.at(i-1);
//merge detail coefficients into three major groups
dwt_a08.resize(dwtCoeff.size()); dwt_d03.resize(dwtCoeff.size());
dwt_d02.resize(dwtCoeff.size()); dwt_d01.resize(dwtCoeff.size());
for (UInt i = 0; i < dwtCoeff.size();i++)
if (i < flag2.at(0))
dwt_a08.at(i) = dwtCoeff.at(i);
else if (i < flag2.at(3))
dwt_d03.at(i) = dwtCoeff.at(i);
else if (i < flag2.at(5))
dwt_d02.at(i) = dwtCoeff.at(i);
else
dwt_d01.at(i) = dwtCoeff.at(i);
//apply inverse DWT, from wavelet domain to time domain
Wavelets::inverseDiscreteWaveletTransform(dwt_d03, wl, flag, idwt_d03, length);
Wavelets::inverseDiscreteWaveletTransform(dwt_d02, wl, flag, idwt_d02, length);
Wavelets::inverseDiscreteWaveletTransform(dwt_d01, wl, flag, idwt_d01, length);
SatelliteTrackingArc d03, d02, d01;
for(UInt i=0; i<posCount; i++)
{
SatelliteTrackingEpoch epoch;
epoch.time = sstArc.at(i).time;
epoch.range = epoch.rangeRate = epoch.rangeAcceleration = idwt_d03.at(i);
d03.push_back(epoch);
}
for(UInt i=0; i<posCount; i++)
{
SatelliteTrackingEpoch epoch;
epoch.time = sstArc.at(i).time;
epoch.range = epoch.rangeRate = epoch.rangeAcceleration = idwt_d02.at(i);
d02.push_back(epoch);
}
for(UInt i=0; i<posCount; i++)
{
SatelliteTrackingEpoch epoch;
epoch.time = sstArc.at(i).time;
epoch.range = epoch.rangeRate = epoch.rangeAcceleration = idwt_d01.at(i);
d01.push_back(epoch);
}
logStatus<<"write decomposed signals"<<Log::endl;
logStatus<<"<"<<fileNameOutHigh<<">"<<Log::endl;
InstrumentFile::write(fileNameOutHigh, d03);
logStatus<<"<"<<fileNameOutMid<<">"<<Log::endl;
InstrumentFile::write(fileNameOutMid, d02);
logStatus<<"<"<<fileNameOutLow<<">"<<Log::endl;
InstrumentFile::write(fileNameOutLow, d01);
}
catch(std::exception &e)
{
GROOPS_RETHROW(e)
}
}
/*************************************************/
|
