SAR自聚焦算法

function fdr = AutoFocusing_Contrast（InputRealFile InputImagFile Range Azimuth FocusDepthPoint INITIAL_STEP MAX_ITERATIVE_NUM FINAL_STEP fdr0f_prf）
%----------------------------最大对比度法自聚焦-----------------------------------%
% 输入：距离压缩后转角存储的I/Q两路数据文件
% 输出：按聚焦深度估计出的多普勒调频率
% 算法：最大对比度法
% 准则：选取每个聚焦深度内对比度最大的距离门为参考信号，搜索使该距离门信号成像结果聚焦比最大的fdr作为估计结果
% 调用本模块需要设定的参数详见源代码
% tic;
ta = （0 : Azimuth-1）*（1/f_prf） ;
T_max = Azimuth*（1/f_prf）;
%--------------------------------进行方位压缩，选取参考信号-------------------------------%
Ndepth = ceil（Range / FocusDepthPoint）;          % 测绘带内的聚焦深度数
% 打开距离压缩后转角存储的I/Q两路数据文件
FidReadReal = fopen（InputRealFile‘r‘）;
FidReadImag = fopen（InputImagFile‘r‘）;
if（（Range - FocusDepthPoint * Ndepth） ~= 0）
    temp_shy = 1;
else
    temp_shy = 0;
end
max_value = zeros（1Ndepth）;
max_pos = zeros（1Ndepth）;
for i = 1 : Ndepth    
    AzimuthRef = conj（exp（j*pi*fdr0（i）* （ ta -T_max/2 ）.^2 ））;
    if（（temp_shy == 1） && （i == Ndepth））
        m_stop = Range - FocusDepthPoint * （Ndepth - 1）;
    else
        m_stop = FocusDepthPoint;
    end
    focus_rate = zeros（1m_stop）;    
    % 取每个聚焦深度内聚焦比最大的距离门作为参考信号
    for m = 1 : m_stop
        TempRe = fread（FidReadReal Azimuth ‘float32‘）;
        TempIm = fread（FidReadImag Azimuth ‘float32‘）;
        Temp = （TempRe + j*TempIm）.‘;
        Hout = Temp .* AzimuthRef;
        out_effective = abs（fft（Hout））;        
        % 计算聚焦比
        focus_rate（m） = var（out_effective） / mean（out_effective）;
   end    
    [max_value（i） max_pos（i）] = max（focus_rate）;
end

%-----------------------------------最大对比度法自聚焦-----------------------------------%

first_validfdr_pos = -1;        % 存放第一个自聚焦成功的聚焦深度，是测绘带内的第几个聚焦深度
update_step = 0;                % 清除步长更新标识
fdr_step = INITIAL_STEP;        % 设定自聚焦估计的初始搜索步长；注意：fdr的步长是在初始步长的基础上，每次减小一倍，直到步长不大于1Hz；
                                % 因此，如果初始步长是2的整数次幂，则fdr的估计精度为1Hz
iterative_num = 1;              % 初始化迭代次数标识
autofocusing_over = 0;          % 初始化自聚焦结束表示
autofocusing_fail = 0;          % 初始化自聚焦失败标识

fdr = zeros（1 Ndepth）;
focus_rate_current_fdr = zeros（13）;
for i = 1 : Ndepth
    % 当前聚焦深度的fdr初值
    current_fdr = fdr0（i）;
    % 读入当前聚焦深度的参考信号
    fseek（FidReadReal （（i - 1） * FocusDepthPoint + max_pos（i） - 1） * Azimuth * 4 -1）;
    fseek（FidReadImag （（i - 1） * FocusDepthPoint + max_pos（i） - 1） * Azimuth * 4 -1）;

    % 得到参考信号频谱
    TempRe = fread（FidReadReal Azimuth ‘float32‘）;
    TempIm = fread（FidReadImag Azimuth ‘float32‘）;
    Temp = （TempRe + j*TempIm）.‘;
% % % % %     HTemp = fft（Temp）;
    
    % 反复迭代，直到估计出当前聚焦深度的fdr
    while（autofocusing_over ~= 1）
        % 反复迭代，直到在当前搜索步长下估计出fdr
        while（update_step ~= 1）
            % 判断当前搜索步长是否为初始步长
            if（fdr_step == INITIAL_STEP）
                % 判断是否为基于初始步长的第一次迭代
                % 是，则需要计算三个聚焦比，即：current_fdr-fdr_step、current_fdr、current_fdr+fdr_step
                if（iterative_num == 1）
                    leftpos_fdr  = -1