Line data Source code
1 : /** \file averagePeriodogram.hpp
2 : * \brief A class to manage calculation of periodograms from time series data.
3 : *
4 : * \author Jared R. Males (jaredmales@gmail.com)
5 : *
6 : * \ingroup signal_processing_files
7 : *
8 : */
9 :
10 : //***********************************************************************//
11 : // Copyright 2020 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 averagePeriodogram_hpp
30 : #define averagePeriodogram_hpp
31 :
32 : #include "psdUtils.hpp"
33 :
34 : namespace mx
35 : {
36 : namespace sigproc
37 : {
38 :
39 : /// Calculate the average periodogram of a time-series.
40 : /**
41 : * Implements the overlapped average periodogram, cf. pp 843 to 845 of \cite 10.5555.1795494. A.k.a. Welch's method
42 : * (https://en.wikipedia.org/wiki/Welch's_method).
43 : *
44 : * This produces
45 : * the variance normalized 1-sided PSD estimate of the periodogram.
46 : *
47 : * Optionally includes a window (by default no window is used, which is equivalent to the rectangle window).
48 : *
49 : * Optionally includes zero-padding, which is applied after windowing. This can be used to increase frequency
50 : * resolution while still allowing overlapping.
51 : *
52 : * Can also be used to calculate the unaveraged non-overlapped periodogram of the entire time series, with or without
53 : * windowing and zero-padding.
54 : *
55 : * Example:
56 : \code
57 : typedef float realT;
58 :
59 : mx::fftwEnvironment<realT> fftwEnv; //for efficient fft planning
60 :
61 : std::vector<realT> ts = getTimeSeries(); //Some function to populate ts.
62 :
63 : math::vectorMeanSub(ts); //The time series should be overall mean subtracted.
64 :
65 : realT dt = 0.1; //The sampling of ts is 0.1 seconds.
66 :
67 : averagePeriodogram<realT> avgPgram(2.0/dt, 4.0/dt, dt); //This sets the averaging length to 2 seconds (20 samples),
68 : zero padding to 4 seconds, with a 1 second overlap (set automatically)
69 :
70 : avgPgram.win(window::hann); //Set the Hann window
71 :
72 : std::vector<realT> pgram = avgPgram(ts); //Calculate the half-overlapped periodogram estimate for ts.
73 :
74 : //Now output the periodogram estimate. This will give two columns, the first frequency, the second power at that
75 : frequency. for(size_t n=0; n< avgPgram.size(); ++n)
76 : {
77 : std::cout << avgPgram[n] << " " << pgram[n] << "\n";
78 : }
79 :
80 : \endcode
81 : *
82 : * \ingroup psds
83 : */
84 : template <typename realT>
85 : class averagePeriodogram
86 : {
87 :
88 : protected:
89 : /** \name Configuration Parameters
90 : *
91 : *@{
92 : */
93 :
94 : size_t m_avgLen{ 0 }; ///< The number of samples in each periodogram estimate.
95 : size_t m_overlap{ 0 }; ///< The number of samples by which to overlap. This should almost always be 0.5*m_avgLen.
96 : ///< Set 0 for the non-overlapped case.
97 :
98 : size_t m_padLen{ 0 }; ///< The size of the periodogram estimate with zero-padding
99 :
100 : protected:
101 : realT m_dt{ 1 }; ///< The time sampling. Only used for normalization and calculation of the frequency scale.
102 :
103 : std::vector<realT> m_win; ///< The window function. By default this is empty, which is equivalent to setting it to
104 : ///< the rectangular window.
105 :
106 : ///@}
107 :
108 : size_t m_size{ 0 }; ///< The size of the periodogram vector, calculated as m_avgLen/2 + 1;
109 :
110 : int m_nOver{ 0 }; ///< The number of overlapping segments. Calculated from m_avgLen and m_overlap;
111 :
112 : realT m_df{ 1 }; ///< The frequency sampling. This is used only for normalization and frequency scale output.
113 :
114 : math::ft::fftT<realT, std::complex<realT>, 1, 0> m_fft;
115 :
116 : realT *m_tsWork{ nullptr };
117 :
118 : std::complex<realT> *m_fftWork{ nullptr };
119 : size_t m_fftWorkSize{ 0 };
120 :
121 : public:
122 : /// C'tor which sets up the optimum overlapped periodogram of the timeseries.
123 : /** Sets m_overlap = 0.5*m_avgLen. If you desire the non-overlapped case use
124 : * the alternate constructor:
125 : * \code
126 : averagePeriogram p( avgLen, 0, dt);
127 : \endcode
128 : */
129 : explicit averagePeriodogram( size_t avgLen /**< [in] The length of averaging in samples.*/ );
130 :
131 : /// C'tor which sets up the optimum overlapped periodogram of the timeseries with zero-padding
132 : /** Sets m_overlap = 0.5*m_avgLen. If you desire the non-overlapped case use
133 : * the alternate constructor:
134 : * \code
135 : averagePeriogram p( avgLen, 0, padLen, dt);
136 : \endcode
137 : */
138 : explicit averagePeriodogram(
139 : size_t avgLen, ///< [in] The length of averaging in samples.
140 : size_t padLen ///< [in] [optional] the length of the padded sample. If 0 or omitted then no padding used.
141 : );
142 :
143 : /// C'tor which sets up the optimum overlapped periodogram of the timeseries with zero-padding and sets the
144 : /// sampling.
145 : /** Sets m_overlap = 0.5*m_avgLen. If you desire the non-overlapped case use
146 : * the alternate constructor:
147 : * \code
148 : averagePeriogram p( avgLen, 0, dt);
149 : \endcode
150 : */
151 : averagePeriodogram( size_t avgLen, ///< [in] The length of averaging in samples.
152 : size_t padLen, ///< [in] the length of the padded sample. Set to 0 for no padding.
153 : realT dt ///< [in] the sampling interval of the time-series
154 : );
155 :
156 : /// C'tor setting up an arbitrary overlap.
157 : /** Set olap to 0 for the unoverlapped case.
158 : */
159 : averagePeriodogram( size_t avgLen, ///< [in] The number of samples in each periodgoram estimate.
160 : size_t padLen, ///< [in] the length of the padded sample. Set to 0 for no padding.
161 : size_t olap, ///< [in] The number of samples by which to overlap. This should almost always be
162 : ///< 0.5*m_avgLen. Set 0 for the non-overlapped case.
163 : realT dt ///< [in] the sampling interval of the time-series
164 : );
165 :
166 : /// D'tor, frees all working memory.
167 : ~averagePeriodogram();
168 :
169 : /// Set the sampling interval of the time-series
170 : /** This also sets the frequency scale of the output.
171 : */
172 : void dt( realT ndt );
173 :
174 : /// Get the sampling interval of the time-series
175 : /** \returns m_dt
176 : */
177 : realT dt();
178 :
179 : /// Get the frequency interval of the periodogram
180 : /** \returns m_df
181 : */
182 : realT df();
183 :
184 : /// Resize the periodogram working memory, setting up the 1/2-overlapped optimum case.
185 : /** This sets the overlap to 0.5*avgLen.
186 : *
187 : * Also performs fft planning.
188 : *
189 : * \returns 0 on success
190 : * \returns -1 on error
191 : */
192 : int resize( size_t avgLen /**< [in] The number of samples in each periodgoram estimate. */ );
193 :
194 : /// Resize the periodogram working memory.
195 : /** Also performs fft planning.
196 : *
197 : * \returns 0 on success
198 : * \returns -1 on error
199 : */
200 : int resize( size_t avgLen, ///< [in] The number of samples in each periodgoram estimate.
201 : size_t padLen, ///< [in] the length of the sample after zero-padding.
202 : size_t olap ///< [in] The number of samples by which to overlap. This should almost always be
203 : ///< 0.5*m_avgLen. Set 0 for the non-overlapped case.
204 : );
205 :
206 : /// Get a reference to the window vector.
207 : /** This allows population of the window with any arbitrary function. You should call
208 : * resizeWin() before using this to make sure that the vector has the correct length.
209 : *
210 : * \returns a reference to m_win.
211 : */
212 : std::vector<realT> &win();
213 :
214 : /// Set the window vector using a function.
215 : /** This will resize m_win, and then call the function with m_win as the argument.
216 : * For example, to use the Hann window from the \ref signal_windows1D "1-D window functions":
217 : * \code
218 : avgPgram.win( window::hann );
219 : \endcode
220 : * which will set the Hann window.
221 : */
222 : void win( void(*winFunc)(std::vector<realT> &) /**<[in] pointer to a function which takes a pre-allocated vector and fills it in with a window function*/);
223 :
224 : /// Resize the window and, unless not desired, initialize to the rectangular window by setting all 1s.
225 : /** If `setRect=false` then the window is not initialized.
226 : */
227 : void resizeWin( bool setRect = true /**< [in] [optional] if false, then the window is not initialized to 1 */ );
228 :
229 : /// Calculate the periodogram for a time-series.
230 : /**
231 : *
232 : */
233 : void operator()( realT *pgram, ///< [out] a pre-allocated to size() array which will be populated with the
234 : ///< periodogram of the time-series
235 : const realT *ts, ///< [in] the time-series
236 : size_t sz ///< [in] the length of the time-series array
237 : );
238 :
239 : /// Calculate the periodogram for a time-series.
240 : /** \overload
241 : *
242 : */
243 : void operator()( std::vector<realT> &pgram, ///< [out] a vector which will be allocated and populated with the
244 : ///< periodogram of the time-series
245 : const std::vector<realT> &ts ///< [in] the time-series
246 : );
247 :
248 : /// Calculate the periodogram for a time-series.
249 : /** \overload
250 : *
251 : * \returns the periodogram as a vector.
252 : */
253 : std::vector<realT> operator()( std::vector<realT> &ts /**< [in] the time-series*/ );
254 :
255 : /// Return the size of the periodogram.
256 : /**
257 : *
258 : * \returns the size of the periodogram estimate for the current setup.
259 : */
260 : size_t size();
261 :
262 : /// Get the frequency at a given point in the periodogram.
263 : /**
264 : * \returns the frequency at point n in the periodogram.
265 : */
266 : realT operator[]( size_t n /**<[in] the point in the periodogram at which the frequency is desired*/ );
267 : };
268 :
269 : template <typename realT>
270 0 : averagePeriodogram<realT>::averagePeriodogram( size_t avgLen )
271 : {
272 0 : resize( avgLen, 0, 0.5 * avgLen );
273 0 : dt( 1 );
274 0 : }
275 :
276 : template <typename realT>
277 0 : averagePeriodogram<realT>::averagePeriodogram( size_t avgLen, size_t padLen )
278 : {
279 0 : resize( avgLen, padLen, 0.5 * avgLen );
280 0 : dt( 1 );
281 0 : }
282 :
283 : template <typename realT>
284 0 : averagePeriodogram<realT>::averagePeriodogram( size_t avgLen, size_t padLen, realT ndt )
285 : {
286 0 : resize( avgLen, padLen, 0.5 * avgLen );
287 0 : dt( ndt );
288 0 : }
289 :
290 : template <typename realT>
291 0 : averagePeriodogram<realT>::averagePeriodogram( size_t avgLen, size_t padLen, size_t olap, realT ndt )
292 : {
293 0 : resize( avgLen, padLen, olap );
294 0 : dt( ndt );
295 0 : }
296 :
297 : template <typename realT>
298 0 : averagePeriodogram<realT>::~averagePeriodogram()
299 : {
300 0 : if( m_tsWork )
301 : {
302 0 : fftw_free( m_tsWork );
303 : }
304 :
305 0 : if( m_fftWork )
306 : {
307 0 : fftw_free( m_fftWork );
308 : }
309 0 : }
310 :
311 : template <typename realT>
312 0 : void averagePeriodogram<realT>::dt( realT ndt )
313 : {
314 0 : m_dt = ndt;
315 0 : m_df = 1.0 / ( m_padLen * m_dt );
316 0 : }
317 :
318 : template <typename realT>
319 0 : realT averagePeriodogram<realT>::dt()
320 : {
321 0 : return m_dt;
322 : }
323 :
324 : template <typename realT>
325 0 : realT averagePeriodogram<realT>::df()
326 : {
327 0 : return m_df;
328 : }
329 :
330 : template <typename realT>
331 0 : int averagePeriodogram<realT>::resize( size_t avgLen )
332 : {
333 0 : return resize( avgLen, 0, 0.5 * avgLen );
334 : }
335 :
336 : template <typename realT>
337 0 : int averagePeriodogram<realT>::resize( size_t avgLen, size_t padLen, size_t olap )
338 : {
339 0 : m_avgLen = avgLen;
340 0 : m_padLen = padLen;
341 :
342 0 : if( m_padLen < m_avgLen )
343 : {
344 0 : m_padLen = m_avgLen;
345 : }
346 :
347 0 : m_size = m_padLen / 2 + 1;
348 :
349 0 : m_df = 1.0 / ( m_padLen * m_dt );
350 :
351 0 : m_overlap = olap;
352 :
353 0 : m_nOver = ( m_avgLen - m_overlap );
354 :
355 0 : m_fft.plan( m_padLen, math::ft::dir::forward, false );
356 :
357 0 : if( m_tsWork )
358 : {
359 0 : fftw_free( m_tsWork );
360 : }
361 :
362 0 : m_tsWork = math::ft::fftw_malloc<realT>( m_padLen );
363 :
364 0 : if( m_fftWork )
365 : {
366 0 : fftw_free( m_fftWork );
367 : }
368 :
369 0 : m_fftWork = math::ft::fftw_malloc<std::complex<realT>>( ( m_padLen / 2 + 1 ) );
370 :
371 0 : return 0;
372 : }
373 :
374 : template <typename realT>
375 0 : std::vector<realT> &averagePeriodogram<realT>::win()
376 : {
377 0 : return m_win;
378 : }
379 :
380 : template <typename realT>
381 0 : void averagePeriodogram<realT>::win( void ( *winFunc )( std::vector<realT> & ) )
382 : {
383 0 : resizeWin( false );
384 0 : winFunc( m_win );
385 0 : }
386 :
387 : template <typename realT>
388 0 : void averagePeriodogram<realT>::resizeWin( bool setRect )
389 : {
390 0 : if( setRect )
391 : {
392 0 : m_win.resize( m_avgLen, static_cast<realT>( 1 ) );
393 : }
394 : else
395 : {
396 0 : m_win.resize( m_avgLen );
397 : }
398 0 : }
399 :
400 : template <typename realT>
401 0 : void averagePeriodogram<realT>::operator()( realT *pgram, const realT *ts, size_t sz )
402 : {
403 0 : if( m_win.size() > 0 && m_win.size() != m_avgLen )
404 : {
405 0 : std::cerr << "averagePeriodogram: Window size not correct.\n";
406 : }
407 :
408 0 : int Navg = sz / m_nOver;
409 :
410 0 : while( Navg * m_nOver + m_avgLen > sz )
411 0 : --Navg;
412 :
413 0 : if( Navg < 1 )
414 0 : Navg = 1; // Always do at least 1!
415 :
416 0 : for( int i = 0; i < Navg; ++i )
417 : {
418 0 : if( m_win.size() == m_avgLen ) // no if inside the for loop
419 : {
420 0 : for( size_t j = 0; j < m_avgLen; ++j )
421 : {
422 0 : m_tsWork[j] = ts[i * m_nOver + j] * m_win[j];
423 : }
424 : }
425 : else
426 : {
427 0 : for( size_t j = 0; j < m_avgLen; ++j )
428 : {
429 0 : m_tsWork[j] = ts[i * m_nOver + j];
430 : }
431 : }
432 :
433 : // Zero-pad
434 0 : for( size_t j = m_avgLen; j < m_padLen; ++j )
435 : {
436 0 : m_tsWork[j] = 0;
437 : }
438 :
439 0 : m_fft( m_fftWork, m_tsWork );
440 :
441 0 : for( size_t j = 0; j < m_size; ++j )
442 0 : pgram[j] += norm( m_fftWork[j] );
443 : }
444 :
445 0 : realT pgramVar = psdVar1sided( m_df, pgram, m_size );
446 :
447 0 : realT tsVar = mx::math::vectorVariance( ts, sz );
448 :
449 0 : for( size_t j = 0; j < m_size; ++j )
450 0 : pgram[j] *= tsVar / pgramVar;
451 0 : }
452 :
453 : template <typename realT>
454 0 : void averagePeriodogram<realT>::operator()( std::vector<realT> &pgram, const std::vector<realT> &ts )
455 : {
456 0 : pgram.resize( m_size );
457 0 : operator()( pgram.data(), ts.data(), ts.size() );
458 0 : }
459 :
460 : template <typename realT>
461 0 : std::vector<realT> averagePeriodogram<realT>::operator()( std::vector<realT> &ts )
462 : {
463 0 : std::vector<realT> pgram;
464 0 : operator()( pgram, ts );
465 0 : return pgram;
466 0 : }
467 :
468 : template <typename realT>
469 0 : size_t averagePeriodogram<realT>::size()
470 : {
471 0 : return m_size;
472 : }
473 :
474 : template <typename realT>
475 0 : realT averagePeriodogram<realT>::operator[]( size_t n )
476 : {
477 0 : return n * m_df;
478 : }
479 :
480 : } // namespace sigproc
481 : } // namespace mx
482 :
483 : #endif // averagePeriodogram_hpp
|
