src/complexmat.cpp

   1 #include "complexmat.hpp"
   2
   3 ComplexMat_::T ComplexMat_::sqr_norm() const
   4 {
   5     assert(n_scales == 1);
   6
   7     int n_channels_per_scale = n_channels / n_scales;
   8     T sum_sqr_norm = 0;
   9     for (int i = 0; i < n_channels_per_scale; ++i) {
  10         for (auto lhs = p_data.hostMem() + i * rows * cols; lhs != p_data.hostMem() + (i + 1) * rows * cols; ++lhs)
  11             sum_sqr_norm += lhs->real() * lhs->real() + lhs->imag() * lhs->imag();
  12     }
  13     sum_sqr_norm = sum_sqr_norm / static_cast<T>(cols * rows);
  14     return sum_sqr_norm;
  15 }
  16
  17 void ComplexMat_::sqr_norm(DynMem_<ComplexMat_::T> &result) const
  18 {
  19     int n_channels_per_scale = n_channels / n_scales;
  20     int scale_offset = n_channels_per_scale * rows * cols;
  21     for (uint scale = 0; scale < n_scales; ++scale) {
  22         T sum_sqr_norm = 0;
  23         for (int i = 0; i < n_channels_per_scale; ++i)
  24             for (auto lhs = p_data.hostMem() + i * rows * cols + scale * scale_offset;
  25                  lhs != p_data.hostMem() + (i + 1) * rows * cols + scale * scale_offset; ++lhs)
  26                 sum_sqr_norm += lhs->real() * lhs->real() + lhs->imag() * lhs->imag();
  27         result.hostMem()[scale] = sum_sqr_norm / static_cast<T>(cols * rows);
  28     }
  29     return;
  30 }
  31
  32 ComplexMat_ ComplexMat_::sqr_mag() const
  33 {
  34     return mat_const_operator([](std::complex<T> &c) { c = c.real() * c.real() + c.imag() * c.imag(); });
  35 }
  36
  37 ComplexMat_ ComplexMat_::conj() const
  38 {
  39     return mat_const_operator([](std::complex<T> &c) { c = std::complex<T>(c.real(), -c.imag()); });
  40 }
  41
  42 ComplexMat_ ComplexMat_::sum_over_channels() const
  43 {
  44     assert(p_data.num_elem == n_channels * rows * cols);
  45
  46     uint n_channels_per_scale = n_channels / n_scales;
  47     uint scale_offset = n_channels_per_scale * rows * cols;
  48
  49     ComplexMat_ result(this->rows, this->cols, 1, n_scales);
  50     for (uint scale = 0; scale < n_scales; ++scale) {
  51         for (uint i = 0; i < rows * cols; ++i) {
  52             std::complex<T> acc = 0;
  53             for (uint ch = 0; ch < n_channels_per_scale; ++ch)
  54                 acc +=  p_data[scale * scale_offset + i + ch * rows * cols];
  55             result.p_data.hostMem()[scale * rows * cols + i] = acc;
  56         }
  57     }
  58     return result;
  59 }
  60
  61 ComplexMat_ ComplexMat_::operator/(const ComplexMat_ &rhs) const
  62 {
  63     return mat_mat_operator([](std::complex<T> &c_lhs, const std::complex<T> &c_rhs) { c_lhs /= c_rhs; }, rhs);
  64 }
  65
  66 ComplexMat_ ComplexMat_::operator+(const ComplexMat_ &rhs) const
  67 {
  68     return mat_mat_operator([](std::complex<T> &c_lhs, const std::complex<T> &c_rhs) { c_lhs += c_rhs; }, rhs);
  69 }
  70
  71 ComplexMat_ ComplexMat_::operator*(const ComplexMat_::T &rhs) const
  72 {
  73     return mat_const_operator([&rhs](std::complex<T> &c) { c *= rhs; });
  74 }
  75
  76 ComplexMat_ ComplexMat_::mul(const ComplexMat_ &rhs) const
  77 {
  78     return matn_mat1_operator([](std::complex<T> &c_lhs, const std::complex<T> &c_rhs) { c_lhs *= c_rhs; }, rhs);
  79 }
  80
  81 ComplexMat_ ComplexMat_::operator+(const ComplexMat_::T &rhs) const
  82 {
  83     return mat_const_operator([&rhs](std::complex<T> &c) { c += rhs; });
  84 }
  85
  86 ComplexMat_ ComplexMat_::operator*(const ComplexMat_ &rhs) const
  87 {
  88     return mat_mat_operator([](std::complex<T> &c_lhs, const std::complex<T> &c_rhs) { c_lhs *= c_rhs; }, rhs);
  89 }
  90
  91 ComplexMat_ ComplexMat_::mat_mat_operator(void (*op)(std::complex<ComplexMat_::T> &, const std::complex<ComplexMat_::T> &), const ComplexMat_ &mat_rhs) const
  92 {
  93     assert(mat_rhs.n_channels == n_channels/n_scales && mat_rhs.cols == cols && mat_rhs.rows == rows);
  94
  95     ComplexMat_ result = *this;
  96     for (uint s = 0; s < n_scales; ++s) {
  97         auto lhs = result.p_data.hostMem() + (s * n_channels/n_scales * rows * cols);
  98         auto rhs = mat_rhs.p_data.hostMem();
  99         for (uint i = 0; i < n_channels/n_scales * rows * cols; ++i)
 100             op(*(lhs + i), *(rhs + i));
 101     }
 102
 103     return result;
 104 }
 105
 106 ComplexMat_ ComplexMat_::matn_mat1_operator(void (*op)(std::complex<ComplexMat_::T> &, const std::complex<ComplexMat_::T> &), const ComplexMat_ &mat_rhs) const
 107 {
 108     assert(mat_rhs.n_channels == 1 && mat_rhs.cols == cols && mat_rhs.rows == rows);
 109
 110     ComplexMat_ result = *this;
 111     for (uint i = 0; i < n_channels; ++i) {
 112         auto lhs = result.p_data.hostMem() + i * rows * cols;
 113         auto rhs = mat_rhs.p_data.hostMem();
 114         for (; lhs != result.p_data.hostMem() + (i + 1) * rows * cols; ++lhs, ++rhs)
 115             op(*lhs, *rhs);
 116     }
 117
 118     return result;
 119 }
 120
 121 ComplexMat_ ComplexMat_::matn_mat2_operator(void (*op)(std::complex<ComplexMat_::T> &, const std::complex<ComplexMat_::T> &), const ComplexMat_ &mat_rhs) const
 122 {
 123     assert(mat_rhs.n_channels == n_channels / n_scales && mat_rhs.cols == cols && mat_rhs.rows == rows);
 124
 125     int n_channels_per_scale = n_channels / n_scales;
 126     int scale_offset = n_channels_per_scale * rows * cols;
 127     ComplexMat_ result = *this;
 128     for (uint i = 0; i < n_scales; ++i) {
 129         for (int j = 0; j < n_channels_per_scale; ++j) {
 130             auto lhs = result.p_data.hostMem() + (j * rows * cols) + (i * scale_offset);
 131             auto rhs = mat_rhs.p_data.hostMem() + (j * rows * cols);
 132             for (; lhs != result.p_data.hostMem() + ((j + 1) * rows * cols) + (i * scale_offset); ++lhs, ++rhs)
 133                 op(*lhs, *rhs);
 134         }
 135     }
 136
 137     return result;
 138 }
 139
 140 ComplexMat_ ComplexMat_::mat_const_operator(const std::function<void (std::complex<ComplexMat_::T> &)> &op) const
 141 {
 142     ComplexMat_ result = *this;
 143     for (uint i = 0; i < n_channels; ++i) {
 144         for (auto lhs = result.p_data.hostMem() + i * rows * cols;
 145              lhs != result.p_data.hostMem() + (i + 1) * rows * cols; ++lhs)
 146             op(*lhs);
 147     }
 148     return result;
 149 }