% Programmed by Chanwoo Kim for the IEEETran Speech, Audio, and Langauge Processing and ICASSP 2012
% chanwook@cs.cmu.edu
%
% Adapted and modified by Mark Harvilla for the DARPA RATS project
% mharvill@cs.cmu.edu
%
% November 10, 2013 
function pncc(szOutFileName,szInFileName,feat_name,fs,ndeltas,appDeltas,del_width,cmn,fileType,addDither,pbs,ppn,aad_H,numcep,skipC0)

if strcmp(fileType,'wav') == 1
   ad_x = wavread(szInFileName);
else 
   ad_x = readsph(szInFileName);
end
ad_x = ad_x * 32768;

dLamda_L = 0.999;
dLamda_S = 0.999;

dSampRate   = fs;
dPowerCoeff = 1 / 15;

dFactor = 2.0;

iM = 2;
iN = 4;

dLamda  = 0.999;
dLamda2 = 0.5;
dDelta1 = 1;

dLamda3 = 0.85;
dDelta2 = 0.2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Flags
%
bPreem    = 1;
bSSF      = pbs;
bPowerLaw = 1;

dFrameLen    = 0.0256;  % 25.6 ms window length, which is the default setting in CMU Sphinx
dFramePeriod = 0.010;   % 10 ms frame period

iFFTSize  = 1024;

iNumFilts = size(aad_H,2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Flags
%
%
% Array Queue Ring-buffer
%
global Queue_aad_P;
global Queue_iHead;
global Queue_iTail;
global Queue_iWindow;
global Queue_iNumElem;

Queue_iWindow  = 2 * iM + 1;
Queue_aad_P    = zeros(Queue_iWindow, iNumFilts);
Queue_iHead    = 0;
Queue_iTail    = 0;
Queue_iNumElem = 0;

iFL = floor(dFrameLen    * dSampRate);
iFP = floor(dFramePeriod * dSampRate);

if length(ad_x) < iFL+iFP %%% need at least 2 frames (M. Harvilla, 4/26/12)
    if size(ad_x,2) == 1
            ad_x = [ad_x;zeros(iFL+iFP-length(ad_x),1)];
    else
            ad_x = [ad_x,zeros(1,iFL+iFP-length(ad_x))];
    end      
end

iNumFrames = floor((length(ad_x) - iFL) / iFP) + 1;
iSpeechLen = length(ad_x);

%%%%%%%%%%%%
% add dither to avoid
% artificial zeros (from UTD CRSS)
% (M. Harvilla, 4/26/12)
if addDither == 1
    dither = rand(length(ad_x), 1) + rand(length(ad_x), 1) - 1;
    spow   = std(ad_x);
    ad_x   = ad_x + 1e-6 * spow * dither;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Pre-emphasis using H(z) = 1 - 0.97 z ^ -1
%
if (bPreem == 1)
    ad_x = filter([1 -0.97], 1, ad_x);
end

i_FI     = 0;
i_FI_Out = 0;

if bSSF == 1
    aad_P_Out = zeros(iNumFrames - 2 * iM, iNumFilts);
else
    aad_P_Out = zeros(iNumFrames, iNumFilts);
end

adSumPower = zeros(1, iNumFrames);
aad_P      = zeros(iNumFrames, iNumFilts);
aad_P_Out2 = zeros(iNumFrames, iNumFilts);
ad_Q_Out   = zeros(1, iNumFilts);
ad_QMVAvg2 = zeros(1, iNumFilts); %%% these pre-allocations may help a little (M. Harvilla, 4/21/12)
ad_QMVAvg3 = zeros(1, iNumFilts);
ad_QMVPeak = zeros(1, iNumFilts);
MEAN_POWER = 1e10;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dPeakVal = 4e+07; % 4.0638e+07 --> Mean from WSJ0-si84  (Important!)
%%%%%%%%%%%
dMean = dPeakVal;
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Obtaining the short-time Power P(i, j)
%
for m = 0 : iFP : iSpeechLen  - iFL 
    ad_x_st = ad_x(m + 1 : m + iFL) .* hamming(iFL);
    adSpec  = fft(ad_x_st, iFFTSize);
    ad_X    = abs(adSpec(1: iFFTSize / 2));
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    % Calculating the Power P(i, j)
    %
    aad_P(i_FI + 1, :) = transpose(ad_X.^2)*(aad_H.^2);
    
    if bSSF == 1 %%% DO MEDIUM-TIME POWER BIAS SUBTRACTION
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %
        % Ring buffer (using a Queue)
        %
        if (i_FI < iM) || (i_FI > iNumFrames - iM - 1) %%% INCLUDE FIRST iM AND LAST iM FRAMES, INSTEAD OF TRUNCATING (M. Harvilla, 4/26/12)
            if ppn == 1            
                adSumPower(i_FI + 1) = sum(aad_P(i_FI + 1, :));

                if adSumPower(i_FI_Out + 1) > dMean
                     dMean = dLamda_S * dMean + (1 - dLamda_S) * adSumPower(i_FI_Out + 1);
                else
                     dMean = dLamda_L * dMean + (1 - dLamda_L) * adSumPower(i_FI_Out + 1);
                end
                aad_P_Out2(i_FI + 1, :) = aad_P(i_FI + 1, :) / (dMean)  * MEAN_POWER;
            else
                aad_P_Out2(i_FI + 1, :) = aad_P(i_FI + 1, :);
            end     
        end
        
        if (i_FI >= 2 * iM + 1)
            Queue_poll();
        end
        
        Queue_offer(aad_P(i_FI + 1, :));
        ad_Q = Queue_avg(iNumFilts);
            
        if (i_FI == 2 * iM)
            ad_PBias = (ad_Q) * 0.9;
        end

        if (i_FI >= 2 * iM)  
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %
            % Bias Update
            %
            for i = 1 : iNumFilts,
                if (ad_Q(i) > ad_PBias(i))
                   ad_PBias(i) = dLamda * ad_PBias(i)  + (1 - dLamda) * ad_Q(i);
                else
                   ad_PBias(i) = dLamda2 * ad_PBias(i) + (1 - dLamda2) * ad_Q(i);
                end
            end

            for i = 1 : iNumFilts,
                ad_Q_Out(i) =   max(ad_Q(i) - ad_PBias(i), 0);

                if (i_FI == 2 * iM)
                    ad_QMVAvg2(i)  =  0.9 * ad_Q_Out(i);
                    ad_QMVAvg3(i)  =  ad_Q_Out(i);
                    ad_QMVPeak(i)  =  ad_Q_Out(i);
                end

                if (ad_Q_Out(i) > ad_QMVAvg2(i))
                     ad_QMVAvg2(i) = dLamda * ad_QMVAvg2(i)  + (1 -  dLamda)  *  ad_Q_Out(i);
                else
                     ad_QMVAvg2(i) = dLamda2 * ad_QMVAvg2(i) + (1 -  dLamda2) *  ad_Q_Out(i);
                end

                dOrg =  ad_Q_Out(i);

                ad_QMVAvg3(i) = dLamda3 * ad_QMVAvg3(i);

                if (ad_Q(i) <  dFactor * ad_PBias(i))
                    ad_Q_Out(i) = ad_QMVAvg2(i);
                else
                     if (ad_Q_Out(i) <= dDelta1 *  ad_QMVAvg3(i))
                        ad_Q_Out(i) = dDelta2 * ad_QMVAvg3(i);
                     end
                end
                ad_QMVAvg3(i) = max(ad_QMVAvg3(i),   dOrg);

                ad_Q_Out(i) =  max(ad_Q_Out(i), ad_QMVAvg2(i));
            end
            ad_w = ad_Q_Out ./ max(ad_Q, eps);

            mvgAvgH = (1/(2*iN+1))*ones(2*iN+1,1);  %%% frequency smoothing re-implemented (M. Harvilla, 4/21/12); it's a filtering operation

            mask1 = 1:2*iN+1; %create a mask to account for edge effects (e.g. due to the min & max operations on the indices in the original implementation)
            mask1 = (2*iN+1)./mask1;
            mask2 = 2*iN:-1:1;
            mask2 = 9./mask2;
            mask  = [mask1,ones(1,length(ad_w)-2*iN-1),mask2];

            ad_w_sm = conv(mvgAvgH,ad_w);   %do the filtering
            ad_w_sm = ad_w_sm.*mask;        %apply the edge-effect mask
            ad_w_sm = ad_w_sm(iN+1:end-iN); %extract only the relevant portion

            aad_P_Out(i_FI_Out + 1, :) = ad_w_sm .* aad_P(i_FI - iM + 1, :);

            if ppn == 1
                %adSumPower(i_FI_Out + 1)   = sum(aad_P_Out(i_FI_Out + 1, :));
                adSumPower(i_FI + 1)   = sum(aad_P_Out(i_FI_Out + 1, :));
                if  adSumPower(i_FI_Out + 1) > dMean
                     dMean = dLamda_S * dMean + (1 - dLamda_S) * adSumPower(i_FI_Out + 1);
                else
                     dMean = dLamda_L * dMean + (1 - dLamda_L) * adSumPower(i_FI_Out + 1);
                end
                aad_P_Out(i_FI_Out + 1, :) = aad_P_Out(i_FI_Out + 1, :) / (dMean)  * MEAN_POWER;
            end
            %i_FI_Out = i_FI_Out + 1;
        end

    else % if not SSF (i.e. NO MEDIUM-TIME POWER BIAS SUBTRACTION)
        if ppn == 1        
            adSumPower(i_FI + 1)   = sum(aad_P(i_FI + 1, :));

            if adSumPower(i_FI_Out + 1) > dMean
                 dMean = dLamda_S * dMean + (1 - dLamda_S) * adSumPower(i_FI_Out + 1);
            else
                 dMean = dLamda_L * dMean + (1 - dLamda_L) * adSumPower(i_FI_Out + 1);
            end
            aad_P_Out(i_FI + 1, :) = aad_P(i_FI + 1, :) / (dMean)  * MEAN_POWER;
        else
            aad_P_Out(i_FI + 1, :) = aad_P(i_FI + 1, :);
        end
    end
    i_FI     = i_FI + 1;
    i_FI_Out = i_FI_Out + 1;
end

if bSSF == 1
    aad_P_Out2(1+iM:iNumFrames-iM,:) = aad_P_Out;
    aad_P_Out = aad_P_Out2;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Apply the nonlinearity
%
if bPowerLaw == 1
    aadSpec = aad_P_Out .^ dPowerCoeff;
else
    aadSpec = log(aad_P_Out + eps);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DCT & CEPSTRAL LIFTERING
%
aadDCT = dct(aadSpec')';

if skipC0 == 1
	aadDCT(:, numcep+2:iNumFilts) = [];
	aadDCT(:, 1) = [];
else
	aadDCT(:, numcep+1:iNumFilts) = [];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% CMN
%
if cmn == 1
    for i = 1 : numcep
           aadDCT( :, i ) = aadDCT( : , i ) - mean(aadDCT( : , i));
    end
end

aadDCT = aadDCT';

if appDeltas == 1
    %%%%%%%%%%%%
    % compute deltas
    %
	ncep = size(aadDCT,1);
	final_feat = zeros(ncep*(1+ndeltas),size(aadDCT,2));
	final_feat(1:ncep,:) = aadDCT;
    for nd = 1:ndeltas
    	d1 = deltas(final_feat((nd-1)*ncep+1:nd*ncep,:),del_width);
		final_feat(nd*ncep+1:(nd+1)*ncep,:) = d1;
    end
else
	final_feat = aadDCT;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Writing the feature in Sphinx format
%
% [iM, iN] = size(aadDCT);
% iNumData = iM * iN;
% fid = fopen(szOutFeatFileName, 'wb');
% fwrite(fid, iNumData, 'int32');
% iCount = fwrite(fid, aadDCT(:), 'float32');
% fclose(fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Writing the feature in SRI format
%
writesri(final_feat,feat_name,szOutFileName);
    
end

function [] = Queue_offer(ad_x)
    global Queue_aad_P;
    global Queue_iTail;
    global Queue_iWindow;
    global Queue_iNumElem;
    
    Queue_aad_P(Queue_iTail + 1, :) = ad_x;
    
    Queue_iTail    = mod(Queue_iTail + 1, Queue_iWindow);
    Queue_iNumElem = Queue_iNumElem + 1;
    
    if Queue_iNumElem > Queue_iWindow
       error ('Queue overflow'); 
    end
    
  
end


function [ad_x] = Queue_poll()
    global Queue_aad_P;
    global Queue_iHead;
    global Queue_iWindow;
    global Queue_iNumElem;
    
   
    
    if Queue_iNumElem <= 0
       error ('No elements'); 
    end
    
    
    ad_x =  Queue_aad_P(Queue_iHead + 1, :);
    
    Queue_iHead    = mod(Queue_iHead + 1, Queue_iWindow);
    Queue_iNumElem = Queue_iNumElem - 1;
 
end


function[adMean] = Queue_avg(iNumFilts)

    global Queue_aad_P;
    global Queue_iHead;
    global Queue_iWindow;
    global Queue_iNumElem;
    
    adMean = zeros(1, iNumFilts);

    
    iPos = Queue_iHead;
    
    
    for i = 1 : Queue_iNumElem
        adMean = adMean + Queue_aad_P(iPos + 1 ,: );
        iPos   = mod(iPos + 1, Queue_iWindow);
    end
    
    adMean = adMean / Queue_iNumElem;

end
