Line data Source code
1 : /** \file autocorrelation.hpp
2 : * \brief Tools for working with autocorrelations
3 : *
4 : * \author Jared R. Males (jaredmales@gmail.com)
5 : *
6 : * \ingroup signal_processing_files
7 : *
8 : */
9 :
10 : //***********************************************************************//
11 : // Copyright 2015, 2016, 2017 Jared R. Males (jaredmales@gmail.com)
12 : //
13 : // This file is part of mxlib.
14 : //
15 : // mxlib is free software: you can redistribute it and/or modify
16 : // it under the terms of the GNU General Public License as published by
17 : // the Free Software Foundation, either version 3 of the License, or
18 : // (at your option) any later version.
19 : //
20 : // mxlib is distributed in the hope that it will be useful,
21 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
22 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23 : // GNU General Public License for more details.
24 : //
25 : // You should have received a copy of the GNU General Public License
26 : // along with mxlib. If not, see <http://www.gnu.org/licenses/>.
27 : //***********************************************************************//
28 :
29 : #ifndef autocorrelation_hpp
30 : #define autocorrelation_hpp
31 :
32 : #include "../math/ft/fftT.hpp"
33 :
34 : namespace mx
35 : {
36 :
37 : namespace sigproc
38 : {
39 :
40 : /// Calculate the autocorrelation of a time-series
41 : /**
42 : * \tparam T is the real type of the data and autocorrelation.
43 : *
44 : * \ingroup signal_processing
45 : */
46 : template <typename T>
47 : void autocorrelation(
48 : T *ac, ///< [out] is the pre-allocated array of length Nac which will contain the autocorrelation on return
49 : size_t Nac, ///< [in] is the length of ac
50 : T *sig, ///< [in] is the input time-series (signal)
51 : size_t Nsig ///< [in] is the length of the input time-series
52 : )
53 : {
54 :
55 : #pragma omp parallel for
56 : for( int i = 0; i < Nac; ++i )
57 : {
58 : ac[i] = 0;
59 : for( int j = 0; j < Nsig - i; ++j )
60 : {
61 : ac[i] += sig[j] * sig[j + i];
62 : }
63 : }
64 :
65 : T norm = ac[0];
66 : for( int i = 0; i < Nac; ++i )
67 : {
68 : ac[i] /= norm;
69 : }
70 : }
71 :
72 : /// Calculate the autocorrelation of a time-series
73 : /**
74 : * \tparam T is the real type of the data and autocorrelation.
75 : *
76 : * \todo this should probably re-allocate if ac.size() != sig.size(), or at least <
77 : *
78 : * \ingroup signal_processing
79 : */
80 : template <typename T>
81 : void autocorrelation( std::vector<T> &ac, ///< [out] will contain the autocorrelation on return. If ac.size()==0 then
82 : ///< it is resized to sig.size().
83 : std::vector<T> &sig ///< [in] is the input time-series (signal)
84 : )
85 : {
86 : if( ac.size() == 0 )
87 : {
88 : ac.resize( sig.size() );
89 : }
90 :
91 : autocorrelation( ac.data(), ac.size(), sig.data(), sig.size() );
92 : }
93 :
94 : /// Functor for calculating the autocorrelation given a PSD
95 : /** Stores the fftT object and related working memory so that
96 : * repeated calls do not re-allocate or re-plan the FFT.
97 : *
98 : * \tparam T is the real type of the PSD and resultig A.C.
99 : *
100 : * \ingroup signal_processing
101 : */
102 : template <typename T>
103 : struct autocorrelationFromPSD
104 : {
105 : std::vector<std::complex<T>> fftOut;
106 : std::vector<std::complex<T>> fftIn;
107 :
108 : math::ft::fftT<std::complex<T>, std::complex<T>, 1, 0> fft;
109 :
110 : /// Calculate the A.C. as the inverse FFT of the PSD
111 : /** This calculates the circular autocorrelation from the PSD.
112 : *
113 : */
114 : void
115 121 : operator()( T *ac, ///< [out] pre-allocated array, on output contains the first Nac points of the autocorrelation
116 : size_t Nac, ///< [in] the allocated size of ac.
117 : T *psd, ///< [in] the 2-sided FFT storage order PSD
118 : size_t Npsd ///< [in] the number of points in the PSD
119 : )
120 : {
121 121 : fft.plan( Npsd, math::ft::dir::forward );
122 :
123 121 : fftOut.resize( Npsd );
124 121 : fftIn.resize( Npsd );
125 :
126 15609 : for( size_t i = 0; i < Npsd; ++i )
127 15488 : fftIn[i] = psd[i];
128 :
129 121 : fft( fftOut.data(), fftIn.data() );
130 :
131 121 : T norm = fftOut[0].real();
132 15609 : for( size_t i = 0; i < Npsd && i < Nac; ++i )
133 15488 : ac[i] = fftOut[i].real() / norm;
134 121 : }
135 :
136 : /// Calculate the A.C. as the inverse FFT of the PSD
137 : /** This calculates the circular autocorrelation from the PSD.
138 : *
139 : */
140 121 : void operator()( std::vector<T> &ac, ///< [out] On output contains the autocorrelation. If ac.size()==0, it is
141 : ///< allocated to psd.size().
142 : std::vector<T> &psd ///< [in] the 2-sided FFT storage order PSD
143 : )
144 : {
145 121 : if( ac.size() == 0 )
146 0 : ac.resize( psd.size() );
147 121 : operator()( ac.data(), ac.size(), psd.data(), psd.size() );
148 121 : }
149 : };
150 :
151 : } // namespace sigproc
152 : } // namespace mx
153 :
154 : #endif // autocorrelation_hpp
|
