COB-53: Only use delay gradient within a block for Doppler correction, the other terms are constant

RTCP/Cobalt/GPUProc/share/gpu/kernels/FIR_Filter.cu

@@ -119,23 +119,19 @@ inline __device__ float sincos_f2f_select(float phi, int pol)
 typedef  const double(*DelaysType)[1][NR_STABS][NR_POLARIZATIONS]; // 2 Polarizations; in seconds
 typedef  const double(*Phase0sType)[NR_STABS][NR_POLARIZATIONS]; // 2 Polarizations; in radians
 
-inline __device__ float correction_factor(unsigned time, double phiAtBegin, double phiAfterEnd, int ri, int pol, unsigned channel)
+inline __device__ float correction_factor(unsigned time, double phiGradient, int pol, unsigned channel)
 {
-    // Offset of this sample between begin and end.
+    // Offset of this sample between begin and end. = t/T fraction
     const double blockOffset = double(time * NR_CHANNELS + channel) / NR_SAMPLES_PER_SUBBAND;
-    //const double blockOffset = double(channel * NR_SAMPLES_PER_CHANNEL + time) / (double)NR_SAMPLES_PER_SUBBAND;
     // Interpolate the required phase rotation for this sample.
     //
     // Single precision angle + sincos is measured to be good enough (2013-11-20).
     // Note that we do the interpolation in double precision still.
     // We can afford to if we keep this kernel memory bound.
-    const float phi = phiAtBegin  * (1.0 - blockOffset) + phiAfterEnd * blockOffset;
-    //if (phi>1e3) {
-    // printf("%d %d %e %e ",time,channel,blockOffset,phi);
-    //}
-    return sincos_f2f_select(-phi, pol);
+    const float phi = phiGradient  * blockOffset;
+    return sincos_f2f_select(phi, pol);
 }
-#define FACTOR(time) correction_factor(time, phiAtBegin, phiAfterEnd, ri, pol, channel)
+#define FACTOR(time) correction_factor(time, phiGradient, pol, channel)
 #endif /* DOPPLER_CORRECTION */
 
 /*!
@@ -221,8 +217,16 @@ __global__ void FIR_filter( void *filteredDataPtr,
   //
   // We need to undo the delay, so we rotate BACK, resulting in a negative constant factor.
   // Note we use the sample frequency=CLOCK_MHZ/1024 (MHz) to normalize the frequency
-  const double phiAtBegin = -2.0 * M_PI * subbandFrequency * delayAtBegin/(CLOCK_MHZ*1e6/1024.0)- phase0;
-  const double phiAfterEnd = -2.0 * M_PI * subbandFrequency * delayAfterEnd/(CLOCK_MHZ*1e6/1024.0)- phase0;
+  // Let tau_0: delay at begin, tau_1: delay after end 
+  // t : time within a block, T: total time (duration) of block,
+  // So, t in [0,T] and  t/T=blockoffset in [0,1]
+  // then delay at t = tau_1 (t/T) + tau_0 (1 - t/T)
+  // exponent at frqeuency f = -j 2 pi (f tau) 
+  // simplifying, exponent = -j 2 pi f tau_0 - j 2 pi f (tau_1 -tau_0) (t/T)
+  // We only need the term that changes with t, so discard the rest
+  // and only keep : 2 pi f (tau_1 - tau_0) (t/T) for the rest of calculations
+  // also replace f with f/f_s where f_s is sample frequency = clock/1024
+  const double phiGradient = 2.0 * M_PI * (subbandFrequency / (CLOCK_MHZ*1e6/1024.0) )*( delayAfterEnd - delayAtBegin );
 #endif
 
   //# const float16 weights = (*weightsData)[channel];