main.c

C语言版,机载SAR回波仿真程序

/*==================================
基于等效斜距模型的SAR回波信号仿真
==================================*/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>

#include "display.h"
#include "Complex.h"
#include "FFT.h"
#include "Parameters.h"

int NextPow2(int x);
void fftshift(comp* Data, int Num);
void main()
{
	//======================定义局部变量=========================
	int i , j;
	int Na;
	int NaPow = 0;
	int NaFFT;
	int Np,Nr;
	int NrPow = 0;
	int NrFFT;
	int StartGate, EndGate;                  //每条距离线上,开始的距离门和结束的距离门数

	float* pfRangeLine;

	double  D2R;                             //deg到rad转换系数
	double t;                                    //方位向时间
	double tau;                                 //距离向时间
	double Phase;                              //临时相位变量
	double Rt;                                 //t时刻斜距
	double ThetaA;                           //天线方向图方位向夹角
    double GainA;                            //方位向天线增益
	double AzimuthBW;                  //方位向波束宽度
	double RangeBW;                     //距离向波束宽度

	double DeltaR;                          //一个距离们所对应的距离
	double Rmin;                            //最近斜距
	double Rref;                             //参考斜距
	double R0;                              //波束中心过场景中心时刻的斜距
	double Wg;                            //地面成像带宽度
	double Ws;                            //地面成像带对应的斜距长度

	double Br ;                               //距离向信号带宽
	double Lambda ;                       //波长
	double fd0;                               //多普勒中心频率
	double fr0;                               //多普勒调频率
	double Ta;                             //合成孔径时间
	double Ba;                             // 多普勒带宽
	double ResA;                         //方位向分辨率
	double ResR;                        //距离向分辨率
	
	comp** Sig;
    comp* pcRangeGate;
	comp* pcRangeGate1;
	comp* pcRangeLine;

	FILE* pFile;

    //========================初始化参数=======================
	D2R = PI/180;
    ThetaL *= D2R;
	ThetaSquint *=D2R;

	Br = Kr*Tp; 
	Lambda = C/f0;
	R0 = Rc/sin(ThetaSquint); 
	fd0 = -2*Vr*cos(ThetaSquint)/Lambda;
	fr0 = 2*Vr*Vr*sin(ThetaSquint)*sin(ThetaSquint)/Lambda/R0;
	AzimuthBW = Lambda/La;
	RangeBW = Lambda/Wa;
	Wg = R0*RangeBW/cos(ThetaL-RangeBW/2);
	Ws = Wg*sin(ThetaL-RangeBW/2);
	Ta = 0.886*R0*AzimuthBW/Vr/sin(ThetaSquint);
	Ba = fr0 * Ta;
	ResA = Vr/Ba;
	ResR = C/Br/2;

	Na = (int)(Ta*PRF);
	NaPow = NextPow2(Na);
	NaFFT = (int)pow(2, NaPow);
	
	DeltaR = C/2/Fs;
	Np = (int)(Tp * Fs);
	Nr = Np + (int)(Ws/DeltaR);
	NrPow = NextPow2(Nr);
	NrFFT = (int)pow(2, NrPow);

	Rref =R0;
	Rmin = Rref - NrFFT*DeltaR/2;

	//=====================显示相关参数==========================
	printf("===================================================\n");
	printf("基于等效斜视距离模型的Stripmap SAR回波仿真!\n");
	printf("===================================================\n");
	printf("波长为：[%f]\n", Lambda);
	printf("合成孔径时间为：[%f]\n", Ta);
	printf("多普勒中心频率为：[%f]\n", fd0);
	printf("多普勒调频率为：[%f]\n", fr0);
	printf("多普勒带宽为：[%f]\n", Ba);
	printf("脉冲重复频率PRF为：[%f]\n\n",PRF);
	printf("方位分辨率为：[%f]\n", ResA);
	printf("距离向分辨率为：[%f]\n\n", ResR);

	printf("距离向点数NrFFT：%d\n", NrFFT);
	printf("方位向点数NaFFT：%d\n\n", NaFFT);

	printf("参考斜距为:[%f]\n", Rref);
	printf("最近斜距：[%f]\n\n", Rmin);

	//======================输出成像所需参数=======================
	pFile = fopen(OutputPara,"w");
	if (NULL == pFile)
	{
         printf("Fail to Open the File\n");
	}
	fprintf(pFile,"%d\n",NaFFT);
	fprintf(pFile,"%d\n",NrFFT);

	fprintf(pFile,"%f\n",Fs);
	fprintf(pFile,"%f\n",Tp);
	fprintf(pFile,"%f\n",Br);
	fprintf(pFile,"%f\n",Lambda);
	fprintf(pFile,"%f\n",PRF);

	fprintf(pFile,"%f\n",Rmin);
	fprintf(pFile,"%f\n",Rref);
	fprintf(pFile,"%f\n",fr0);
	fprintf(pFile,"%f\n",fd0);
	fprintf(pFile,"%f\n",Ta);

	fclose(pFile);

	//========================回波仿真============================
	//分配内存
	Sig = (comp **)calloc(NaFFT, sizeof(comp *));
	for (i = 0; i < NaFFT; i ++)
	{
		Sig[i] = (comp*)calloc(NrFFT, sizeof(comp));
	}

	pFile = fopen(ViewRange, "w");

	for ( i = 0; i < NaFFT; i++)
	{
		t = (i - NaFFT/2)/PRF;
		Rt = sqrt(R0*R0+(Vr*t)*(Vr*t)-2*R0*Vr*t*cos(ThetaSquint));
		fprintf(pFile,"%f\n",Rt);

		ThetaA = fabs(ThetaSquint - asin(R0*sin(ThetaSquint)/Rt));
		if (ThetaA <= 1e-5)
		{
			GainA = 1.0;
		}
		else if (ThetaA > 0.886*AzimuthBW/2)
		{
			GainA = 0.0;
		}
		else
		{
			GainA = sin(PI*ThetaA/AzimuthBW)/(PI*ThetaA/AzimuthBW);
		}

		StartGate = (int)((Rt-DeltaR*Np/2-Rmin)/DeltaR) ;
		EndGate = (int)((Rt+DeltaR*Np/2-Rmin)/DeltaR) ;

		for ( j = StartGate; j <= EndGate; j++)
		{
			tau = 2*Rmin/C+ j/Fs;
			Phase = -4*PI*Rt/Lambda - PI*Kr*pow((tau - 2*Rt/C),2);

			Sig[i][j].re = Sig[i][j].re + (float)(pow(GainA,2)*cos(Phase));
			Sig[i][j].im = Sig[i][j].im + (float)(pow(GainA,2)*sin(Phase));
		}

		printf("Executed:%f  percent.\r", (float)i*100/(NaFFT-1));
	}

	fclose(pFile);
	printf("SAR Echo Simulation  is Totally Accomplished.\n");
	
	//==================回波数据文件写入======================
	pFile = fopen(EchoData, "wb");
	for (i = 0; i<NaFFT; i++)
	{
		fwrite(Sig[i], sizeof(comp), NrFFT, pFile);
	}
	fclose(pFile);
	printf("完成回波数据文件的写入！\n");


	//==================回波数据相位文件写入======================
	pfRangeLine = (float *)calloc(NrFFT,sizeof(float));
	pFile = fopen(EchoPhase,"wb");
	for (i = 0; i < NaFFT; i++)
	{
		for ( j = 0; j< NrFFT; j++)
		{
			pfRangeLine[j] = arg(Sig[i][j]);
		}
		fwrite(pfRangeLine, sizeof(float), NrFFT, pFile);
	}
	fclose(pFile);
	printf("完成回波数据相位文件的写入！\n");

	display(EchoPhase, PhaseBMP,NrFFT, NaFFT, NaFFT);


	//==================方位谱测试文件写入======================
	pcRangeGate = (comp*)calloc(NaFFT,sizeof(comp));
	pcRangeGate1 = (comp*)calloc(NaFFT,sizeof(comp));
	pFile = fopen(AzimuthSpectrum,"w");
	
	for (j = 0; j<NrFFT; j++)
	{
		for ( i = 0; i < NaFFT; i++)
		{
			pcRangeGate[i]=Sig[i][j];
		}
		FFT(pcRangeGate, NaFFT, 1);
		fftshift(pcRangeGate, NaFFT);
		for (i = 0; i < NaFFT; i++)
		{
			pcRangeGate1[i].re = pcRangeGate1[i].re+pcRangeGate[i].re;
            pcRangeGate1[i].im = pcRangeGate1[i].im+pcRangeGate[i].im;
		}
	}
	for (i = 0; i < NaFFT; i++ )
	{
		fprintf(pFile,"%f\n", cabs1(pcRangeGate1[i]));
	}
	fclose(pFile);

	//==================距离谱测试文件写入======================
	pcRangeLine = (comp*)calloc(NrFFT, sizeof(comp));
	pfRangeLine = (float *)calloc(NrFFT,sizeof(float));
	pFile = fopen(RangeSpectrum,"w");
	for (i = 0; i<NaFFT; i++)
	{
		for ( j = 0; j < NrFFT; j++)
		{
			pcRangeLine[j] = Sig[i][j];
		}
		FFT(pcRangeLine, NrFFT, 1);
		fftshift(pcRangeLine,NrFFT);
		for (j = 0; j < NrFFT; j++)
		{
			pfRangeLine[j] = pfRangeLine[j] + cabs1(pcRangeLine[j]);
		}
	}

	for (i = 0; i < NrFFT; i++)
	{
		fprintf(pFile,"%f\n", pfRangeLine[i]);
	}
	fclose(pFile);
	
	printf("完成距离谱测试文件的写入！\n");
	
	//==================释放内存===========================
	free(pcRangeLine);
	free(pcRangeGate);
	free(pfRangeLine);
	free(pcRangeGate1);
	for (i = 0; i<NaFFT; i++)
	{
		free(Sig[i]);
	}

}


int NextPow2(int x)
{
	int i = 0;
	do 
	{
		i++;
	} while ((0x0001<<i)<x);
	return i;
}

void fftshift(comp *Data, int Num)
{
	comp cTemp;
	int i;
	for (i = 0; i < Num/2; i++)
	{
		cTemp = Data[i];
		Data[i] = Data[i+Num/2];
		Data[i+Num/2] = cTemp;
	}
}