Line data Source code
1 : /** \file speckleAmpPSD.hpp
2 : * \author Jared R. Males (jaredmales@gmail.com)
3 : * \brief Calculates the PSD of speckle intensity given a modified Fourier mode amplitude PSD.
4 : * \ingroup mxAO_files
5 : *
6 : */
7 :
8 : #ifndef speckleAmpPSD_hpp
9 : #define speckleAmpPSD_hpp
10 :
11 : #include <vector>
12 :
13 : #include <Eigen/Dense>
14 :
15 : #include "../../math/constants.hpp"
16 : #include "../../math/randomT.hpp"
17 : #include "../../math/vectorUtils.hpp"
18 :
19 : #include "../../sigproc/signalWindows.hpp"
20 : #include "../../sigproc/psdUtils.hpp"
21 : #include "../../sigproc/psdFilter.hpp"
22 : #include "../../sigproc/psdVarMean.hpp"
23 : #include "../../sigproc/averagePeriodogram.hpp"
24 :
25 : #include "jincFuncs.hpp"
26 :
27 : namespace mx
28 : {
29 : namespace AO
30 : {
31 : namespace analysis
32 : {
33 :
34 : /// Calculate the PSD of the speckle intensity given the PSD of Fourier mode amplitude
35 : /**
36 : * Does so by generating a time-series from the PSD using Fourier-domain convolution, and calculating
37 : * the average periodogram.
38 : *
39 : * A Hann window is used.
40 : *
41 : * Statistics are not currently normalized. You need to normalize to match expected contrast if desired.
42 : *
43 : * Will also calculate the binned-variances in the generated time-series, if the arguments vars and bins are not null
44 : * pointers.
45 : *
46 : * \todo Figure out what's going on with psdFilter normalization!
47 : * \todo probably need to not do overlapped averaging in periodogram. Could drop mean-sub then.
48 : */
49 : template <typename realT>
50 0 : int speckleAmpPSD(
51 : std::vector<realT> &spFreq, ///< [out] The frequency grid of the output PSD
52 : std::vector<realT>
53 : &spPSD, ///< [out] The speckle amplitude PSD corresponding to the freq coordinates. Will be resized.
54 : const std::vector<realT> &freq, ///< [in] The Frequency grid of the input PSD
55 : const std::vector<realT> &fmPSD, ///< [in] The Fourier mode PSD.
56 : const std::vector<std::complex<realT>>
57 : &fmXferFxn, ///< [in] The complex error transfer function, as a function of freq
58 : const std::vector<realT> &nPSD, ///< [in] The open-loop noise PSD
59 : const std::vector<std::complex<realT>> &nXferFxn, ///< [in] The noise transfer function, as a function of freq
60 : int N, ///< [in] The number of trials to use in calculating the amplitude PSD.
61 : std::vector<realT> *vars = nullptr, ///< [out] [optional] The binned variances of the time series generated.
62 : std::vector<realT> *bins = nullptr, ///< [in] [optional] The bin sizes to use in calculating vars.
63 : bool noPSD = false ///< [in] [optional] if true then the PSD is not actually calculated, only the binned variances
64 : )
65 : {
66 :
67 0 : std::vector<realT> psd2, npsd2;
68 0 : std::vector<std::complex<realT>> xfer2, nxfer2;
69 :
70 0 : bool hasZero = false;
71 0 : if( freq[0] == 0 )
72 0 : hasZero = true;
73 :
74 : // First augment to two-sided DFT form
75 0 : mx::sigproc::augment1SidedPSD( psd2, fmPSD, !hasZero, 0.5 );
76 0 : mx::sigproc::augment1SidedPSD( npsd2, nPSD, !hasZero, 0.5 );
77 :
78 : // And augment the xfer fxns to two sided form
79 0 : sigproc::augment1SidedPSD( xfer2, fmXferFxn, !hasZero, 1.0 );
80 0 : sigproc::augment1SidedPSD( nxfer2, nXferFxn, !hasZero, 1.0 );
81 :
82 : // The time sampling
83 0 : realT dt = 1. / ( 2 * freq.back() );
84 :
85 : // Indices for getting the middle half
86 0 : int Nwd = 1.0 * psd2.size();
87 0 : int NwdStart = 0.5 * psd2.size() - 0.5 * Nwd;
88 :
89 0 : int Nsamp = 1.0 * psd2.size();
90 0 : int NsampStart = 0.5 * Nwd - 0.5 * Nsamp;
91 :
92 0 : sigproc::averagePeriodogram<realT> globalAvgPgram( Nsamp * 0.1, dt );
93 0 : spPSD.resize( globalAvgPgram.size() );
94 0 : for( size_t n = 0; n < spPSD.size(); ++n )
95 0 : spPSD[n] = 0;
96 :
97 0 : std::vector<std::vector<realT>> means;
98 0 : if( bins != nullptr && vars != nullptr )
99 0 : means.resize( bins->size() );
100 :
101 : // Calculate the Fourier Mode variance for normalization
102 0 : realT fmVar = sigproc::psdVar( freq, fmPSD );
103 : // and the noise variance
104 0 : realT nVar = sigproc::psdVar( freq, nPSD );
105 :
106 0 : #pragma omp parallel
107 : {
108 : // Filters for imposing the PSDs
109 : mx::sigproc::psdFilter<realT, 1> filt;
110 : filt.psd( psd2, freq[1] - freq[0] );
111 :
112 : mx::sigproc::psdFilter<realT, 1> nfilt;
113 : nfilt.psd( npsd2, freq[1] - freq[0] );
114 :
115 : // FFTs for going to Fourier domain and back to time domain.
116 : math::ft::fftT<realT, std::complex<realT>, 1, 0> fft( psd2.size() );
117 : math::ft::fftT<std::complex<realT>, realT, 1, 0> fftB( psd2.size(), math::ft::dir::backward );
118 :
119 : // Fourier transform working memmory
120 : std::vector<std::complex<realT>> tform1( psd2.size() );
121 : std::vector<std::complex<realT>> tform2( psd2.size() );
122 :
123 : std::vector<std::complex<realT>> Ntform1( psd2.size() );
124 : std::vector<std::complex<realT>> Ntform2( psd2.size() );
125 :
126 : // Normally distributed random numbers
127 : math::normDistT<realT> normVar;
128 : normVar.seed();
129 :
130 : // The two modal amplitude series
131 : std::vector<realT> fm_n( psd2.size() );
132 : std::vector<realT> fm_nm( psd2.size() );
133 :
134 : // The two noise amplitude series
135 : std::vector<realT> N_n( psd2.size() );
136 : std::vector<realT> N_nm( psd2.size() );
137 :
138 : // The speckle amplitude
139 : std::vector<realT> vnl( Nwd );
140 : std::vector<realT> vn( Nsamp );
141 :
142 : // Periodogram averager
143 : sigproc::averagePeriodogram<realT> avgPgram( Nsamp * 0.1, dt ); //, 0, 1);
144 : avgPgram.win( sigproc::window::hann );
145 :
146 : // The temporary periodogram
147 : std::vector<realT> tpgram( avgPgram.size() ); // spPSD.size());
148 :
149 : #pragma omp for
150 : for( int k = 0; k < N; ++k )
151 : {
152 : // Generate the time-series
153 : for( int i = 0; i < psd2.size(); ++i )
154 : {
155 : fm_n[i] = normVar;
156 : N_n[i] = normVar;
157 : N_nm[i] = normVar;
158 : }
159 :
160 : // Filter and normalize the fourier mode time series
161 : filt( fm_n );
162 : math::vectorMeanSub( fm_n );
163 : realT actvar = math::vectorVariance( fm_n );
164 : realT norm = sqrt( fmVar / actvar );
165 : for( size_t q = 0; q < fm_n.size(); ++q )
166 : fm_n[q] *= norm;
167 :
168 : // And move it to the Fourier domain
169 : fft( tform1.data(), fm_n.data() );
170 :
171 : // Filter and normalize the measurement mode time series
172 : nfilt.filter( N_n );
173 : nfilt.filter( N_nm );
174 :
175 : realT Nactvar = 0.5 * ( math::vectorVariance( N_n ) + math::vectorVariance( N_nm ) );
176 : norm = sqrt( nVar / Nactvar );
177 : for( size_t q = 0; q < fm_n.size(); ++q )
178 : N_n[q] *= norm;
179 : for( size_t q = 0; q < fm_n.size(); ++q )
180 : N_nm[q] *= norm;
181 :
182 : // And move them to the Fourier domain
183 : fft( Ntform1.data(), N_n.data() );
184 : fft( Ntform2.data(), N_nm.data() );
185 :
186 : std::complex<realT> scale = std::complex<realT>( ( tform1.size() ), 0 );
187 :
188 : // Apply the modal phase shift, and apply the measurement noise.
189 : for( size_t m = 0; m < tform1.size(); ++m )
190 : {
191 : // Apply the phase shift to form the 2nd time series
192 : tform2[m] = tform1[m] * exp( std::complex<realT>( 0, math::half_pi<realT>() ) );
193 :
194 : // Apply the augmented ETF to two time-series
195 : tform1[m] *= xfer2[m] / scale;
196 : tform2[m] *= xfer2[m] / scale;
197 :
198 : // Ntform2[m] = Ntform1[m]*exp( std::complex<realT>(0, math::half_pi<realT>() ));
199 :
200 : Ntform1[m] *= nxfer2[m] / scale;
201 : Ntform2[m] *= nxfer2[m] / scale;
202 : }
203 :
204 : //<<<<<<<<****** Transform back to the time domain.
205 : fftB( fm_n.data(), tform1.data() );
206 : fftB( fm_nm.data(), tform2.data() );
207 : fftB( N_n.data(), Ntform1.data() );
208 : fftB( N_nm.data(), Ntform2.data() );
209 :
210 : // Calculate the speckle amplitude and mean-subtract
211 : realT mn = 0;
212 : for( int i = 0; i < Nwd; ++i )
213 : {
214 : realT h1 = fm_n[i + NwdStart] + N_n[i + NwdStart];
215 : realT h2 = fm_nm[i + NwdStart] + N_nm[i + NwdStart];
216 :
217 : vnl[i] = ( pow( h1, 2 ) + pow( h2, 2 ) );
218 : mn += vnl[i];
219 : }
220 : mn /= vnl.size();
221 :
222 : // Extract middle sample and mean subtract
223 : for( int i = 0; i < Nsamp; ++i )
224 : vn[i] = vnl[i + NsampStart] - mn;
225 :
226 : // Calculate PSD of the speckle amplitude
227 : if( !noPSD )
228 : avgPgram( tpgram, vn );
229 :
230 : // Accumulate
231 : #pragma omp critical
232 : {
233 : if( bins != nullptr && vars != nullptr )
234 : math::vectorBinMeans( means, *bins, vn );
235 : for( size_t i = 0; i < spPSD.size(); ++i )
236 : spPSD[i] += tpgram[i];
237 : }
238 : }
239 : } // pragma omp parallel
240 :
241 0 : if( !noPSD )
242 : {
243 : //-- Calculate frequency grid (which is maybe not identical to input freq)
244 0 : realT df = ( 1.0 ) / ( 2.0 * spPSD.size() * dt );
245 0 : spFreq.resize( spPSD.size() );
246 0 : for( size_t i = 0; i < spFreq.size(); ++i )
247 : {
248 0 : spFreq[i] = globalAvgPgram[i];
249 0 : spPSD[i] /= N;
250 : }
251 : // realT spVar = sigproc::psdVar1sided(df, spPSD.data(), spPSD.size());
252 : // for(size_t i=0; i< spPSD.size(); ++i) spPSD[i] *= (fmVar*fmVar/spVar);
253 : }
254 :
255 : // Calculate binned variances.
256 0 : if( bins != nullptr && vars != nullptr )
257 : {
258 0 : vars->resize( bins->size() );
259 0 : for( size_t i = 0; i < bins->size(); ++i )
260 : {
261 0 : ( *vars )[i] = math::vectorVariance( means[i] );
262 : }
263 : }
264 :
265 0 : return 0;
266 0 : }
267 :
268 : /// Calculate the variance of the mean vs bin size in a speckle time-series
269 : /**
270 : * Does so by generating a time-series from the PSD using Fourier-domain convolution, then
271 : * calculates the binned-variances in the generated time-series.
272 : *
273 : * This is just a wrapper for speckleAmpPSD, with no actual PSD calculation.
274 : *
275 : */
276 : template <typename realT>
277 : int speckleAmpVarMean(
278 : std::vector<realT> &vars, ///< [out] The binned variances of the time series generated.
279 : std::vector<realT> &bins, ///< [in] The bin sizes to use to calculate the variances.
280 : std::vector<realT> &freq, ///< [in] The Frequency grid of the input PSD
281 : std::vector<realT> &fmPSD, ///< [in] The open-loop Fourier mode PSD.
282 : std::vector<std::complex<realT>> &fmXferFxn, ///< [in] The complex error transfer function, as a function of freq
283 : std::vector<realT> &nPSD, ///< [in] The open-loop noise PSD
284 : std::vector<std::complex<realT>> &nXferFxn, ///< [in] The noise transfer function, as a function of freq
285 : int N ///< [in] The number of trials to use in calculating the amplitude time-series.
286 : )
287 : {
288 : std::vector<realT> spFreq, spPSD;
289 :
290 : return speckleAmpPSD( spFreq, spPSD, freq, fmPSD, fmXferFxn, nPSD, nXferFxn, N, &vars, &bins, true );
291 : }
292 :
293 : } // namespace analysis
294 : } // namespace AO
295 : } // namespace mx
296 :
297 : #endif // speckleAmpPSD_hpp
|
