function [dataOut, fdTarget, samplePosTarget] = suqareShapeFreqEst(dataIn, interp, fs, fmin, fmax, fd)
%suqareShapeFreqEst - 利用星座图的集中程度在一定频率范围内使用频率步进估计方波成形信号的频偏
%
% data - input complex data, row vector
% dataLen - data length
% interp - interpolation rate
% fs - sampling rate, in Herts
% fmin - minimum estimation frequency, in Herts
% fmax - maximum estimation frequency, in Herts
% fd - estimation frequency step, in Herts
  
% dimension
[dataLen,channelNum] = size(dataIn);

freqList = fmin:fd:fmax; 
freqNum = length(freqList);
varRecord = zeros(freqNum,interp);

dataOut = zeros(dataLen/interp,channelNum);
fdTarget = zeros(channelNum,1);
samplePosTarget = zeros(channelNum,1);

for channel = 1:channelNum

  for freq = 1:freqNum

    % carrier
    cxr = exp(-1i*2*pi*freqList(freq)/fs*(0:dataLen-1)');

    % apply frequency shift onto input data
    seqRx = dataIn(:,channel).*cxr;
  
    for samplePos = 0:interp-1
      seqRxDown = downsample(seqRx,interp,samplePos); % 下采样
      seqRxDownQuaterOne = abs(real(seqRxDown)) + 1i*abs(imag(seqRxDown)); % 集中到第一象限
      varRecord(freq,samplePos+1) = var(seqRxDownQuaterOne);
    end
  
  end

  [fdIndex, samplePosIndex] = find(varRecord==min(min(varRecord)));
  fdTarget(channel) = freqList(fdIndex);
  samplePosTarget(channel) = samplePosIndex-1;

  cxrTarget = exp(-1i*2*pi*fdTarget(channel)/fs*(0:dataLen-1)');
  seqRxTarget = dataIn(:,channel).*cxrTarget;
  dataOut(:,channel) = downsample(seqRxTarget,interp,samplePosTarget(channel));
end

end