sar.cpp

上传一个自己开发的合成孔径雷达（SAR)目标成像仿真模拟软件。

	}
   
}

void CSAR::RCMC()
{
   /* 接收回波变换到距离多普勒域 */
	int wa=1,wr=1,win=1;
    int wap=0,wrp=0,winp=0;
	int na,nr;
	int i;
	double Rat;//SAR与目标距离随SAR运动而改变
	double N_cR0,N_cRat;
	while(wa * 2 <= Na)
	{
		wa *= 2;
		wap++;
	}
    while(wr * 2 <= Nr)
	{
		wr *= 2;
		wrp++;
	}
    while(win * 2 <= (RCM_in + 1)*Nr)
	{
		win *= 2;
		winp++;
	}
	complex<double> *TD = new complex<double>[wa];
    complex<double> *FD = new complex<double>[wa];
    complex<double> *LpTemFD =new complex<double>[win];//用来存放插值后距离向频域信号
    complex<double> *LpTemTD =new complex<double>[win];//用来存放插值后距离向时域信号

	for (nr = 0; nr < Nr;nr++)
	{
		for (na = 0;na < Na;na++) TD[na] = complex<double>(*(LpRecSignalRe + na*Nr + nr),*(LpRecSignalIm + na*Nr + nr));
         FFT(TD, FD,wap);
        for (na = 0;na < Na;na++)
		{
          *(LpRecSignalRe + na*Nr + nr) = FD[na].real();
          *(LpRecSignalIm + na*Nr + nr) = FD[na].imag();
		}
	}
	delete TD;
	delete FD;
    TD = new complex<double>[wr];
    FD = new complex<double>[wr];
    /*计算距离徙动，在插值域通过计算Rt处和Rat处的像素个数偏移差得到*/
     N_cR0 = floor((RCM_in + 1)*Rt/Rbin);

	 for( na =0 ;na < Na/2;na++)
	 {
		 if (na >(Na/2)*(B_d/PRF)) 
			 Rat = sqrt((V_ami*Na/(2*PRF))*(V_ami*Na/(2*PRF)) + Rt*Rt);
		 else 
			  Rat = sqrt( (V_ami*(na)/B_d)*(V_ami*(na)/B_d) +Rt*Rt);   
		 N_cRat = floor( (RCM_in + 1)*Rat/Rbin);
		 for (nr = 0;nr < Nr;nr++) TD[nr] = complex<double>(*(LpRecSignalRe + na*Nr + nr),*(LpRecSignalIm + na*Nr + nr));
         FFT(TD, FD,wrp);
         for (i = 0;i < Nr/2;i++)
		 {
			 LpTemFD[i] = FD[i];
             LpTemFD[RCM_in*Nr + Nr/2 + i] = FD[Nr/2 +i];
		 }
		 for (i =  Nr/2; i< (RCM_in*Nr +Nr/2);i++) LpTemFD[i] = complex<double>(0,0);
         IFFT(LpTemFD,LpTemTD,winp);
		 /*将时域信号徙动进行校正*/
		 for (nr = 0;nr < Nr -int(floor((N_cRat -N_cR0)/(RCM_in+1)));nr++) 
		 {
          *(LpRecSignalRe + na*Nr + nr) = LpTemTD[int(N_cRat -N_cR0) + nr*(RCM_in+1)].real();
          *(LpRecSignalIm + na*Nr + nr) = LpTemTD[int(N_cRat -N_cR0) + nr*(RCM_in+1)].imag(); 
		 }
         for (nr = Nr - int(floor((N_cRat -N_cR0)/(RCM_in+1)));nr < Nr;nr++)
		 {
           *(LpRecSignalRe + na*Nr + nr) = 0;
           *(LpRecSignalIm + na*Nr + nr) = 0;
		 }
		 
		 /*对Na-na进行距离徙动校正*/
         for (nr =0;nr < Nr;nr++) TD[nr] = complex<double>(*(LpRecSignalRe + (Na-1-na)*Nr + nr),*(LpRecSignalIm + (Na-1 - na)*Nr + nr));
         FFT(TD, FD,wrp);
         for (i = 0;i < Nr/2;i++)
		 {
			 LpTemFD[i] = FD[i];
             LpTemFD[RCM_in*Nr + Nr/2 + i] = FD[Nr/2 +i];
		 }
		 for (i =  Nr/2; i< (RCM_in*Nr +Nr/2);i++) LpTemFD[i] = complex<double>(0,0);
         IFFT(LpTemFD,LpTemTD,winp);
		 /*将时域信号徙动进行校正*/
		 for (nr = 0;nr < Nr -int(floor((N_cRat -N_cR0)/(RCM_in+1)));nr++) 
		 {
          *(LpRecSignalRe + (Na-1-na)*Nr + nr) = LpTemTD[int(N_cRat -N_cR0) + nr*(RCM_in+1)].real();
          *(LpRecSignalIm + (Na-1-na)*Nr + nr) = LpTemTD[int(N_cRat -N_cR0) + nr*(RCM_in+1)].imag(); 
		 }
         for (nr = Nr - int(floor((N_cRat -N_cR0)/(RCM_in+1)));nr < Nr;nr++)
		 {
           *(LpRecSignalRe + (Na-1-na)*Nr + nr) = 0;
           *(LpRecSignalIm + (Na-1-na)*Nr + nr) = 0;
		 }

	 }
	 delete TD;
     delete FD;
	 delete LpTemTD;
     delete LpTemFD;


}

void CSAR::EchoAzimuthCompress()
{
  if (LpAizm_Time ==NULL) return;
  int na,nr,i,j;
  double angle;
  int w=1;
  int wp=0;
  int win = 1;
  int winp = 0;
  while(w * 2 <= Na)
	{
		w *= 2;
		wp++;
	}
  while(win * 2 <= (A_in+1)*Na)
	{
		win *= 2;
		winp++;
	}

  complex<double> *TD = new complex<double>[w];
  complex<double> *FD = new complex<double>[w];
  complex<double> *LpTemFilter =new complex<double>[win];//用来存放插值后的方位向匹配滤波器
  complex<double> *LpTemFD =new complex<double>[win];//用来存放插值后的方位向频域信号
  complex<double> *LpTemTD =new complex<double>[win];//用来存放插值后方位向时域信号
  complex<double> *Temp =new complex<double>[win];
  memset(TD,0,sizeof(double)*Na);
  /*产生方位向匹配滤波器*/
  for (na = 0;na < NaEffective;na++)
  {
      angle = -PI*FR*LpAizm_Time[na]*LpAizm_Time[na];
      TD[na] = complex<double>(cos(angle),sin(angle));
  }
    FFT(TD, FD,wp);
	 //下面代码完成匹配滤波器中间插0
   for(i = 0;i<Na/2;i++)
   {
       LpTemFilter[i] = FD[i];
	   LpTemFilter[A_in*Na + Na/2 + i] = FD[Na/2+i];
	}
   for (j = Na/2;j < A_in*Na + Na/2;j++)  LpTemFilter[j] = complex<double>(0,0);//在滤波器中间插入A_in*Nr个0 
   for (nr = 0;nr <Nr; nr++)
   {
       //下面代码完成回波方位向多普勒信号中间插0
	   for (na = 0;na<Na/2;na++)
	   {
	    LpTemFD[na] = complex<double>(LpRecSignalRe[na*Nr+nr], LpRecSignalIm[na*Nr+nr]);
        LpTemFD[A_in*Na+Na/2+na] = complex<double>(LpRecSignalRe[(Na/2+na)*Nr+nr], LpRecSignalIm[(Na/2+na)*Nr+nr]);
	   }
	   for (na = Na/2;na<A_in*Na+Na/2;na++) LpTemFD[na] = complex<double>(0,0);
	   //下面完成多普勒域方位向，通过与方位匹配滤波器点乘
       for ( i = 0;i < (A_in+1)*Na;i++) LpTemFD[i] = LpTemFD[i]*LpTemFilter[i];
	   IFFT(LpTemFD,LpTemTD,winp);
       for (j = int((A_in+1)*(Na-floor(NaEffective/2)));j < (A_in+1)*Na;j++) Temp[j-int((A_in+1)*(Na-floor(NaEffective/2)))] = LpTemTD[j];  
	   for (j = 1;j < int((A_in+1)*(Na-floor(NaEffective/2)));j++) Temp[j + int((A_in+1)*(floor(NaEffective/2)))] = LpTemTD[j];
       for (na = 0;na< Na;na++) 
	   {
         LpRecSignalRe[na*Nr+nr] = Temp[A_in*int(floor(NaEffective/2)) +na].real(); 
         LpRecSignalIm[na*Nr+nr] = Temp[A_in*int(floor(NaEffective/2)) +na].imag();
	   }
    }
   delete TD;
   delete FD;
   delete LpTemFilter;
   delete LpTemFD;
   delete LpTemTD;
   delete Temp;
   
}

void CSAR::CSAParaInit()
{
  int na,nr;
  if (LpAizm_Frequency != NULL) delete LpAizm_Frequency;
     LpAizm_Frequency = new double [Na];
  
  if ( LpAizm_Frequency != NULL)
  {
	  for (na = 0;na < Na;na++) 
	  {
        LpAizm_Frequency[na] = -PRF/2.0 + na*PRF/Na;
	  }
  }

  if (LpRange_Frequency != NULL) delete LpRange_Frequency;
     LpRange_Frequency = new double [Nr];
  if (LpRange_Frequency != NULL) 
  {
     for (nr = 0;nr < Nr;nr++) 
	  {
        LpRange_Frequency[nr] = -FS/2.0 + (nr + 1)*FS/Nr;
	  }
  }

  if (LpVTao != NULL) delete LpVTao;
     LpVTao = new double [Nr];
  if (LpVTao != NULL) 
  {
     for (nr = 0;nr < Nr;nr++) 
	  {
        LpVTao[nr] = ((nr+1 - Nr/2.0)*Rbin + Rt)*2.0/C;
	  }
  } 

    if(LpBeta != NULL) delete LpBeta;
	   LpBeta =  new double [Na];
    if (LpBeta != NULL)
	{
      double Data;
	  for (na = 0;na < Na;na++)
	  {
	  Data = Wavelenth*LpAizm_Frequency[na]/(2*V_ami);
      LpBeta[na] = sqrt(1- Data*Data);
      }
	}

  if( LpKsrcInv != NULL) delete LpKsrcInv;
     LpKsrcInv = new double [Na];
  if(( LpKsrcInv !=NULL)&&(LpBeta != NULL))
  {
     for (na = 0;na < Na;na++) 
	  {
	    LpKsrcInv[na] = 2.0*Wavelenth*Ref*(1-pow(LpBeta[na],2))/C/C/(pow(LpBeta[na],3));
	  }
  }

  if(LpKm !=NULL) delete LpKm;
     LpKm = new double [Na];
  if ((LpKm != NULL)&&(LpKsrcInv !=NULL))
  {
   for (na = 0;na < Na;na++) 
   {
     LpKm[na] = fabs(KR)/(1 + fabs(KR)*LpKsrcInv[na]);
   }
  }


  if (LpVTaoRef != NULL) delete LpVTaoRef;
     LpVTaoRef = new double [Na];
  if ((LpVTaoRef != NULL)&&(LpBeta != NULL)) 
  {   
      for (na = 0;na < Na;na++) 
	  {
	    LpVTaoRef[na] = 2.0*Ref/(LpBeta[na]*C);
	  }
  }
  
}

void CSAR::CSAChirpSacle()
{
    int na,nr;
	double angle;
	complex<double> ChirpScaleData;
	complex<double> ChirpScaleF;
    for (nr = 0;nr< Nr;nr++)
	{
     for (na = 0;na < Na;na++)
	 {
       angle = -PI*LpKm[na]*(1/LpBeta[na]-1)*(LpVTao[nr]-LpVTaoRef[na])*(LpVTao[nr]-LpVTaoRef[na]);
	   ChirpScaleF = complex<double>(cos(angle),sin(angle)); 
       ChirpScaleData = complex<double>(LpRecSignalRe[na*Nr+nr], LpRecSignalIm[na*Nr+nr]); 
       ChirpScaleData = ChirpScaleData*ChirpScaleF;
       LpRecSignalRe[na*Nr+nr] = ChirpScaleData.real();
       LpRecSignalIm[na*Nr+nr] = ChirpScaleData.imag();
	 }
    }
	
}

void CSAR::CSAAzimuthFFT()
{
  int w=1;
  int wp=0;
  int na,nr;
  while(w * 2 <= Na)
  {
	w *= 2;
	wp++;
  }
  complex<double> *TD = new complex<double>[w];
  complex<double> *FD = new complex<double>[w];
  memset(TD,0,sizeof(double)*w);
  for (nr = 0;nr < Nr;nr++)
  {
    for (na = 0;na < Na;na++)
		TD[na]= complex<double>(LpRecSignalRe[na*Nr+nr], LpRecSignalIm[na*Nr+nr]);
        FFT(TD, FD,wp);
    for (na = 0;na < Na;na++)
	{
	 LpRecSignalRe[na*Nr+nr] = FD[na<w/2? na+w/2: na-w/2].real(); //保存方位向FFT变换的实部并作fftshift
     LpRecSignalIm[na*Nr+nr] = FD[na<w/2? na+w/2: na-w/2].imag(); //保存方位向FFT变换的虚部并作fftshift
	}
  }
  delete TD;
  delete FD;

}

void CSAR::CSARangeCompress()
{
  //首先要保证进行方位向FFT变换，变换到二维频域
   int na,nr;
   int w=1;
   int wp=0;
   double phase1,phase2,angle;
   complex<double> CSARangeCompressF;
   while(w * 2 <= Nr)
   {
    w *= 2;
	wp++;
   }
   complex<double> *TD = new complex<double>[w];
   complex<double> *FD = new complex<double>[w];
   memset(FD,0,sizeof(double)*w);
   memset(TD,0,sizeof(double)*w);
   for(na = 0;na<Na;na++)
   {
    phase1 = -PI*(1/fabs(KR) + LpKsrcInv[na])*LpBeta[na];
    phase2 = 4.0*PI*Ref*(1.0/LpBeta[na] -1)/C;
	for(nr = 0;nr<Nr;nr++)
	{
      angle = (phase1*LpRange_Frequency[nr] + phase2)*LpRange_Frequency[nr];
      CSARangeCompressF = complex<double>(cos(angle),sin(angle));
	  TD[nr] = complex<double>(LpRecSignalRe[na*Nr+nr],LpRecSignalIm[na*Nr+nr]);
	  FD[nr<w/2? nr+w/2:nr-w/2] = TD[nr]*CSARangeCompressF;
	}
     IFFT(FD, TD,wp);
	 for(nr = 0;nr<Nr;nr++)
	 {
	  LpRecSignalRe[na*Nr+nr] = TD[nr].real();
      LpRecSignalIm[na*Nr+nr] = TD[nr].imag();
     }
   }
    delete FD;
	delete TD;
   
}

void CSAR::CSAAzimuthCompress()
{
   int na,nr;
   int w=1;
   int wp=0;
   double phase1,phase2,angle;
   complex<double> CSAAzimuthCompressF;
   while(w * 2 <= Na)
   {
    w *= 2;
	wp++;
   }
   complex<double> *TD = new complex<double>[w];
   complex<double> *FD = new complex<double>[w];
   memset(FD,0,sizeof(double)*w);
   for (nr = 0;nr<Nr;nr++)
   {
	 for(na = 0;na < Na;na++)
	 {
	   phase1 = -4.0*PI*(1-LpBeta[na])/Wavelenth;
       phase2 = 4.0*PI*LpKm[na]/C/C*(1/LpBeta[na])*(1/LpBeta[na]-1);
       angle = phase1*LpVTao[nr]*C/2 +phase2*(LpVTao[nr]*C/2 -Ref)*(LpVTao[nr]*C/2 -Ref);
       CSAAzimuthCompressF = complex<double>(cos(angle),sin(angle));
       FD[na] = complex<double>(LpRecSignalRe[na*Nr+nr], LpRecSignalIm[na*Nr+nr]);
	   FD[na]= FD[na]* CSAAzimuthCompressF;
	 }
     IFFT(FD, TD,wp);
     for(na = 0;na < Na;na++)
	 {
      LpRecSignalRe[na*Nr+nr] = TD[na].real();
      LpRecSignalIm[na*Nr+nr] = TD[na].imag();
	 }
   }
   delete TD;
   delete FD;

}

void CSAR::CSARangeFFT()
{
  int w=1;
  int wp=0;
  int na,nr;
  while(w * 2 <= Nr)
  {
	w *= 2;
	wp++;
  }
  complex<double> *TD = new complex<double>[w];
  complex<double> *FD = new complex<double>[w];
  memset(TD,0,sizeof(double)*w);
  for (na = 0;na < Na;na++)
  {
    for (nr = 0;nr < Nr;nr++)
		TD[nr]= complex<double>(LpRecSignalRe[na*Nr+nr], LpRecSignalIm[na*Nr+nr]);
        FFT(TD, FD,wp);
    for (nr = 0;nr < Nr;nr++)
	{
	 LpRecSignalRe[na*Nr+nr] = FD[nr<w/2? nr+w/2: nr-w/2].real(); //保存距离向FFT变换的实部并作fftshift
     LpRecSignalIm[na*Nr+nr] = FD[nr<w/2? nr+w/2: nr-w/2].imag(); //保存距离向FFT变换的虚部并作fftshift
	}
  }
  delete TD;
  delete FD;
}
/* 该函数通过打开面目标图像文件FileName,得到面目标图像的宽高及文件数据*/
/* 输入参数FileName, 输出参数 w,h, lpAreaImage                        */
void CSAR::GetAreaInformation(CString *FileName, int* w, int* h, BYTE **lpAreaImage)
{
    CFile cf;
    BITMAPFILEHEADER bfh;
    BITMAPINFOHEADER bih;
	if(!cf.Open(*FileName,CFile::modeReadWrite))
	{
		AfxMessageBox("Read AreaTarget Failure");
		return;
	}

	cf.Read(&bfh,sizeof(BITMAPFILEHEADER));
	cf.Read(&bih,sizeof(BITMAPINFOHEADER));	
	if(bih.biBitCount!=0x0008)
	{
	    AfxMessageBox("Area target must be  8bit scale bitmap");
	    return;
	}
	*h = bih.biHeight;
	*w = bih.biWidth;
	*lpAreaImage = new BYTE[(*w)*(*h)];
	if (*lpAreaImage == NULL) 
	{
		AfxMessageBox("Can't allocate memory");
		return;
	}
	cf.Seek(1024,CFile::current);/*隔过调色版不读*/   
	cf.ReadHuge(*lpAreaImage,(*w)*(*h));
}