针对新ArrayFire版本的问题自适应MatchFilter算法

Question

我正在尝试将匹配的过滤器算法（给定为here）用于arrayfire版本3.6.4。这是我最终得到的结果：

#include <arrayfire.h>

using namespace af;

struct SAR_data { //! SAR data structure format
    double deltaF;  //! step size of frequency data (Hz)
    double maxWr;   //! maximum scene size in range direction (m)
    double maxWx;   //! maximum scene size in cross-range direction (m)

    array minF;   //! (Np x 1): Vector containing the start frequency of 
                  //! each pulse (Hz)
    array x_mat;  //! (Sx x Sy): the x-positions of each pixel (m)
    array y_mat;  //! (Sx x Sy): the y-positions of each pixel (m)
    array z_mat;  //! (Sx x Sy): the z-positions of each pixel (m)

    array AntX;   //! (Np x 1): the x-position of the sensor at each pulse (m)
    array AntY;   //! (Np x 1): the y-position of the sensor at each pulse (m)
    array AntZ;   //! (Np x 1): the z-position of the sensor at each pulse (m)

    array R0;     //! (Np x 1): the range to scene center (m)
    array phdata; //! (K x Np): phase history data (frequency domain),
                  //! fast time in rows, slow time in columns
    array im_final; //! (Sx x Sy): the complex image value at each pixel
};      

// matched filter realization using ArrayFire
void matched_filter_ArrayFire(SAR_data& data) {
   double c = 299792458.0; // speed of light

   // Determine the size of the phase history data
   int K  = data.phdata.dims(0);  // # of frequency bins per pulse
   int Np = data.phdata.dims(1);  // # of pulses

   // Initialize the image with all zero values (complex)
   data.im_final = constant<cdouble>(0, data.x_mat.dims(), c64);
   array im_slices = constant<cdouble>(0, K, data.x_mat.dims(0),
                              data.x_mat.dims(1), c64);
   array fspan = array(seq(0.0, K-1)) * data.deltaF;

   cdouble unit = {0, 1};
   for (int ii = 0; ii < Np; ii++) {
      // compute differential range for each image pixel (m)
      array dx = data.AntX(ii) - data.x_mat;
      array dy = data.AntY(ii) - data.y_mat;
      array dz = data.AntZ(ii) - data.z_mat;
      array dR = sqrt(dx*dx + dy*dy + dz*dz) - data.R0(ii);

      // calculate the frequency of each sample in the pulse (Hz)
      array freq = data.minF(ii) + fspan;
      array tt = data.phdata(span,ii);

      // perform the Matched Filter operation
      gfor (array jj, K) {
         im_slices(jj,span,span) = 
            tt(jj)*exp(unit*((double)(4.0*Pi/c)*freq(jj)*dR)); 
      }
      tt = sum(im_slices,0);
      data.im_final = data.im_final + moddims(tt, data.x_mat.dims());
   }
}

代码编译良好，但是当我尝试运行它时，它会引发异常：

What() is:ArrayFire Exception (Invalid input size:203):
In function class af::dim4 __cdecl getOutDims(const class af::dim4 &,const class af::dim4 &,bool)
In file src\backend\common\ArrayInfo.cpp:130
Invalid dimension for argument 1
Expected: ldims == rdims

In function class af::array __cdecl af::operator -(const class af::array &,const class af::array &)
In file src\api\cpp\array.cpp:838

我发现，它抱怨以下几行：

      array dx = data.AntX(ii) - data.x_mat;
      array dy = data.AntY(ii) - data.y_mat;
      array dz = data.AntZ(ii) - data.z_mat;

因为data.AntX（ii）（即1x1）和data.x_mat（即Sx x Sy）的尺寸不匹配。不过，在旧版本的ArrayFire中，它过去可以正常工作。显然，我需要某种方法告诉ArrayFire为data.x_mat的整个矩阵维“扩展” data.AntX（ii）的值。但是怎么做？

Answer 1

使用@pradeep提示的算法的工作版本：

// matched filter realization using ArrayFire
void matched_filter_ArrayFire(SAR_data& data) {
   double c = 299792458.0; // speed of light

   // Determine the size of the phase history data
   int K  = data.phdata.dims(0);  // # of frequency bins per pulse
   int Np = data.phdata.dims(1);  // # of pulses

   // Initialize the image with all zero values (complex)
   data.im_final = constant<cdouble>({0,0}, data.x_mat.dims(), c64);
   array im_slices = constant<cdouble>({0,0}, data.x_mat.dims(0),
                              data.x_mat.dims(1), K, c64);
   array fspan = array(seq(0.0, K-1)) * data.deltaF;

   cdouble unit = {0, 1};
   for (int ii = 0; ii < Np; ii++) {
      // compute differential range for each image pixel (m)
      array dx = tile(data.AntX(ii), data.x_mat.dims()) - data.x_mat;
      array dy = tile(data.AntY(ii), data.y_mat.dims()) - data.y_mat;
      array dz = tile(data.AntZ(ii), data.x_mat.dims()) - data.z_mat;
      array dR = sqrt(dx*dx + dy*dy + dz*dz) - tile(data.R0(ii), dx.dims());

      // calculate the frequency of each sample in the pulse (Hz)
      array freq = tile(data.minF(ii), fspan.dims()) + fspan;

      // perform the Matched Filter operation
      gfor(seq jj, K) {
          im_slices(span,span,jj) = tile(data.phdata(jj,ii), dR.dims())*
                exp(unit*(double)(4.0*Pi/c)*tile(freq(jj), dR.dims())*dR); 
      }
      // sum over the last dimension
      data.im_final += moddims(sum(im_slices,2),data.x_mat.dims());
   }
}