log_taper.cc 4.02 KB
Newer Older
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
/*! \file log_taper.cc
 * \brief provide taper for averaging on a logarithmic scale (implementation)
 * 
 * ----------------------------------------------------------------------------
 * 
 * \author Thomas Forbriger
 * \date 12/01/2019
 * 
 * provide taper for averaging on a logarithmic scale (implementation)
 * 
 * Copyright (c) 2019 by Thomas Forbriger (BFO Schiltach) 
 * 
 * ----
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version. 
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program. If not, see <http://www.gnu.org/licenses/>.
 * ----
 *
 * REVISIONS and CHANGES 
 *  - 12/01/2019   V1.0   Thomas Forbriger
 * 
 * ============================================================================
 */
#define PSDXX_LOG_TAPER_CC_VERSION \
  "PSDXX_LOG_TAPER_CC   V1.0"

#include <psdxx/helper.h>
37
#include <psdxx/error.h>
38
#include <psdxx/debug.h>
39
#include <cmath>
40
41
42
43
44

namespace psd {

  namespace helper {

45
46
    LogTaper::LogTaper(const psd::TDseries::Tcoc& f,
                       const unsigned int& i): Msum(0)
47
48
49
    {
      PSDXX_assert(f.size()>2,
                   "at least three values of frequency are required");
50
51
52
      PSDXX_assert((i>=f.f()) && (i<=f.l()),
                   "index is out of range");
      if (i==f.f())
53
54
      {
        this->Mloc=Flower_end;
55
        this->Mli=new psd::helper::LogIndex(f(i), f(i+1), f(i+2));
56
      }
57
      else if (i==f.l())
58
59
      {
        this->Mloc=Fupper_end;
60
        this->Mli=new psd::helper::LogIndex(f(i-2), f(i-1), f(i));
61
62
63
64
      }
      else
      {
        this->Mloc=Fmiddle;
65
        this->Mli=new psd::helper::LogIndex(f(i-1), f(i), f(i+1));
66
      }
67
68
    } // LogTaper::LogTaper(const psd::TDseries::Tcoc& f,
      //                    const unsigned int& i)
69
70
71

  /* ---------------------------------------------------------------------- */

72
    LogTaper::~LogTaper()
73
74
    {
      delete Mli;
75
    } // LogTaper::~LogTaper()
76
77
78

  /* ---------------------------------------------------------------------- */

79
    double LogTaper::operator()(const double& f) 
80
81
    {
      double retval=0;
82
83
84
85
86
87
88
89
90
91
92
93
      double rindex=this->Mli->operator()(f);

      // taper function
      // 0.5*(1.-cos(ri*pi))

      if (this->Mloc==Flower_end)
      {
        // index of center frequency is 0
        // index range is from -1 to +1
        // taper is 1 for frequency index <= 0
        if ((rindex > 0.) && (rindex < 1.))
        {
94
          retval=0.5*(1.+std::cos(rindex*M_PI));
95
96
97
98
99
100
101
102
103
104
105
        }
        else if (rindex <= 0.)
        {
          retval=1.;
        }
      }
      else if (this->Mloc==Fupper_end)
      {
        // index of center frequency is 2
        // index range is from +1 to +3
        // taper is 1 for frequency index >= 2
106
        if ((rindex >= 1.) && (rindex < 2.))
107
        {
108
109
110
111
112
113
          retval=0.5*(1.+std::cos(rindex*M_PI));
          /*
          std::cerr << PSDXX_value(f) << " "
            << PSDXX_value(rindex) << " "
            << PSDXX_value(retval) << std::endl;
            */
114
115
116
117
118
119
120
121
122
123
        }
        else if (rindex >= 2.)
        {
          retval=1.;
        }
      }
      else if (this->Mloc==Fmiddle)
      {
        // index of center frequency is 1
        // index range is from 0 to +2
124
        if ((rindex > -0.) && (rindex < 2.))
125
126
127
128
129
130
131
132
133
134
        {
          retval=0.5*(1.-std::cos(rindex*M_PI));
        }
      }
      else
      {
        PSDXX_abort("illegal mode; this should never happen; "
                    "report this as a bug");
      }

135
136
      this->Msum += retval;
      return(retval);
137
    } // double LogTaper::operator()(const double& f) const
138

139
140
141
142
143
  } // namespace helper

} // namespace psd

/* ----- END OF log_taper.cc ----- */