function shifts = pmces(ctl,maxShift,fs)
%PMCES Dr. Robert J. Dick's 1981 SSB detector
% shifts = PMCES(inCtl,maxShift,fs) detects an SSB shift for each input
% file specified in inCtl, compensates for the shift, and then writes a 
% .wav output of the compensated audio. shifts is a vector containing the
% estimated shift for each file.
%
% inCtl is a list of input and output files, comma separated (in the same
% format as FE control files). Each input file is read, an SSB shift is
% detected for that file, and then the compensated audio is written to the
% output file 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}; shifts = zeros(length(inList),1);
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); 
 %}
 [x_in,~] = audioread(tline);
 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
 figure,bar(histVals,manHist);
 axis([-350 350 0 0.6]);
 xlabel('Frequency shift (Hz)','fontsize', 16,'fontweight','bold','fontname','helvetica');
 set(gca,'PlotBoxAspectRatio',[5,3,1],'FontSize',16,'fontweight','bold','fontname','helvetica');
 
 [~,HmaxInd] = max(manHist);
 fshift      = histVals(HmaxInd); shifts(cnt) = -fshift;
 y           = freqshift(x_in,fs,fshift);
 
 audiowrite(outFile,0.9*y/max(abs(y)),fs);
 
 %if mod(cnt,100) == 0
    fprintf('%d files processed\n',cnt);
 %end
end
