diff --git a/GPUTool.cu b/GPUTool.cu index 81aad62..58fa1da 100644 --- a/GPUTool.cu +++ b/GPUTool.cu @@ -19,14 +19,141 @@ #define LAMP_CUDA_PI 3.141592653589793238462643383279 - // 定义参数 __device__ cuComplex cuCexpf(cuComplex x) { float factor = exp(x.x); return make_cuComplex(factor * cos(x.y), factor * sin(x.y)); } - + +// 定义仿真所需参数 + + __device__ float GPU_getSigma0dB(CUDASigmaParam param,float theta) { + return param.p1 + param.p2 * exp(-param.p3 * theta) + param.p4 * cos(param.p5 * theta + param.p6); +} + + __device__ CUDAVector GPU_VectorAB(CUDAVector A, CUDAVector B) { + CUDAVector C; + C.x = B.x - A.x; + C.y = B.y - A.y; + C.z = B.z - A.z; + return C; +} + + __device__ float GPU_VectorNorm2(CUDAVector A) { + return sqrtf(A.x * A.x + A.y * A.y + A.z * A.z); + } + + __device__ float GPU_dotVector(CUDAVector A, CUDAVector B) { + return A.x * B.x + A.y * B.y + A.z * B.z; +} + + __device__ float GPU_CosAngle_VectorA_VectorB(CUDAVector A, CUDAVector B) { + return GPU_dotVector(A, B) / (GPU_VectorNorm2(A)*GPU_VectorNorm2(B)); +} + + __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}; + float Xst = -1 * RstX; // 卫星 --> 地面 + 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; + // 天线指向在天线坐标系下的值 + float Xant = (Xst * (AntYaxisY * AntZaxisZ - AntYaxisZ * AntZaxisY) + Xst * (AntXaxisZ * AntZaxisY - AntXaxisY * AntZaxisZ) + Xst * (AntXaxisY * AntYaxisZ - AntXaxisZ * AntYaxisY)) / (AntXaxisX * (AntYaxisY * AntZaxisZ - AntZaxisY * AntYaxisZ) - AntYaxisX * (AntXaxisY * AntZaxisZ - AntXaxisZ * AntZaxisY) + AntZaxisX * (AntXaxisY * AntYaxisZ - AntXaxisZ * AntYaxisY)); + float Yant = (Yst * (AntYaxisZ * AntZaxisX - AntYaxisX * AntZaxisZ) + Yst * (AntXaxisX * AntZaxisZ - AntXaxisZ * AntZaxisX) + Yst * (AntYaxisX * AntXaxisZ - AntXaxisX * AntYaxisZ)) / (AntXaxisX * (AntYaxisY * AntZaxisZ - AntZaxisY * AntYaxisZ) - AntYaxisX * (AntXaxisY * AntZaxisZ - AntXaxisZ * AntZaxisY) + AntZaxisX * (AntXaxisY * AntYaxisZ - AntXaxisZ * AntYaxisY)); + float Zant = (Zst * (AntYaxisX * AntZaxisY - AntYaxisY * AntZaxisX) + Zst * (AntXaxisY * AntZaxisX - AntXaxisX * AntZaxisY) + Zst * (AntXaxisX * AntYaxisY - AntYaxisX * AntXaxisY)) / (AntXaxisX * (AntYaxisY * AntZaxisZ - AntZaxisY * AntYaxisZ) - AntYaxisX * (AntXaxisY * AntZaxisZ - AntXaxisZ * AntZaxisY) + AntZaxisX * (AntXaxisY * AntYaxisZ - AntXaxisZ * AntYaxisY)); + // 计算theta 与 phi + float Norm = sqrtf(Xant * Xant + Yant * Yant + Zant * Zant); // 计算 pho + float ThetaAnt = acosf(Zant / Norm); // theta 与 Z轴的夹角 + float YsinTheta = Yant / sinf(ThetaAnt); + float PhiAnt = (YsinTheta / abs(YsinTheta)) * acosf(Xant / (Norm * sinf(ThetaAnt))); + result.theta = ThetaAnt; + result.phi = PhiAnt; + result.pho = Norm; + return result; +} + + /** +天线方向图插值方法，以双线性插值算法为基础，由theta与phi组合得到的矩阵图为基础数据，通过插值计算的方法获取目标点的数据。 +其中行是theta、列是phi +*/ + __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; + 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); + 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 r = R; + float amp = Pt * TransAnt * ReciveAnt; + amp = amp * sigma0; + amp = amp / (powf(4 * LAMP_CUDA_PI, 2) * powf(r, 4)); // 反射强度 + float phi = (-4 * LAMP_CUDA_PI / lamda) * r; + cuComplex echophi = make_cuComplex(0, phi); + cuComplex echophiexp = cuCexpf(echophi); + cuComplex echo; + echo.x = echophiexp.x * amp; + echo.y = echophiexp.y * amp; + return echo; + } + __global__ void CUDA_DistanceAB(float* Ax, float* Ay, float* Az, float* Bx, float* By, float* Bz,float *R, long len) { long idx = blockIdx.x * blockDim.x + threadIdx.x; if (idx < len) { @@ -144,6 +271,127 @@ __global__ void CUDA_calculationEcho(float* sigma0, float* TransAnt, float* Reci } } +__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_RTPC_Kernel( + float antPx, float antPy, float antPz,// 天线坐标 + 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* demx, float* demy, float* demz, long* demcls, // 地面坐标 + float* demslopex, float* demslopey, float* demslopez, float* demslopeangle,// 地面坡度 + float* Tantpattern, float Tstarttheta, float Tstartphi, float Tdtheta, float Tdphi, long Tthetapoints, long Tphipoints,// 天线方向图相关 + float* Rantpattern, float Rstarttheta, float Rstartphi, float Rdtheta, float Rdphi, long Rthetapoints, long Rphipoints,// 天线方向图相关 + float lamda, float fs, float nearrange, float Pt, long Freqnumbers, // 参数 + CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen,// 地表覆盖类型-sigma插值对应函数-ulaby + cuComplex* outecho, long* d_echoAmpFID, + long len + ) { + long idx = blockIdx.x * blockDim.x + threadIdx.x; + if (idx < len) { + long clsid = demcls[idx]; + CUDAVector Rs{ antPx,antPy,antPz }; + CUDAVector Rt{ demx[idx],demy[idx],demz[idx] }; + CUDAVector Rst{ Rs.x - Rt.x,Rs.y - Rt.y,Rs.z - Rt.z }; + CUDAVector Vslope{ demslopex[idx],demslopey[idx],demslopez[idx] }; + float R = GPU_VectorNorm2(Rst); // 斜距 + + CUDAVectorEllipsoidal Rtanttheta = GPU_SatelliteAntDirectNormal( // 地面目标在天线的位置 + Rst.x, Rst.y, Rst.z, + antXaxisX, antXaxisY, antXaxisZ, + antYaxisX, antYaxisY, antYaxisZ, + antZaxisX, antZaxisY, antZaxisZ, + antDirectX, antDirectY, antDirectZ); + + float localangle=GPU_CosAngle_VectorA_VectorB(Rst, Vslope); // 距地入射角 + float sigma = GPU_getSigma0dB(sigma0Paramslist[clsid], localangle); + sigma = powf(10.0, sigma / 10.0);// 后向散射系数 + + // 发射方向图 + float transPattern = GPU_BillerInterpAntPattern(Tantpattern, + Tstarttheta, Tstartphi, Tdtheta, Tdphi,Tthetapoints, Tphipoints, + Rtanttheta.theta, Rtanttheta.phi); + + // 接收方向图 + float receivePattern = GPU_BillerInterpAntPattern(Rantpattern, + Rstarttheta, Rstartphi, Rdtheta, Rdphi, Rthetapoints, Rphipoints, + Rtanttheta.theta, Rtanttheta.phi); + // 计算振幅、相位 + float amp = Pt * transPattern * receivePattern * sigma / (powf(4 * LAMP_CUDA_PI, 2) * powf(R, 4)); + 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 >= Freqnumbers) { + timeID = 0; + amp = 0; + + } + else {} + + cuComplex echo; + echo.x = echophiexp.x * amp; + echo.y = echophiexp.y * amp; + outecho[idx] = echo; + d_echoAmpFID[idx] = timeID; + } +} //错误提示 @@ -156,8 +404,10 @@ void checkCudaError(cudaError_t err, const char* msg) { } // 主机参数内存声明 -extern "C" void mallocCUDAHost(void* ptr, long memsize) { +extern "C" void* mallocCUDAHost( long memsize) { + void* ptr; cudaMallocHost(&ptr, memsize); + return ptr; } // 主机参数内存释放 @@ -166,8 +416,10 @@ extern "C" void FreeCUDAHost(void* ptr) { } // GPU参数内存声明 -extern "C" void mallocCUDADevice(void* ptr, long memsize) { +extern "C" void* mallocCUDADevice( long memsize) { + void* ptr; cudaMalloc(&ptr, memsize); + return ptr; } // GPU参数内存释放 @@ -268,12 +520,51 @@ extern "C" void calculationEcho(float* sigma0,float* TransAnt,float* ReciveAnt, echoAmp, FreqID, len); cudaDeviceSynchronize(); +} + + + + + +extern "C" void CUDA_RTPC_SiglePRF( + float antPx, float antPy, float antPZ,// 天线坐标 + 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* demx, float* demy, float* demz, long* demcls, // 地面坐标 + float* demslopex, float* demslopey, float* demslopez, float* demslopeangle,// 地面坡度 + float* Tantpattern, float Tstarttheta, float Tstartphi, float Tdtheta, float Tdphi, long Tthetapoints, long Tphipoints,// 天线方向图相关 + float* Rantpattern, float Rstarttheta, float Rstartphi, float Rdtheta, float Rdphi, long Rthetapoints, long Rphipoints,// 天线方向图相关 + float lamda, float fs, float nearrange, float Pt, long Freqnumbers, // 参数 + CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen,// 地表覆盖类型-sigma插值对应函数-ulaby + cuComplex* outecho, long* d_echoAmpFID, + long len +) { + int blockSize = 256; // 每个块的线程数 + int numBlocks = (len + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小 + // 调用 CUDA 核函数 + CUDA_RTPC_Kernel<<>>( + antPx, antPy, antPZ,// 天线坐标 + antXaxisX, antXaxisY, antXaxisZ, // 天线坐标系 + antYaxisX, antYaxisY, antYaxisZ, // + antZaxisX, antZaxisY, antZaxisZ, + antDirectX, antDirectY, antDirectZ,// 天线指向 + demx, demy, demz, demcls, // 地面坐标 + demslopex, demslopey, demslopez, demslopeangle,// 地面坡度 + Tantpattern, Tstarttheta, Tstartphi, Tdtheta, Tdphi, Tthetapoints, Tphipoints,// 天线方向图相关 + Rantpattern, Rstarttheta, Rstartphi, Rdtheta, Rdphi, Rthetapoints, Rphipoints,// 天线方向图相关 + lamda, fs, nearrange, Pt, Freqnumbers, // 参数 + sigma0Paramslist, sigmaparamslistlen,// 地表覆盖类型-sigma插值对应函数-ulaby + outecho, d_echoAmpFID, + len + ); + cudaDeviceSynchronize(); } - #endif \ No newline at end of file diff --git a/GPUTool.cuh b/GPUTool.cuh index 505fd09..64d38fa 100644 --- a/GPUTool.cuh +++ b/GPUTool.cuh @@ -17,14 +17,33 @@ enum LAMPGPUDATETYPE { }; +extern "C" struct CUDASigmaParam { + float p1; + float p2; + float p3; + float p4; + float p5; + float p6; +}; +extern "C" struct CUDAVector { + float x; + float y; + float z; +}; + +extern "C" struct CUDAVectorEllipsoidal { + float theta; + float phi; + float pho; +}; // GPU 内存函数 -extern "C" void mallocCUDAHost(void* ptr, long memsize); // 主机内存声明 +extern "C" void* mallocCUDAHost( long memsize); // 主机内存声明 extern "C" void FreeCUDAHost(void* ptr); -extern "C" void mallocCUDADevice(void* ptr, long memsize); // GPU内存声明 +extern "C" void* mallocCUDADevice( long memsize); // GPU内存声明 extern "C" void FreeCUDADevice(void* ptr); extern "C" void HostToDevice(void* hostptr, void* deviceptr, long memsize);//GPU 内存数据转移 设备 -> GPU extern "C" void DeviceToHost(void* hostptr, void* deviceptr, long memsize);//GPU 内存数据转移 GPU -> 设备 @@ -38,11 +57,35 @@ extern "C" void Norm_Vector(float* Vx, float* Vy, float* Vz, float* R, long memb extern "C" void cosAngle_VA_AB(float* Ax, float* Ay, float* Az, float* Bx, float* By, float* Bz, float* anglecos, long len); extern "C" void 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); extern "C" void calculationEcho(float* sigma0, float* TransAnt, float* ReciveAnt,float* localangle, float* R, float* slopeangle,float nearRange, float Fs, float pt, float lamda, long FreqIDmax,cuComplex* echoAmp, long* FreqID, long len); +extern "C" void CUDA_RTPC_SiglePRF( + float antPx, float antPy, float antPZ,// 天线坐标 + 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* demx, float* demy, float* demz,long* demcls, // 地面坐标 + float* demslopex, float* demslopey, float* demslopez, float* demslopeangle,// 地面坡度 + float* Tantpattern, float Tstarttheta, float Tstartphi, float Tdtheta, float Tdphi, long Tthetapoints, long Tphipoints,// 天线方向图相关 + float* Rantpattern, float Rstarttheta, float Rstartphi, float Rdtheta, float Rdphi, long Rthetapoints, long Rphipoints,// 天线方向图相关 + float lamda, float fs, float nearrange, float Pt, long Freqnumbers, // 参数 + CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen,// 地表覆盖类型-sigma插值对应函数-ulaby + cuComplex* outecho,long* d_echoAmpFID, + long len +); + + #endif #endif + + + + + + + /** *