|
| 1 | +#include "ElNullAnalyses.h" |
| 2 | + |
| 3 | +#include <algorithm> |
| 4 | +#include <cmath> |
| 5 | +#include <numeric> |
| 6 | + |
| 7 | +#include <isce3/core/Interp1d.h> |
| 8 | +#include <isce3/core/Kernels.h> |
| 9 | +#include <isce3/except/Error.h> |
| 10 | +#include <isce3/focus/RangeComp.h> |
| 11 | + |
| 12 | +namespace isce3 { namespace antenna { |
| 13 | + |
| 14 | +Eigen::ArrayXcd linearInterpComplex1d( |
| 15 | + const Eigen::Ref<const Eigen::ArrayXd>& x0, const Linspace_t& x, |
| 16 | + const Eigen::Ref<const Eigen::ArrayXcd>& y) |
| 17 | +{ |
| 18 | + // check the input arguments |
| 19 | + if (!(x.spacing() > 0.0)) |
| 20 | + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), |
| 21 | + "The argument 'x' must be monotonically increasing!"); |
| 22 | + if (x.size() != y.size()) |
| 23 | + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), |
| 24 | + "The arguments 'x' and 'y' must have the same size!"); |
| 25 | + if ((x0.minCoeff() < x.first()) || (x0.maxCoeff() > x.last())) |
| 26 | + throw isce3::except::InvalidArgument( |
| 27 | + ISCE_SRCINFO(), "The 'x0' values are out of range of 'x'!"); |
| 28 | + |
| 29 | + // form the linear kernel |
| 30 | + auto linear_kernel = isce3::core::LinearKernel<double>(); |
| 31 | + // allocate interpolated output vector |
| 32 | + Eigen::ArrayXcd y0(x0.size()); |
| 33 | +#pragma omp parallel for |
| 34 | + for (Eigen::Index idx = 0; idx < x0.size(); ++idx) { |
| 35 | + auto x_int = (x0(idx) - x.first()) / x.spacing(); |
| 36 | + y0(idx) = isce3::core::interp1d( |
| 37 | + linear_kernel, y.data(), y.size(), 1, x_int, false); |
| 38 | + } |
| 39 | + return y0; |
| 40 | +} |
| 41 | + |
| 42 | +tuple_ant genAntennaPairCoefs( |
| 43 | + const Eigen::Ref<const Eigen::ArrayXcd>& el_cut_left, |
| 44 | + const Eigen::Ref<const Eigen::ArrayXcd>& el_cut_right, |
| 45 | + double el_ang_start, double el_ang_step, |
| 46 | + std::optional<double> el_res_max) |
| 47 | +{ |
| 48 | + // check input arguments |
| 49 | + if (!(el_ang_step > 0.0)) |
| 50 | + throw isce3::except::InvalidArgument( |
| 51 | + ISCE_SRCINFO(), "EL angle step must be a positive value!"); |
| 52 | + |
| 53 | + if (el_res_max) { |
| 54 | + if (!(*el_res_max > 0.0)) |
| 55 | + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), |
| 56 | + "Max EL angle resolution must be a positive value!"); |
| 57 | + } else |
| 58 | + *el_res_max = el_ang_step; |
| 59 | + |
| 60 | + const auto cut_size {el_cut_left.size()}; |
| 61 | + if (el_cut_right.size() != cut_size) |
| 62 | + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), |
| 63 | + "Size mismatch between two EL-cut vectors left and right!"); |
| 64 | + |
| 65 | + // pick EL cut patterns within peak-to-peak power gains of two adjacent |
| 66 | + // beams, resample/interpolate the picked values into finer el-resolution if |
| 67 | + // necessary based on min (_el_res_max, el_ang_step) and store their |
| 68 | + // conjugate versions as final weighting coefs as a function EL angle vector |
| 69 | + // "el_ang_vec" in (rad). |
| 70 | + Eigen::Index idx_peak_left; |
| 71 | + Eigen::Index idx_peak_right; |
| 72 | + el_cut_left.abs().maxCoeff(&idx_peak_left); |
| 73 | + el_cut_right.abs().maxCoeff(&idx_peak_right); |
| 74 | + if (idx_peak_left >= idx_peak_right) |
| 75 | + throw isce3::except::InvalidArgument(ISCE_SRCINFO(), |
| 76 | + "Wrong order or incorrect EL cuts of left and right beams!"); |
| 77 | + |
| 78 | + // take out one sample on each side (that is to exclude the peaks) |
| 79 | + idx_peak_left += 1; |
| 80 | + idx_peak_right -= 1; |
| 81 | + // check to have at least 3 points! |
| 82 | + auto num_idx_p2p = idx_peak_right - idx_peak_left + 1; |
| 83 | + if (num_idx_p2p < 3) |
| 84 | + throw isce3::except::RuntimeError(ISCE_SRCINFO(), |
| 85 | + "There is not enough samples/separation " |
| 86 | + "between peaks of EL cuts!"); |
| 87 | + |
| 88 | + // form output uniform EL angle vector within [idx_peak_left, |
| 89 | + // idx_peak_right] |
| 90 | + auto el_space = std::min(*el_res_max, el_ang_step); |
| 91 | + auto el_start = idx_peak_left * el_ang_step + el_ang_start; |
| 92 | + auto el_stop = idx_peak_right * el_ang_step + el_ang_start; |
| 93 | + auto num_idx = |
| 94 | + static_cast<Eigen::Index>((el_stop - el_start) / el_space + 1); |
| 95 | + |
| 96 | + Eigen::ArrayXd el_vec = |
| 97 | + Eigen::ArrayXd::LinSpaced(num_idx, el_start, el_stop); |
| 98 | + |
| 99 | + // interpolate the complex coeff within [idx_peak_left, idx_peak_right] if |
| 100 | + // necessary |
| 101 | + if (*el_res_max < el_ang_step) |
| 102 | + // perform linear interpolation of complex coeffs |
| 103 | + { |
| 104 | + auto el_ang_linspace = Linspace_t(el_ang_start, el_ang_step, cut_size); |
| 105 | + Eigen::ArrayXcd coef_left = Eigen::conj( |
| 106 | + linearInterpComplex1d(el_vec, el_ang_linspace, el_cut_left)); |
| 107 | + Eigen::ArrayXcd coef_right = Eigen::conj( |
| 108 | + linearInterpComplex1d(el_vec, el_ang_linspace, el_cut_right)); |
| 109 | + return std::make_tuple(coef_left, coef_right, el_vec); |
| 110 | + } |
| 111 | + // no interpolation |
| 112 | + return std::make_tuple( |
| 113 | + Eigen::conj(el_cut_left.segment(idx_peak_left, num_idx)), |
| 114 | + Eigen::conj(el_cut_right.segment(idx_peak_left, num_idx)), el_vec); |
| 115 | +} |
| 116 | + |
| 117 | +std::tuple<double, Eigen::Index, double> locateAntennaNull( |
| 118 | + const Eigen::Ref<const Eigen::ArrayXcd>& coef_left, |
| 119 | + const Eigen::Ref<const Eigen::ArrayXcd>& coef_right, |
| 120 | + const Eigen::Ref<const Eigen::ArrayXd>& el_ang_vec) |
| 121 | +{ |
| 122 | + // antenna null power pattern = |left**2 - right**2|/(left**2 + right**2) |
| 123 | + // Its min is the null location. |
| 124 | + auto pow_ant_left = coef_left.abs2(); |
| 125 | + auto pow_ant_right = coef_right.abs2(); |
| 126 | + Eigen::ArrayXd pow_null_ant = (pow_ant_left - pow_ant_right).abs() / |
| 127 | + (pow_ant_left + pow_ant_right); |
| 128 | + // locate the null from power pattern |
| 129 | + Eigen::Index idx_null_ant; |
| 130 | + double min_val = pow_null_ant.minCoeff(&idx_null_ant); |
| 131 | + double max_val = pow_null_ant.maxCoeff(); |
| 132 | + if (!(max_val > min_val)) |
| 133 | + throw isce3::except::RuntimeError(ISCE_SRCINFO(), |
| 134 | + "The antenna null pattern does not have a dip!"); |
| 135 | + // get the peak-normalized null value in linear scale |
| 136 | + double val_null {min_val / max_val}; |
| 137 | + // get EL angle of the null in (rad) |
| 138 | + double el_null_ant = el_ang_vec(idx_null_ant); |
| 139 | + return {el_null_ant, idx_null_ant, val_null}; |
| 140 | +} |
| 141 | + |
| 142 | +tuple_echo formEchoNull(const std::vector<std::complex<float>>& chirp_ref, |
| 143 | + const Eigen::Ref<const RowMatrixXcf>& echo_left, |
| 144 | + const Eigen::Ref<const RowMatrixXcf>& echo_right, double sr_start, |
| 145 | + double sr_spacing, const Eigen::Ref<const Eigen::ArrayXcd>& coef_left, |
| 146 | + const Eigen::Ref<const Eigen::ArrayXcd>& coef_right, |
| 147 | + const Eigen::Ref<const Eigen::ArrayXd>& sr_coef) |
| 148 | +{ |
| 149 | + // form rangecomp obj used for all range comp in forming echo null pattern |
| 150 | + auto rgc_obj = isce3::focus::RangeComp(chirp_ref, echo_left.cols(), 1, |
| 151 | + isce3::focus::RangeComp::Mode::Valid); |
| 152 | + // final number of range bins for the echo after range comp |
| 153 | + const auto num_rgb_echo = rgc_obj.outputSize(); |
| 154 | + // form uniform slant range (m) vector for only valid part of final |
| 155 | + // rangecomp echo |
| 156 | + const auto sr_stop = sr_start + (num_rgb_echo - 1) * sr_spacing; |
| 157 | + Eigen::ArrayXd sr_echo_valid = |
| 158 | + Eigen::ArrayXd::LinSpaced(num_rgb_echo, sr_start, sr_stop); |
| 159 | + |
| 160 | + // get limited [first, last] indices by intersecting uniform slant ranges |
| 161 | + // from only valid-mode rangecomp echo with that of non-uniform one from |
| 162 | + // antenna weighting coefs |
| 163 | + auto [idx_coef_first, idx_coef_last, idx_echo_first, idx_echo_last] = |
| 164 | + detail::intersect_nearest(sr_coef, sr_echo_valid); |
| 165 | + const auto size_coef = idx_coef_last - idx_coef_first + 1; |
| 166 | + const auto num_rgb_null = idx_echo_last - idx_echo_first + 1; |
| 167 | + // now get array of relative indices for mapping uniform sr_echo[first, |
| 168 | + // last] to non-uniform sr_coeff[first, last] |
| 169 | + auto idx_coef_vec = |
| 170 | + detail::locate_nearest(sr_coef.segment(idx_coef_first, size_coef), |
| 171 | + sr_echo_valid.segment(idx_echo_first, num_rgb_null)); |
| 172 | + // add the first index to all values for absolute mapping from echo[first, |
| 173 | + // last] to antenna coeff. |
| 174 | + idx_coef_vec += idx_coef_first; |
| 175 | + |
| 176 | + // get limited coefs left/right within non-uniform indices "idx_coef_vec" |
| 177 | + // related to uniform range bins of echo. |
| 178 | + Eigen::ArrayXcd coef_left_limit(num_rgb_null); |
| 179 | + Eigen::ArrayXcd coef_right_limit(num_rgb_null); |
| 180 | + for (Eigen::Index idx = 0; idx < idx_coef_vec.size(); ++idx) { |
| 181 | + coef_left_limit(idx) = coef_left(idx_coef_vec(idx)); |
| 182 | + coef_right_limit(idx) = coef_right(idx_coef_vec(idx)); |
| 183 | + } |
| 184 | + |
| 185 | + // initialize the peak-normalized averaged echo null power |
| 186 | + Eigen::ArrayXd pow_null_avg = Eigen::ArrayXd::Zero(num_rgb_null); |
| 187 | + // allocate rangeline vector for range compression of left/right echoes |
| 188 | + Eigen::ArrayXcf rgc_left(num_rgb_echo); |
| 189 | + Eigen::ArrayXcf rgc_right(num_rgb_echo); |
| 190 | + // allocate double precision lines for left and right weighted rangecomp |
| 191 | + // echo within null formation part only |
| 192 | + Eigen::ArrayXcd line_left(num_rgb_null); |
| 193 | + Eigen::ArrayXcd line_right(num_rgb_null); |
| 194 | + // loop over range lines /pulses |
| 195 | + for (Eigen::Index pulse = 0; pulse < echo_left.rows(); ++pulse) { |
| 196 | + // range compression of echoes left/right |
| 197 | + rgc_obj.rangecompress(rgc_left.data(), echo_left.row(pulse).data()); |
| 198 | + rgc_obj.rangecompress(rgc_right.data(), echo_right.row(pulse).data()); |
| 199 | + |
| 200 | + // get the power and add up its double precision version |
| 201 | + line_left = rgc_left.segment(idx_echo_first, num_rgb_null) |
| 202 | + .cast<std::complex<double>>() * |
| 203 | + coef_left_limit; |
| 204 | + line_right = rgc_right.segment(idx_echo_first, num_rgb_null) |
| 205 | + .cast<std::complex<double>>() * |
| 206 | + coef_right_limit; |
| 207 | + |
| 208 | + // form the null power to be averaged over range lines |
| 209 | + pow_null_avg += |
| 210 | + (line_left - line_right).abs() / (line_left + line_right).abs(); |
| 211 | + } |
| 212 | + auto max_pow_null = pow_null_avg.maxCoeff(); |
| 213 | + if (!(max_pow_null > pow_null_avg.minCoeff())) |
| 214 | + throw isce3::except::RuntimeError(ISCE_SRCINFO(), |
| 215 | + "The echo null power pattern does not have a dip!"); |
| 216 | + // peak-normalized power pattern |
| 217 | + pow_null_avg /= max_pow_null; |
| 218 | + |
| 219 | + // power (linear), slant range (m) , index for EL angles (-) |
| 220 | + return std::make_tuple(pow_null_avg, |
| 221 | + sr_coef.segment(idx_coef_vec(0), num_rgb_null), idx_coef_vec); |
| 222 | +} |
| 223 | + |
| 224 | +}} // namespace isce3::antenna |
0 commit comments