doppler_centroid_estimate.m

利用能量均衡法对实测的合成孔径雷达回波数据估计出其多普勒中心

function Doppler_centroid=Doppler_Centroid_Estimate(TempFile_In,SysMem)
%函数Doppler_centroid=Doppler_Centroid_Estimate(TempFile_In)估计多普勒中心
%输入参数：TempFile_In为输入回波文件，未进行处理的文件（EIQ或ECO）,SysMem系统内存
%输出参数：Doppler_centroid为多普勒中心位置
%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

[path,name,ext]=fileparts(TempFile_In);   %获得回波文件的地址、文件名和后缀
fid=fopen(TempFile_In,'r');                              %读取回波头文件信息
fread(fid,[1,18*4+98],'int8');
Head=fread(fid,[1,2],'int32');
fclose(fid);
slicelen=Head(1);rawdatawidth=Head(2);
%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
sub_rang=min(floor(SysMem*16384/rawdatawidth),slicelen); %在内存范围内距离向分段读取点数
sub_rang_Num=ceil(slicelen/sub_rang);                    %按方位向读取时距离向分段数
Hwaitbar=waitbar(0,'多普勒中心估计...');  
Sr_average=zeros(1,rawdatawidth);
for i=1:sub_rang_Num
    %对ECO文件的读取
    if ext=='.eco'
       fid=fopen(TempFile_In,'r');fread(fid,256+2*sub_rang*(i-1),'float32');
       if i==sub_rang_Num
           sub_rang=slicelen-(i-1)*sub_rang;%考虑到分段读取时最后一段的长度和前面的不同，需要重新计算最后一段的长度
       end
       Sr=fread(fid,[2*sub_rang,rawdatawidth],[num2str(2*sub_rang),'*float32'],(2*slicelen-2*sub_rang));
       Sr=complex(Sr(1:2:2*sub_rang,:)-mean(mean(Sr(1:2:2*sub_rang,:))),Sr(2:2:2*sub_rang,:)-mean(mean(Sr(2:2:2*sub_rang,:))));
    %对EIQ文件的读取
    else 
        fid=fopen(TempFile_In,'r');fread(fid,1024+sub_rang*(i-1),'uint8');
        if i==sub_rang_Num
            sub_rang=slicelen-(i-1)*sub_rang;%考虑到分段读取时最后一段的长度和前面的不同，需要重新计算最后一段的长度
        end
        Sr=fread(fid,[2*sub_rang,rawdatawidth],[num2str(sub_rang),'*uint8'],(slicelen-sub_rang)*4);
        Sr_Real=Sr(1:sub_rang,:)-128;Sr_Imag=-1*Sr(sub_rang+1:2*sub_rang,:)-128;%EIQ文件的实部和虚部
        Sr=complex(Sr_Real,Sr_Imag);
    end
    Sr=abs(fftshift(fft(Sr,[],2),2));%回波变换到多普勒域中
    Sr_average=sum(Sr,1);%将频谱进行叠加    
    waitbar(i/sub_rang_Num);                                 %进度条进度
end
close(Hwaitbar);%关闭回波文件和进度条
Sr_average=Sr_average/slicelen;%将频谱求平均
%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
%去除坏值点
for i=2:rawdatawidth
    Dif_Sr(i)=Sr_average(i)-Sr_average(i-1);
end
mean_Dif_Sr=mean(Dif_Sr);
for i=2:rawdatawidth
    if Dif_Sr(i)>2*mean(Dif_Sr)
        Sr_average(i)=Sr_average(i-1);
    end
end
%求多普勒中心
sum_Sr_average=sum(Sr_average);%频谱能量
for i=1:rawdatawidth
   if (sum(Sr_average(1:i))<sum_Sr_average/2)&...
           (sum(Sr_average(1:i+1))>sum_Sr_average/2)
      Doppler_centroid=i;%多普勒中心所在的位置
   end
end
figure;plot(Sr_average);