Line data Source code
1 : /** \file ftUtils.hpp
2 : * \brief Fourier Transform Utilities
3 : * \ingroup ft_files
4 : * \author Jared R. Males (jaredmales@gmail.com)
5 : *
6 : */
7 :
8 : //***********************************************************************//
9 : // Copyright 2025 Jared R. Males (jaredmales@gmail.com)
10 : //
11 : // This file is part of mxlib.
12 : //
13 : // mxlib is free software: you can redistribute it and/or modify
14 : // it under the terms of the GNU General Public License as published by
15 : // the Free Software Foundation, either version 3 of the License, or
16 : // (at your option) any later version.
17 : //
18 : // mxlib is distributed in the hope that it will be useful,
19 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
20 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 : // GNU General Public License for more details.
22 : //
23 : // You should have received a copy of the GNU General Public License
24 : // along with mxlib. If not, see <http://www.gnu.org/licenses/>.
25 : //***********************************************************************//
26 :
27 : #ifndef ftUtils_hpp
28 : #define ftUtils_hpp
29 :
30 : namespace mx
31 : {
32 : namespace math
33 : {
34 : namespace ft
35 : {
36 :
37 : /// Add padding to the output of a real-to-complex (R2C) FFT
38 : /**
39 : * The output array should be pre-sized according to:
40 : * \verbatim
41 : * aout.cols() == ain.cols() + padCols
42 : * aout.rows() == ain.rows() + padRows
43 : * \endverbatim
44 : * where `padRows` and `padCols are the desired amount of padding. \p aout will be in
45 : * FFT order and suitable for input to the inverse transform (complext-to-real).
46 : *
47 : * \ingroup ft
48 : */
49 : template <class eigenArrOutT, class eigenArrInT>
50 24 : void padR2CFFTOutput( eigenArrOutT &aout, /**< [out] The padded array. Must be pre-allocated. The pad region will
51 : not be set to zero. Size of the pad is determined by the
52 : difference in size between \p aout and \p ain*/
53 : eigenArrInT &ain /**< [in] The array to be paddded, in real-to-complex FFTW storage order.*/ )
54 : {
55 24 : int ci = 0.5 * ain.cols() + 1;
56 24 : int cf = ain.cols() - ci;
57 :
58 24 : aout.block( 0, 0, ain.rows(), ci ) = ain.block( 0, 0, ain.rows(), ci );
59 24 : aout.block( 0, aout.cols() - cf, ain.rows(), cf ) = ain.block( 0, ain.cols() - cf, ain.rows(), cf );
60 24 : }
61 :
62 : /// Add padding to the output of a complex-to-complex (C2C) FFT
63 : /** The output array should be pre-sized according to:
64 : * \verbatim
65 : * aout.row() == ain.rows() + padRows
66 : * aout.cols() == ain.cols() + padCols
67 : * \endverbatim
68 : * where `padRows` and `padCols are the desired amount of padding. \p aout will be in
69 : * FFT order and suitable for input to the inverse transform (complex-to-complex).
70 : *
71 : * \ingroup ft
72 : */
73 : template <class eigenArrOutT, class eigenArrInT>
74 : void padC2CFFTOutput( eigenArrOutT &aout, /**< [out] The padded array. Must be pre-allocated. The pad region will
75 : not be set to zero. Size of the pad is determined by the
76 : difference in size between \p aout and \p ain*/
77 : eigenArrInT &ain /**< [in] The array to be paddded, in complex-to-complex
78 : FFTW storage order.*/ )
79 : {
80 : int ri = 0.5 * ain.rows() + 1;
81 : int ci = 0.5 * ain.cols() + 1;
82 :
83 : int rf = ain.rows() - ri;
84 : int cf = ain.cols() - ci;
85 :
86 : aout.topLeftCorner( ri, ci ) = ain.topLeftCorner( ri, ci );
87 : aout.bottomLeftCorner( rf, ci ) = ain.bottomLeftCorner( rf, ci );
88 : aout.topRightCorner( ri, cf ) = ain.topRightCorner( ri, cf );
89 : aout.bottomRightCorner( rf, cf ) = ain.bottomRightCorner( rf, cf );
90 : }
91 :
92 : /// Augment and pad the output of a real-to-complex (R2C) forward FFT for use with a complex-to-complex (C2C) inverse
93 : /// FFT
94 : /** The output array should be pre-sized according to:
95 : * \verbatim
96 : * aout.row() == 2*(ain.rows()-1) + padRows
97 : * aout.cols() == ain.cols() + padCols
98 : * \endverbatim
99 : * where padRows and padCols can be 0. \p aout will be in
100 : * FFT order and suitable for input to the inverse transform (complex-to-complex).
101 : *
102 : * \ingroup ft
103 : */
104 : template <class eigenArrOutT, class eigenArrInT>
105 0 : void augmentR2CFFTOutput( eigenArrOutT &aout, eigenArrInT &ain )
106 : {
107 :
108 0 : int padRows = aout.rows() - ain.rows();
109 0 : int padCols = aout.cols() - ain.cols();
110 :
111 : // aout.resize( 2 * ( ain.rows() - 1 ) + padRows, ain.cols() + padCols );
112 :
113 0 : int r1 = ain.rows();
114 0 : int c1 = 0.5 * ain.cols() + 1;
115 :
116 0 : int r2 = aout.rows() - ( r1 + padRows );
117 0 : int c2 = aout.cols() - ( c1 + padCols );
118 :
119 0 : aout.block( 0, 0, r1, c1 ) = ain.block( 0, 0, r1, c1 );
120 :
121 0 : aout.block( 0, c1 + padCols, r1, c2 ) = ain.block( 0, c1, r1, c2 );
122 :
123 : // Upper Right corner
124 0 : for( int cc = 1; cc < c1 - 1; ++cc )
125 : {
126 0 : for( int rr = 1; rr < r1 - 1; ++rr )
127 : {
128 0 : aout( aout.rows() - rr, aout.cols() - cc ) = conj( ain( rr, cc ) );
129 : }
130 : }
131 :
132 : // Lower Right corner
133 0 : for( int rr = 1; rr < ain.rows() - 1; ++rr )
134 : {
135 0 : aout( aout.rows() - rr, 0 ) = conj( ain( rr, 0 ) );
136 :
137 : // After ain.rows()+1, for column 1 and greater, the columns are reversed and conjugate
138 0 : aout.block( aout.rows() - rr, 1, 1, aout.cols() - padCols - c1 + 1 ) =
139 0 : conj( ain.block( rr, c1 - 1, 1, ain.cols() - c1 + 1 ).reverse() );
140 : }
141 0 : }
142 :
143 : } // namespace ft
144 : } // namespace math
145 : } // namespace mx
146 :
147 : #endif // ftUtils_hpp
|
