function shift = pmces_batch(ctl,maxShift,fs)
%PMCES Dr. Robert J. Dick's 1981 SSB detector (BATCH version)
% shifts = PMCES(inCtl,maxShift,fs) detects a single SSB shift for all 
% input files specified in inCtl, compensates for the shift, and then 
% writes a .wav output of the compensated audio. shift is the SSB shift
% that was detected.
%
% inCtl is a list of input and output files, comma separated (in the same
% format as FE control files). Each input file is read and processed, then 
% a single SSB shift estimate is produced for all the files. Subsequently, 
% the compensated audio is written to the output files on the right side of 
% the comma. If the output directories do not exist, they are created. If 
% the output files already exist, they are overwritten.
%
% maxShift should be set equal to the absolute maximum shift necessary to be
% detected. Setting maxShift larger than necessary decreases the precision.
%
% Input files should be either .wav or .sph uncompressed format with sampling
% rate fs. The file format is automatically deduced from the file extension. 
%
% Reference: "Co-channel interference separation"
% Pattern analysis and recognition corporation
% Dr. Robert J. Dick
% RADC-TR-80-365 Final Technical Report, December 1980
%
% Rome Air Development Center
% Air Force Systems Command
% Griffis Air Force Base, NY 13441
%
% Adapted and implemented by Mark Harvilla
% Carnegie Mellon University
% December 2, 2011

T_wind   = 0.08;
T_shift  = T_wind/2;
%fs       = 16000;
iFFTSize = floor(T_wind*fs);

minFreq  = 55;
maxFreq  = 300;
minPtVal = floor(fs/maxFreq);
maxPtVal = floor(fs/minFreq);

freqRes  = 5;
histVals = -(maxShift-freqRes):freqRes:maxShift;

fid = fopen(ctl);
L1  = textscan(fid,'%s%s','delimiter',',');
fclose(fid);
inList = L1{1}; outList = L1{2};

Hout = 0;

for cnt = 1:length(inList)
  tline = strtrim(inList{cnt});
  fid2  = fopen(tline,'r');
  if strcmp(tline(end-2:end),'wav') == 1 % <-- MS-WAV format
      fseek(fid2, 44, 'bof');
  else % <-- else assumed to be SPHERE format
      fseek(fid2, 1024, 'bof');
  end
  x_in = fread(fid2, 'int16');
  x_in = x_in/32768;
  fclose(fid2);
  
  X = extract_spec(x_in,T_wind,T_shift,fs,iFFTSize);
 
  iNumWindows = size(X,2);
  cCorr       = zeros(iFFTSize,iNumWindows);
  %cCorr       = zeros(iFFTSize/2 + 1,iNumWindows);
  F0          = zeros(iNumWindows,1);
  df          = zeros(iNumWindows,1);
  fundShift   = zeros(iNumWindows,1);
  %mV1         = zeros(iNumWindows,1);
  %nums        = zeros(1,10000);
 
  manHist = zeros(1,length(histVals));
  for k = 1:iNumWindows
      X_frame    = X(:,k);
      X_frame    = [X_frame;zeros(iFFTSize-size(X,1),1)];
      cCorr(:,k) = ifft(X_frame);
      corrMag    = abs(cCorr(:,k));
      corrMagSqr = corrMag.^2;
      corrRe     = real(cCorr(:,k));
      corrIm     = imag(cCorr(:,k));
      [mV,mInd]  = max(corrMagSqr(minPtVal:maxPtVal));
     
      mInd = mInd + minPtVal - 1;
     
      if mInd == minPtVal || mInd == maxPtVal %if endpoints, don't interpolate
        F0(k) = fs/mInd;
      
        C_Tb_I = corrIm(mInd);
        C_Tb_R = corrRe(mInd);
      else %if not endpoints, interpolate
        mV1    = corrMagSqr(mInd-1); 
        mV2    = corrMagSqr(mInd+1);
        C1     = 0.5*(mV2 - mV1); 
        C2     = 0.5*(mV1 + mV2) - mV;
        avgInd = mInd - C1/(2*C2);
        F0(k)  = fs/avgInd;
 
        R1     = 0.5*(corrRe(mInd+1) - corrRe(mInd-1));
        R2     = 0.5*(corrRe(mInd-1) + corrRe(mInd+1)) - corrRe(mInd);
        I1     = 0.5*(corrIm(mInd+1) - corrIm(mInd-1));
        I2     = 0.5*(corrIm(mInd-1) + corrIm(mInd+1)) - corrIm(mInd);
        C_Tb_R = corrRe(mInd) + (avgInd - mInd)*(R1 + (avgInd - mInd)*R2);
        C_Tb_I = corrIm(mInd) + (avgInd - mInd)*(I1 + (avgInd - mInd)*I2);
      end
      d     = atan2(-C_Tb_I,C_Tb_R);
      df(k) = d*F0(k)/(2*pi);
     
      fundShift(k) = df(k);
      nums1        = fundShift(k):F0(k):maxShift;
      nums2        = fliplr((fundShift(k)-F0(k)):-F0(k):-maxShift);
      nums         = [nums2,nums1];
     
      indx            = quantiz(nums,-(maxShift-freqRes):freqRes:maxShift-freqRes);
      manHist(indx+1) = manHist(indx+1) + mV;
  %     A = (k == 1) || (k == iNumWindows);
  %     if mod(k,50) == 0 || A
  %      figure,bar(histVals,manHist);
  %     end
  end
  Hout = Hout + manHist;
end
[~,HmaxInd] = max(Hout);
fshift      = histVals(HmaxInd);
shift       = -fshift;

for cnt = 1:length(inList)
  tline = strtrim(inList{cnt});
  fid2  = fopen(tline,'r');

  outFile  = strtrim(outList{cnt});
  file_dir = fileparts(outFile);
  if ~exist(file_dir,'dir')
        mkdir(file_dir);
  end

  if strcmp(tline(end-2:end),'wav') == 1 % <-- MS-WAV format
      fseek(fid2, 44, 'bof');
  else % <-- else assumed to be SPHERE format
      fseek(fid2, 1024, 'bof');
  end
  x_in = fread(fid2, 'int16');
  x_in = x_in/32768;
  fclose(fid2);

  y = freqshift(x_in,fs,fshift);
  wavwrite(0.9*y/max(abs(y)),fs,outFile);
end
