log_taper.cc 3.79 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 <cmath>
39
40
41
42
43

namespace psd {

  namespace helper {

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
    Taper::Taper(const psd::TDseries::Tcoc& f,
                 const unsigned int& i): Msum(0)
    {
      PSDXX_assert(f.size()>2,
                   "at least three values of frequency are required");
      if (i==0)
      {
        this->Mloc=Flower_end;
        this->Mli=new psd::helper::LogIndex(f(f.f()), f(f.f()+1), f(f.f()+2));
      }
      else if (i==(f.size()-1))
      {
        this->Mloc=Fupper_end;
        this->Mli=new psd::helper::LogIndex(f(f.l()-2), f(f.l()-1), f(f.l()));
      }
      else
      {
        this->Mloc=Fmiddle;
        this->Mli=new psd::helper::LogIndex(f(f.f()+i-1), 
                                            f(f.f()+i), f(f.f()+i+1));
      }
    } // Taper::Taper(const psd::TDseries::Tcoc& f,
      //              const unsigned int& i)

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

    Taper::~Taper()
    {
      delete Mli;
    } // Taper::~Taper()

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

    double Taper::operator()(const double& f) 
    {
      double retval=0;
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
      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.))
        {
          retval=0.5*(1.-std::sin(rindex*M_PI));
        }
        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
        if ((rindex < 2.) && (rindex > 1.))
        {
          retval=0.5*(1.-std::sin(rindex*M_PI));
        }
        else if (rindex >= 2.)
        {
          retval=1.;
        }
      }
      else if (this->Mloc==Fmiddle)
      {
        // index of center frequency is 1
        // index range is from 0 to +2
        if ((rindex > -2.) && (rindex < 2.))
        {
          retval=0.5*(1.-std::cos(rindex*M_PI));
        }
      }
      else
      {
        PSDXX_abort("illegal mode; this should never happen; "
                    "report this as a bug");
      }

128
129
130
131
      this->Msum += retval;
      return(retval);
    } // double Taper::operator()(const double& f) const

132
133
134
135
136
  } // namespace helper

} // namespace psd

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