RasterProcessTool/LAMPSARProcessProgram/ToolBox/SimulationSAR/GPURFPCKernel.cu

599 lines
19 KiB
Plaintext
Raw Normal View History

#include <iostream>
#include <memory>
#include <cmath>
#include <complex>
#include <device_launch_parameters.h>
#include <cuda_runtime.h>
#include <cublas_v2.h>
#include <cuComplex.h>
#include "BaseConstVariable.h"
#include "GPURFPC.cuh"
#ifdef __CUDANVCC___
__device__ float GPU_getSigma0dB(CUDASigmaParam param, float theta) {//<2F><><EFBFBD><EFBFBD>ֵ
float sigma = param.p1 + param.p2 * exp(-param.p3 * theta) + param.p4 * cos(param.p5 * theta + param.p6);
return sigma;
}
__device__ CUDAVectorEllipsoidal GPU_SatelliteAntDirectNormal(
float RstX, float RstY, float RstZ,
float antXaxisX, float antXaxisY, float antXaxisZ,
float antYaxisX, float antYaxisY, float antYaxisZ,
float antZaxisX, float antZaxisY, float antZaxisZ,
float antDirectX, float antDirectY, float antDirectZ
) {
CUDAVectorEllipsoidal result{ 0,0,-1 };
// <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
float Xst = -1 * RstX; // <20><><EFBFBD><EFBFBD> --> <20><><EFBFBD><EFBFBD>
float Yst = -1 * RstY;
float Zst = -1 * RstZ;
float AntXaxisX = antXaxisX;
float AntXaxisY = antXaxisY;
float AntXaxisZ = antXaxisZ;
float AntYaxisX = antYaxisX;
float AntYaxisY = antYaxisY;
float AntYaxisZ = antYaxisZ;
float AntZaxisX = antZaxisX;
float AntZaxisY = antZaxisY;
float AntZaxisZ = antZaxisZ;
// <20><>һ<EFBFBD><D2BB>
float RstNorm = sqrtf(Xst * Xst + Yst * Yst + Zst * Zst);
float AntXaxisNorm = sqrtf(AntXaxisX * AntXaxisX + AntXaxisY * AntXaxisY + AntXaxisZ * AntXaxisZ);
float AntYaxisNorm = sqrtf(AntYaxisX * AntYaxisX + AntYaxisY * AntYaxisY + AntYaxisZ * AntYaxisZ);
float AntZaxisNorm = sqrtf(AntZaxisX * AntZaxisX + AntZaxisY * AntZaxisY + AntZaxisZ * AntZaxisZ);
float Rx = Xst / RstNorm;
float Ry = Yst / RstNorm;
float Rz = Zst / RstNorm;
float Xx = AntXaxisX / AntXaxisNorm;
float Xy = AntXaxisY / AntXaxisNorm;
float Xz = AntXaxisZ / AntXaxisNorm;
float Yx = AntYaxisX / AntYaxisNorm;
float Yy = AntYaxisY / AntYaxisNorm;
float Yz = AntYaxisZ / AntYaxisNorm;
float Zx = AntZaxisX / AntZaxisNorm;
float Zy = AntZaxisY / AntZaxisNorm;
float Zz = AntZaxisZ / AntZaxisNorm;
float Xant = (Rx * Yy * Zz - Rx * Yz * Zy - Ry * Yx * Zz + Ry * Yz * Zx + Rz * Yx * Zy - Rz * Yy * Zx) / (Xx * Yy * Zz - Xx * Yz * Zy - Xy * Yx * Zz + Xy * Yz * Zx + Xz * Yx * Zy - Xz * Yy * Zx);
float Yant = -(Rx * Xy * Zz - Rx * Xz * Zy - Ry * Xx * Zz + Ry * Xz * Zx + Rz * Xx * Zy - Rz * Xy * Zx) / (Xx * Yy * Zz - Xx * Yz * Zy - Xy * Yx * Zz + Xy * Yz * Zx + Xz * Yx * Zy - Xz * Yy * Zx);
float Zant = (Rx * Xy * Yz - Rx * Xz * Yy - Ry * Xx * Yz + Ry * Xz * Yx + Rz * Xx * Yy - Rz * Xy * Yx) / (Xx * Yy * Zz - Xx * Yz * Zy - Xy * Yx * Zz + Xy * Yz * Zx + Xz * Yx * Zy - Xz * Yy * Zx);
// <20><><EFBFBD><EFBFBD>theta <20><> phi
float Norm = sqrtf(Xant * Xant + Yant * Yant + Zant * Zant); // <20><><EFBFBD><EFBFBD> pho
float ThetaAnt = acosf(Zant / Norm); // theta <20><> Z<><5A><EFBFBD>ļн<C4BC>
float PhiAnt = atanf(Yant / Xant); // -pi/2 ~pi/2
if (abs(Yant) < PRECISIONTOLERANCE) { // X<><58><EFBFBD><EFBFBD>
PhiAnt = 0;
}
else if (abs(Xant) < PRECISIONTOLERANCE) { // Y<><59><EFBFBD>ϣ<EFBFBD>ԭ<EFBFBD><D4AD>
if (Yant > 0) {
PhiAnt = PI / 2;
}
else {
PhiAnt = -PI / 2;
}
}
else if (Xant < 0) {
if (Yant > 0) {
PhiAnt = PI + PhiAnt;
}
else {
PhiAnt = -PI + PhiAnt;
}
}
else { // Xant>0 X <20><><EFBFBD><EFBFBD>
}
if (isnan(PhiAnt)) {
printf("V=[%f,%f,%f];norm=%f;thetaAnt=%f;phiAnt=%f;\n", Xant, Yant, Zant, Norm, ThetaAnt, PhiAnt);
}
result.theta = ThetaAnt;
result.phi = PhiAnt;
result.Rho = Norm;
return result;
}
__device__ float GPU_BillerInterpAntPattern(float* antpattern,
float starttheta, float startphi, float dtheta, float dphi,
long thetapoints, long phipoints,
float searththeta, float searchphi) {
float stheta = searththeta;
float sphi = searchphi;
if (stheta > 90) {
return 0;
}
else {}
float pthetaid = (stheta - starttheta) / dtheta;//
float pphiid = (sphi - startphi) / dphi;
long lasttheta = floorf(pthetaid);
long nextTheta = lasttheta + 1;
long lastphi = floorf(pphiid);
long nextPhi = lastphi + 1;
if (lasttheta < 0 || nextTheta < 0 || lastphi < 0 || nextPhi < 0 ||
lasttheta >= thetapoints || nextTheta >= thetapoints || lastphi >= phipoints || nextPhi >= phipoints)
{
return 0;
}
else {
float x = stheta;
float y = sphi;
float x1 = lasttheta * dtheta + starttheta;
float x2 = nextTheta * dtheta + starttheta;
float y1 = lastphi * dphi + startphi;
float y2 = nextPhi * dphi + startphi;
float z11 = antpattern[lasttheta * phipoints + lastphi];
float z12 = antpattern[lasttheta * phipoints + nextPhi];
float z21 = antpattern[nextTheta * phipoints + lastphi];
float z22 = antpattern[nextTheta * phipoints + nextPhi];
//z11 = powf(10, z11 / 10); // dB-> <20><><EFBFBD><EFBFBD>
//z12 = powf(10, z12 / 10);
//z21 = powf(10, z21 / 10);
//z22 = powf(10, z22 / 10);
float GainValue = (z11 * (x2 - x) * (y2 - y)
+ z21 * (x - x1) * (y2 - y)
+ z12 * (x2 - x) * (y - y1)
+ z22 * (x - x1) * (y - y1));
GainValue = GainValue / ((x2 - x1) * (y2 - y1));
return GainValue;
}
}
__device__ cuComplex GPU_calculationEcho(float sigma0, float TransAnt, float ReciveAnt,
float localangle, float R, float slopeangle, float Pt, float lamda) {
float amp = Pt * TransAnt * ReciveAnt;
amp = amp * sigma0;
amp = amp / (powf(4 * LAMP_CUDA_PI, 2) * powf(R, 4)); // <20><><EFBFBD><EFBFBD>ǿ<EFBFBD><C7BF>
float phi = (-4 * LAMP_CUDA_PI / lamda) * R;
cuComplex echophi = make_cuComplex(0, phi);
cuComplex echophiexp = cuCexpf(echophi);
cuComplex echo = make_cuComplex(echophiexp.x * amp, echophiexp.y * amp);
return echo;
}
__global__ void CUDA_SatelliteAntDirectNormal(float* RstX, float* RstY, float* RstZ,
float antXaxisX, float antXaxisY, float antXaxisZ,
float antYaxisX, float antYaxisY, float antYaxisZ,
float antZaxisX, float antZaxisY, float antZaxisZ,
float antDirectX, float antDirectY, float antDirectZ,
float* thetaAnt, float* phiAnt
, long len) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < len) {
float Xst = -1 * RstX[idx]; // <20><><EFBFBD><EFBFBD> --> <20><><EFBFBD><EFBFBD>
float Yst = -1 * RstY[idx];
float Zst = -1 * RstZ[idx];
float AntXaxisX = antXaxisX;
float AntXaxisY = antXaxisY;
float AntXaxisZ = antXaxisZ;
float AntYaxisX = antYaxisX;
float AntYaxisY = antYaxisY;
float AntYaxisZ = antYaxisZ;
float AntZaxisX = antZaxisX;
float AntZaxisY = antZaxisY;
float AntZaxisZ = antZaxisZ;
// <20><>һ<EFBFBD><D2BB>
float RstNorm = sqrtf(Xst * Xst + Yst * Yst + Zst * Zst);
float AntXaxisNorm = sqrtf(AntXaxisX * AntXaxisX + AntXaxisY * AntXaxisY + AntXaxisZ * AntXaxisZ);
float AntYaxisNorm = sqrtf(AntYaxisX * AntYaxisX + AntYaxisY * AntYaxisY + AntYaxisZ * AntYaxisZ);
float AntZaxisNorm = sqrtf(AntZaxisX * AntZaxisX + AntZaxisY * AntZaxisY + AntZaxisZ * AntZaxisZ);
float Rx = Xst / RstNorm;
float Ry = Yst / RstNorm;
float Rz = Zst / RstNorm;
float Xx = AntXaxisX / AntXaxisNorm;
float Xy = AntXaxisY / AntXaxisNorm;
float Xz = AntXaxisZ / AntXaxisNorm;
float Yx = AntYaxisX / AntYaxisNorm;
float Yy = AntYaxisY / AntYaxisNorm;
float Yz = AntYaxisZ / AntYaxisNorm;
float Zx = AntZaxisX / AntZaxisNorm;
float Zy = AntZaxisY / AntZaxisNorm;
float Zz = AntZaxisZ / AntZaxisNorm;
float Xant = (Rx * Yy * Zz - Rx * Yz * Zy - Ry * Yx * Zz + Ry * Yz * Zx + Rz * Yx * Zy - Rz * Yy * Zx) / (Xx * Yy * Zz - Xx * Yz * Zy - Xy * Yx * Zz + Xy * Yz * Zx + Xz * Yx * Zy - Xz * Yy * Zx);
float Yant = -(Rx * Xy * Zz - Rx * Xz * Zy - Ry * Xx * Zz + Ry * Xz * Zx + Rz * Xx * Zy - Rz * Xy * Zx) / (Xx * Yy * Zz - Xx * Yz * Zy - Xy * Yx * Zz + Xy * Yz * Zx + Xz * Yx * Zy - Xz * Yy * Zx);
float Zant = (Rx * Xy * Yz - Rx * Xz * Yy - Ry * Xx * Yz + Ry * Xz * Yx + Rz * Xx * Yy - Rz * Xy * Yx) / (Xx * Yy * Zz - Xx * Yz * Zy - Xy * Yx * Zz + Xy * Yz * Zx + Xz * Yx * Zy - Xz * Yy * Zx);
// <20><><EFBFBD><EFBFBD>theta <20><> phi
float Norm = sqrtf(Xant * Xant + Yant * Yant + Zant * Zant); // <20><><EFBFBD><EFBFBD> pho
float ThetaAnt = acosf(Zant / Norm); // theta <20><> Z<><5A><EFBFBD>ļн<C4BC>
float PhiAnt = atanf(Yant / Xant); // -pi/2 ~pi/2
if (abs(Yant) < PRECISIONTOLERANCE) { // X<><58><EFBFBD><EFBFBD>
PhiAnt = 0;
}
else if (abs(Xant) < PRECISIONTOLERANCE) { // Y<><59><EFBFBD>ϣ<EFBFBD>ԭ<EFBFBD><D4AD>
if (Yant > 0) {
PhiAnt = PI / 2;
}
else {
PhiAnt = -PI / 2;
}
}
else if (Xant < 0) {
if (Yant > 0) {
PhiAnt = PI + PhiAnt;
}
else {
PhiAnt = -PI + PhiAnt;
}
}
else { // Xant>0 X <20><><EFBFBD><EFBFBD>
}
if (isnan(PhiAnt)) {
printf("V=[%f,%f,%f];norm=%f;thetaAnt=%f;phiAnt=%f;\n", Xant, Yant, Zant, Norm, ThetaAnt, PhiAnt);
}
//if (abs(ThetaAnt - 0) < PRECISIONTOLERANCE) {
// PhiAnt = 0;
//}
//else {}
thetaAnt[idx] = ThetaAnt * r2d;
phiAnt[idx] = PhiAnt * r2d;
//printf("Rst=[%f,%f,%f];AntXaxis = [%f, %f, %f];AntYaxis=[%f,%f,%f];AntZaxis=[%f,%f,%f];phiAnt=%f;thetaAnt=%f;\n", Xst, Yst, Zst
// , AntXaxisX, AntXaxisY, AntXaxisZ
// , AntYaxisX, AntYaxisY, AntYaxisZ
// , AntZaxisX, AntZaxisY, AntZaxisZ
// , phiAnt[idx]
// , thetaAnt[idx]
//);
}
}
__global__ void CUDA_BillerInterpAntPattern(float* antpattern,
float starttheta, float startphi, float dtheta, float dphi,
long thetapoints, long phipoints,
float* searththeta, float* searchphi, float* searchantpattern,
long len) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < len) {
float stheta = searththeta[idx];
float sphi = searchphi[idx];
float pthetaid = (stheta - starttheta) / dtheta;//
float pphiid = (sphi - startphi) / dphi;
long lasttheta = floorf(pthetaid);
long nextTheta = lasttheta + 1;
long lastphi = floorf(pphiid);
long nextPhi = lastphi + 1;
if (lasttheta < 0 || nextTheta < 0 || lastphi < 0 || nextPhi < 0 ||
lasttheta >= thetapoints || nextTheta >= thetapoints || lastphi >= phipoints || nextPhi >= phipoints)
{
searchantpattern[idx] = 0;
}
else {
float x = stheta;
float y = sphi;
float x1 = lasttheta * dtheta + starttheta;
float x2 = nextTheta * dtheta + starttheta;
float y1 = lastphi * dphi + startphi;
float y2 = nextPhi * dphi + startphi;
float z11 = antpattern[lasttheta * phipoints + lastphi];
float z12 = antpattern[lasttheta * phipoints + nextPhi];
float z21 = antpattern[nextTheta * phipoints + lastphi];
float z22 = antpattern[nextTheta * phipoints + nextPhi];
z11 = powf(10, z11 / 10);
z12 = powf(10, z12 / 10);
z21 = powf(10, z21 / 10);
z22 = powf(10, z22 / 10);
float GainValue = (z11 * (x2 - x) * (y2 - y)
+ z21 * (x - x1) * (y2 - y)
+ z12 * (x2 - x) * (y - y1)
+ z22 * (x - x1) * (y - y1));
GainValue = GainValue / ((x2 - x1) * (y2 - y1));
searchantpattern[idx] = GainValue;
}
}
}
__global__ void CUDA_calculationEcho(float* sigma0, float* TransAnt, float* ReciveAnt,
float* localangle, float* R, float* slopeangle,
float nearRange, float Fs, float Pt, float lamda, long FreqIDmax,
cuComplex* echoArr, long* FreqID,
long len) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < len) {
float r = R[idx];
float amp = Pt * TransAnt[idx] * ReciveAnt[idx];
amp = amp * sigma0[idx];
amp = amp / (powf(4 * LAMP_CUDA_PI, 2) * powf(r, 4)); // <20><><EFBFBD><EFBFBD>ǿ<EFBFBD><C7BF>
// <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>λ
float phi = (-4 * LAMP_CUDA_PI / lamda) * r;
cuComplex echophi = make_cuComplex(0, phi);
cuComplex echophiexp = cuCexpf(echophi);
float timeR = 2 * (r - nearRange) / LIGHTSPEED * Fs;
long timeID = floorf(timeR);
//if (timeID < 0 || timeID >= FreqIDmax) {
// timeID = 0;
// amp = 0;
//}
cuComplex echo = make_cuComplex(echophiexp.x, echophiexp.y);
echoArr[idx] = echo;
FreqID[idx] = timeID;
}
}
__global__ void CUDA_AntPatternInterpGain(float* anttheta, float* antphi, float* gain,
float* antpattern, float starttheta, float startphi, float dtheta, float dphi, int thetapoints, int phipoints, long len) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < len) {
float temptheta = anttheta[idx];
float tempphi = antphi[idx];
float antPatternGain = GPU_BillerInterpAntPattern(antpattern,
starttheta, startphi, dtheta, dphi, thetapoints, phipoints,
temptheta, tempphi);
gain[idx] = antPatternGain;
}
}
__global__ void CUDA_InterpSigma(
long* demcls, float* sigmaAmp, float* localanglearr, long len,
CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < len) {
long clsid = demcls[idx];
float localangle = localanglearr[idx];
CUDASigmaParam tempsigma = sigma0Paramslist[clsid];
if (localangle < 0 || localangle >= LAMP_CUDA_PI / 2) {
sigmaAmp[idx] = 0;
}
else {}
if (abs(tempsigma.p1) < PRECISIONTOLERANCE &&
abs(tempsigma.p2) < PRECISIONTOLERANCE &&
abs(tempsigma.p3) < PRECISIONTOLERANCE &&
abs(tempsigma.p4) < PRECISIONTOLERANCE &&
abs(tempsigma.p5) < PRECISIONTOLERANCE &&
abs(tempsigma.p6) < PRECISIONTOLERANCE
) {
sigmaAmp[idx] = 0;
}
else {
float sigma = GPU_getSigma0dB(tempsigma, localangle);
sigma = powf(10.0, sigma / 10.0);// <20><><EFBFBD><EFBFBD>ɢ<EFBFBD><C9A2>ϵ<EFBFBD><CFB5>
//printf("cls:%d;localangle=%f;sigma0=%f;\n", clsid, localangle, sigma);
sigmaAmp[idx] = sigma;
}
}
}
__global__ void CUDAKernel_RFPC_Caluation_R_Gain(
float antX, float antY, float antZ, // <20><><EFBFBD>ߵ<EFBFBD><DFB5><EFBFBD><EFBFBD><EFBFBD>
float* targetX, float* targetY, float* targetZ, long len, // <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
long* demCls,
float* demSlopeX, float* demSlopeY, float* demSlopeZ, // <20>ر<EFBFBD><D8B1><EFBFBD>ʸ<EFBFBD><CAB8>
float antXaxisX, float antXaxisY, float antXaxisZ, // <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ<EFBFBD><CFB5>X<EFBFBD><58>
float antYaxisX, float antYaxisY, float antYaxisZ,// <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ<EFBFBD><CFB5>Y<EFBFBD><59>
float antZaxisX, float antZaxisY, float antZaxisZ,// <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ<EFBFBD><CFB5>Z<EFBFBD><5A>
float antDirectX, float antDirectY, float antDirectZ,// <20><><EFBFBD>ߵ<EFBFBD>ָ<EFBFBD><D6B8>
float Pt,// <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
double refPhaseRange,
float* TransAntpattern, float Transtarttheta, float Transstartphi, float Transdtheta, float Transdphi, int Transthetapoints, int Transphipoints, // <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>߷<EFBFBD><DFB7><EFBFBD>ͼ
float* ReceiveAntpattern, float Receivestarttheta, float Receivestartphi, float Receivedtheta, float Receivedphi, int Receivethetapoints, int Receivephipoints,//<2F><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>߷<EFBFBD><DFB7><EFBFBD>ͼ
float NearR, float FarR, // <20><><EFBFBD>Χ
CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen,// <20><>ֵͼ
float* factorj, long freqnum,
double* outR, // <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
float* outAmp // <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < len) {
double tx = targetX[idx];
double ty = targetY[idx];
double tz = targetZ[idx];
double RstX = antX - tx; // <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʸ<EFBFBD><CAB8>
double RstY = antY - ty;
double RstZ = antZ - tz;
float slopeX = demSlopeX[idx];
float slopeY = demSlopeY[idx];
float slopeZ = demSlopeZ[idx];
double RstR2 = RstX * RstX + RstY * RstY + RstZ * RstZ;
double RstR = sqrt(RstR2); // ʸ<><CAB8><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
//printf("antX=%f;antY=%f;antZ=%f;targetX=%f;targetY=%f;targetZ=%f;RstR=%.6f;diffR=%.6f;\n",antX,antY,antZ,targetX,targetY,targetZ,RstR, RstR - 9.010858499003178e+05);
if (RstR<NearR || RstR>FarR) {
outAmp[idx] = 0;
outR[idx] = 0;
}
else {
// <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
float slopR = sqrtf(slopeX * slopeX + slopeY * slopeY + slopeZ * slopeZ); //
float dotAB = RstX * slopeX + RstY * slopeY + RstZ * slopeZ;
float localangle = acosf(dotAB / (RstR * slopR)); // <20>ֵ<EFBFBD><D6B5><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
float ampGain = 0;
// <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>߷<EFBFBD><DFB7><EFBFBD>ͼָ<CDBC><D6B8>
CUDAVectorEllipsoidal antVector = GPU_SatelliteAntDirectNormal(
RstX, RstY, RstZ,
antXaxisX, antXaxisY, antXaxisZ,
antYaxisX, antYaxisY, antYaxisZ,
antZaxisX, antZaxisY, antZaxisZ,
antDirectX, antDirectY, antDirectZ
);
if (antVector.Rho > 0) {
// <20><><EFBFBD><EFBFBD><E4B7BD>ͼ
float temptheta = antVector.theta * r2d;
float tempphi = antVector.phi * r2d;
float TansantPatternGain =
GPU_BillerInterpAntPattern(
TransAntpattern,
Transtarttheta, Transstartphi, Transdtheta, Transdphi, Transthetapoints, Transphipoints,
temptheta, tempphi);
// <20><><EFBFBD>շ<EFBFBD><D5B7><EFBFBD>ͼ
float antPatternGain = GPU_BillerInterpAntPattern(
ReceiveAntpattern,
Receivestarttheta, Receivestartphi, Receivedtheta, Receivedphi, Receivethetapoints, Receivephipoints,
temptheta, tempphi);
// <20><><EFBFBD><EFBFBD>
float sigma0 = 0;
{
long clsid = demCls[idx];
//printf("clsid=%d\n", clsid);
CUDASigmaParam tempsigma = sigma0Paramslist[clsid];
if (localangle < 0 || localangle >= LAMP_CUDA_PI / 2) {
sigma0 = 0;
}
else {}
if (abs(tempsigma.p1) < PRECISIONTOLERANCE &&
abs(tempsigma.p2) < PRECISIONTOLERANCE &&
abs(tempsigma.p3) < PRECISIONTOLERANCE &&
abs(tempsigma.p4) < PRECISIONTOLERANCE &&
abs(tempsigma.p5) < PRECISIONTOLERANCE &&
abs(tempsigma.p6) < PRECISIONTOLERANCE
) {
sigma0 = 0;
}
else {
float sigma = GPU_getSigma0dB(tempsigma, localangle);
sigma0 = powf(10.0, sigma / 10.0);// <20><><EFBFBD><EFBFBD>ɢ<EFBFBD><C9A2>ϵ<EFBFBD><CFB5>
}
}
ampGain = TansantPatternGain * antPatternGain;
ampGain = ampGain / (powf(4 * LAMP_CUDA_PI, 2) * powf(RstR, 4)); // <20><><EFBFBD><EFBFBD>ǿ<EFBFBD><C7BF>
outAmp[idx] = ampGain * Pt * sigma0;
outR[idx] = RstR - refPhaseRange;
}
else {
outAmp[idx] = 0;
outR[idx] = 0;
}
}
}
}
__global__ void CUDAKernel_PRF_CalFreqEcho(
double* Rarr, float* ampArr, long pixelcount,
float* factorj, long freqnum,
double dx, double nearR,
cuComplex* PRFEcho, long prfid) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < freqnum) {
float fatorj = factorj[idx];
float phi = 0;
float amptemp = 0;
cuComplex tempfreqEcho = PRFEcho[prfid * freqnum + idx];
for (long i = 0; i < pixelcount; i++) { // <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
//phi = (R = R - (floor(R / lamda) - 1) * lamda)* fatorj; // <20><>λ
float phi = Rarr[i] * factorj[idx]; // <20><>λ
amptemp = ampArr[i];
//printf("amp=%f\n", amptemp);
// Eular; exp(ix)=cos(x)+isin(x)
tempfreqEcho.x = tempfreqEcho.x + amptemp * cos(phi); // ʵ<><CAB5>
tempfreqEcho.y = tempfreqEcho.y + amptemp * sin(phi); // <20>鲿
//printf("freqid=%d;fatorj=%.12f;d_R=%.10f;phi=%.10f;echo=complex(%.5f,%.5f)\n", idx, fatorj, Rarr[i], phi, tempfreqEcho.x, tempfreqEcho.y);
}
PRFEcho[prfid * freqnum + idx] = tempfreqEcho;
}
}
__global__ void CUDAKernel_PRFSumEcho_Rows(
double* Rarr,float* ampArr,long Rows,long Cols,
long startRid,
float* factorj, long freqnum,
cuComplex* freqRowsbuffer, long tempRows
){
long idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < Rows) { // <20><><EFBFBD>л<EFBFBD><D0BB><EFBFBD>
double R = 0;
double tempamp = 0;
float phi = 0;
long rid = idx + startRid;
float factor = 0;
for (long jj = 0; jj < freqnum; jj++) {
tempamp = ampArr[rid * Cols + jj];
R = Rarr[rid * Cols + jj];
for (long ii = 0; ii < freqnum; ii++) {
phi = R * factorj[ii];
freqRowsbuffer[idx * freqnum + ii].x = freqRowsbuffer[idx * freqnum + ii].x + tempamp * cos(phi); // ʵ<><CAB5>
freqRowsbuffer[idx * freqnum + ii].y = freqRowsbuffer[idx * freqnum + ii].y + tempamp * sin(phi); // <20>鲿
}
//freqRowsbuffer[idx * freqnum + ii] = tempfreqEcho;
}
}
}
__global__ void CUDAKernel_PRFSumEcho_Freq(
cuComplex* freqRowsbuffer, long tempRows,long freqnum,
cuComplex* PRFEcho, long prfid
) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < freqnum) { // <20><><EFBFBD>л<EFBFBD><D0BB><EFBFBD>
cuComplex tempfreqEcho = freqRowsbuffer[prfid * freqnum + idx];
cuComplex temp = tempfreqEcho;
for (long ii = 0; ii < tempRows; ii++) { // <20><><EFBFBD>ͻ<EFBFBD><CDBB><EFBFBD>
temp = freqRowsbuffer[ii * freqnum + idx];
tempfreqEcho.x = tempfreqEcho.x + temp.x;
tempfreqEcho.y = tempfreqEcho.y + temp.y;
}
freqRowsbuffer[prfid * freqnum + idx] = tempfreqEcho;
}
}
#endif