MAP

function Exp_Ray_1(isMasterNode)

% Some description of the experiment goes here

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Set up the basic folders

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

expName = 'scrap6_with_eff';

dataFolderPrefix = 'featR';

if isunix

    expFolderPrefix = '/scratch/nrclark/exps/';

else

    expFolderPrefix = 'D:\Exps';

end

expFolder = fullfile(expFolderPrefix,expName);

hmmFolder = fullfile(expFolder,'hmm');

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Sort out the training (LEARNING) condition

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

learnFolder = fullfile(expFolder,'featL');

xL = cJob('L', learnFolder);

xL.participant = 'Normal';

xL.MAPparamChanges= {'DRNLParams.rateToAttenuationFactorProb=0;', 'OMEParams.rateToAttenuationFactorProb=0;' };

xL.MAPparamChanges= {};

xL.noiseLevToUse   =  -200;

xL.noiseLevToUse   =  20;

xL.speechLevToUse  =  50;

xL.MAPopHSR = 1;

xL.MAPopMSR = 0;

xL.MAPopLSR = 0;

xL.numCoeff = 14;

xL.removeEnergyStatic = 0;

%%%%% Group of params that will influence simulation run time %%%%%%%

xL.numWavs = 512; %MAX=8440

testWavs = 48; %MAX = 358

nzLevel = [-200 40 60];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

xL.noisePreDur = 1;

xL.noisePostDur = 0.1;

xL.truncateDur  = xL.noisePreDur-0.1; 

xL.noiseName = '20TalkerBabble';

if isMasterNode && ~isdir(xL.opFolder)

    mkdir(xL.opFolder);

    xL = xL.assignFiles;

    xL.storeSelf; % This is a call to a member function and is not a pointless line of code!

end

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Sort out the testing (RECOGNITION) conditions

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

recConditions = numel(nzLevel);

tmpIdx=0;

for nn = 0*recConditions+1:1*recConditions    

    tmpIdx=tmpIdx+1;

    xR{nn} = xL; %simply copy the "Learn" object and change it a bit below

    recFolder = fullfile(expFolder,[dataFolderPrefix num2str(nn)]);

    xR{nn}.opFolder = recFolder;    

    %These are the interesting differences between training and testing

    xR{nn}.numWavs = testWavs; %MAX = 358

    xR{nn}.noiseLevToUse = nzLevel(tmpIdx);

    %xR{nn}.MAPparamChanges= {'DRNLParams.rateToAttenuationFactorProb=0;'};

    %Now just to wrap it up ready for processing

    if isMasterNode && ~isdir(xR{nn}.opFolder)

        mkdir(xR{nn}.opFolder);

        xR{nn} = xR{nn}.assignWavPaths('R');

        xR{nn} = xR{nn}.assignFiles;

        xR{nn}.storeSelf;

    end

end

% tmpIdx=0;

% for nn = 1*recConditions+1:2*recConditions    

%     tmpIdx=tmpIdx+1;

%     xR{nn} = xL; %simply copy the "Learn" object and change it a bit below

%     recFolder = fullfile(expFolder,[dataFolderPrefix num2str(nn)]);

%     xR{nn}.opFolder = recFolder;    

%     

%     %These are the interesting differences between training and testing

%     xR{nn}.numWavs = testWavs; %MAX = 358

%     xR{nn}.noiseLevToUse = nzLevel(tmpIdx);

%     xR{nn}.MAPparamChanges= {'DRNLParams.rateToAttenuationFactorProb=-10^(-10/20);'};

%     

%     

%     %Now just to wrap it up ready for processing

%     if isMasterNode && ~isdir(xR{nn}.opFolder)

%         mkdir(xR{nn}.opFolder);

%         xR{nn} = xR{nn}.assignWavPaths('R');

%         xR{nn} = xR{nn}.assignFiles;

%         xR{nn}.storeSelf;

%     end

% end

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ** Generate features **

% This is the time consuming, processing intensive portion of the program.

% Nodes that are not the master node are only interested in the opFolder

% member of the jobjects.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

worker(xL.opFolder);

maxConds = nn;

if ~isMasterNode %dont bother wasting master node effort on generating testing features (for now)

    for nn = 1:maxConds

        worker(xR{nn}.opFolder);

    end

end

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Train and test the recogniser - a job for the master node only

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if isMasterNode    

    while(~all(xL.todoStatus==2))        

        disp('Waiting on straggler nodes to complete their jobs before HMM is trained . . .')

        pause(30); %Wait for 30 seconds before looking again

        xL.lockJobList;

        xL = xL.loadSelf; %Reload incase changed

        xL.unlockJobList;

    end

    y = cHMM(hmmFolder);    

    y.numCoeff = (xL.numCoeff-logical(xL.removeEnergyStatic)) * 3;

    y.createSCP(xL.opFolder)

    y.createMLF(xL.opFolder)

    y.train(xL.opFolder) %This node can be busy training, even if other jobs are being processed for testing

    % ALLOW MASTER NODE TO MUCK IN WITH GENERATING TESTING FEATURES ONCE

    % HMM HAS BEEN TRAINED

    for nn = 1:maxConds

        worker(xR{nn}.opFolder);

    end    

    xR{end}.lockJobList;

    xR{end} = xR{end}.loadSelf; %Reload changes

    xR{end}.unlockJobList;

    while(~all(xR{end}.todoStatus==2))        

        disp('Waiting on straggler nodes to complete their jobs before HMM is tested . . .')

        pause(30); %Wait for 30 seconds before looking again

        xR{end}.lockJobList;

        xR{end} = xR{end}.loadSelf; %Reload incase changed

        xR{end}.unlockJobList;

    end

    for nn = 1:maxConds

        y.createSCP(xR{nn}.opFolder);

        y.test(xR{nn}.opFolder);

    end

    %Show all of the scores in the command window at the end

    for nn = 1:maxConds

        y.score(xR{nn}.opFolder);

    end

end

function Exp_Tutorial_1(isMasterNode)

% Some description of the experiment goes here

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Set up the basic folders

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

expName = 'Tutorial';

dataFolderPrefix = 'featR';

if isunix

    expFolderPrefix = '/scratch/nrclark/exps/';

else

    expFolderPrefix = 'D:\Exps';

end

expFolder = fullfile(expFolderPrefix,expName);

hmmFolder = fullfile(expFolder,'hmm');

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Sort out the training (LEARNING) condition

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

learnFolder = fullfile(expFolder,'featL');

xL = cJob('L', learnFolder);

xL.participant = 'Normal';

xL.MAPparamChanges= {'DRNLParams.rateToAttenuationFactorProb=0;', 'OMEParams.rateToAttenuationFactorProb=0;' };

xL.noiseLevToUse   =  -200;

xL.speechLevToUse  =  60;

xL.MAPopHSR = 1;

xL.MAPopMSR = 0;

xL.MAPopLSR = 0;

xL.numCoeff = 14;

xL.removeEnergyStatic = 0;

%%%%% Group of params that will influence simulation run time %%%%%%%

xL.numWavs = 12; %MAX=8440

testWavs = 6; %MAX = 358

nzLevel = [-200 40:10:70];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

xL.noisePreDur = 1;

xL.noisePostDur = 0.1;

xL.truncateDur  = xL.noisePreDur-0.1; 

xL.noiseName = 'pink_demo';

if isMasterNode && ~isdir(xL.opFolder)

    mkdir(xL.opFolder);

    xL = xL.assignFiles;

    xL.storeSelf; % This is a call to a member function and is not a pointless line of code!

end

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Sort out the testing (RECOGNITION) conditions

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

recConditions = numel(nzLevel);

tmpIdx=0;

for nn = 0*recConditions+1:1*recConditions    

    tmpIdx=tmpIdx+1;

    xR{nn} = xL; %simply copy the "Learn" object and change it a bit below

    recFolder = fullfile(expFolder,[dataFolderPrefix num2str(nn)]);

    xR{nn}.opFolder = recFolder;    

    %These are the interesting differences between training and testing

    xR{nn}.numWavs = testWavs; %MAX = 358

    xR{nn}.noiseLevToUse = nzLevel(tmpIdx);

    xR{nn}.MAPparamChanges= {'DRNLParams.rateToAttenuationFactorProb=0;'};

    %Now just to wrap it up ready for processing

    if isMasterNode && ~isdir(xR{nn}.opFolder)

        mkdir(xR{nn}.opFolder);

        xR{nn} = xR{nn}.assignWavPaths('R');

        xR{nn} = xR{nn}.assignFiles;

        xR{nn}.storeSelf;

    end

end

tmpIdx=0;

for nn = 1*recConditions+1:2*recConditions    

    tmpIdx=tmpIdx+1;

    xR{nn} = xL; %simply copy the "Learn" object and change it a bit below

    recFolder = fullfile(expFolder,[dataFolderPrefix num2str(nn)]);

    xR{nn}.opFolder = recFolder;    

    %These are the interesting differences between training and testing

    xR{nn}.numWavs = testWavs; %MAX = 358

    xR{nn}.noiseLevToUse = nzLevel(tmpIdx);

    xR{nn}.MAPparamChanges= {'DRNLParams.rateToAttenuationFactorProb=-10^(-10/20);'};

    %Now just to wrap it up ready for processing

    if isMasterNode && ~isdir(xR{nn}.opFolder)

        mkdir(xR{nn}.opFolder);

        xR{nn} = xR{nn}.assignWavPaths('R');

        xR{nn} = xR{nn}.assignFiles;

        xR{nn}.storeSelf;

    end

end

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ** Generate features **

% This is the time consuming, processing intensive portion of the program.

% Nodes that are not the master node are only interested in the opFolder

% member of the jobjects.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

worker(xL.opFolder);

maxConds = nn;

if ~isMasterNode %dont bother wasting master node effort on generating testing features (for now)

    for nn = 1:maxConds

        worker(xR{nn}.opFolder);

    end

end

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Train and test the recogniser - a job for the master node only

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if isMasterNode    

    while(~all(xL.todoStatus==2))        

        disp('Waiting on straggler nodes to complete their jobs before HMM is trained . . .')

        pause(30); %Wait for 30 seconds before looking again

        xL.lockJobList;

        xL = xL.loadSelf; %Reload incase changed

        xL.unlockJobList;

    end

    y = cHMM(hmmFolder);    

    y.numCoeff = (xL.numCoeff-logical(xL.removeEnergyStatic)) * 3;

    y.createSCP(xL.opFolder)

    y.createMLF(xL.opFolder)

    y.train(xL.opFolder) %This node can be busy training, even if other jobs are being processed for testing

    % ALLOW MASTER NODE TO MUCK IN WITH GENERATING TESTING FEATURES ONCE

    % HMM HAS BEEN TRAINED

    for nn = 1:maxConds

        worker(xR{nn}.opFolder);

    end    

    xR{end}.lockJobList;

    xR{end} = xR{end}.loadSelf; %Reload changes

    xR{end}.unlockJobList;

    while(~all(xR{end}.todoStatus==2))        

        disp('Waiting on straggler nodes to complete their jobs before HMM is tested . . .')

        pause(30); %Wait for 30 seconds before looking again

        xR{end}.lockJobList;

        xR{end} = xR{end}.loadSelf; %Reload incase changed

        xR{end}.unlockJobList;

    end

    for nn = 1:maxConds

        y.createSCP(xR{nn}.opFolder);

        y.test(xR{nn}.opFolder);

    end

    %Show all of the scores in the command window at the end

    for nn = 1:maxConds

        y.score(xR{nn}.opFolder);

    end

end

function Exp_Tutorial_2(isMasterNode)

% This tutorial recycles a HMM

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Set up the basic experiment parameters

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

expName = 'Tutorial';

dataFolderPrefix = 'recycle_featR';

if isunix

    expFolderPrefix = '/scratch/nrclark/exps/';

else

    expFolderPrefix = 'D:\Exps';

end

% expFolderPrefix = pwd;

expFolder = fullfile(expFolderPrefix,expName);

hmmFolder = fullfile(expFolder,'hmm');

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Sort out the training (LEARNING) condition

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

learnFolder = fullfile(expFolder,'featL');

xL = cJob('L', learnFolder);

xL.participant = 'Normal';

xL.MAPparamChanges= {'DRNLParams.rateToAttenuationFactorProb=0;', 'OMEParams.rateToAttenuationFactorProb=0;' };

xL.noiseLevToUse   =  -200;

xL.speechLevToUse  =  60;

xL.MAPopHSR = 1;

xL.MAPopMSR = 0;

xL.MAPopLSR = 0;

xL.numCoeff = 14;

xL.removeEnergyStatic = 0;

%%%%% Group of params that will influence simulation run time %%%%%%%

xL.numWavs = 12; %MAX=8440

testWavs = 6; %MAX = 358

nzLevel = [-200 40:10:70];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

xL.noisePreDur = 1;

xL.noisePostDur = 0.1;

xL.truncateDur  = xL.noisePreDur-0.1; 

xL.noiseName = 'pink_demo';

% if isMasterNode && ~isdir(xL.opFolder)

%     mkdir(xL.opFolder);

%     xL = xL.assignFiles;

%     xL.storeSelf;

% end

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Sort out the testing (RECOGNITION) conditions

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

recConditions = numel(nzLevel);

tmpIdx=0;

for nn = 0*recConditions+1:1*recConditions    

    tmpIdx=tmpIdx+1;

    xR{nn} = xL; %simply copy the "Learn" object and change it a bit below

    recFolder = fullfile(expFolder,[dataFolderPrefix num2str(nn)]);

    xR{nn}.opFolder = recFolder;    

    %These are the interesting differences between training and testing

    xR{nn}.numWavs = testWavs; %MAX = 358

    xR{nn}.noiseLevToUse = nzLevel(tmpIdx); 

    xR{nn}.MAPparamChanges= {'DRNLParams.a=400;'};

    %Now just to wrap it up ready for processing

    if isMasterNode && ~isdir(xR{nn}.opFolder)

        mkdir(xR{nn}.opFolder);

        xR{nn} = xR{nn}.assignWavPaths('R');

        xR{nn} = xR{nn}.assignFiles;

        xR{nn}.storeSelf;

    end

end

tmpIdx=0;

for nn = 1*recConditions+1:2*recConditions    

    tmpIdx=tmpIdx+1;

    xR{nn} = xL; %simply copy the "Learn" object and change it a bit below

    recFolder = fullfile(expFolder,[dataFolderPrefix num2str(nn)]);

    xR{nn}.opFolder = recFolder;    

    %These are the interesting differences between training and testing

    xR{nn}.numWavs = testWavs; %MAX = 358

    xR{nn}.noiseLevToUse = nzLevel(tmpIdx); 

    xR{nn}.MAPparamChanges= {'DRNLParams.a=400;'};

    xR{nn}.mainGain = [27.2013;   26.0797;   26.0939;   26.7997;   26.0520];  

    xR{nn}.TCdBO    = [37;   37;   37;   37;   37];      %Compression thresholds (in dB OUTPUT from 2nd filt)

    xR{nn}.TMdBO    = [20;   20;   20;   20;   20];      %MOC thresholds (in dB OUTPUT from 2nd filt)

    xR{nn}.ARthresholddB = 85;       % dB SPL (input signal level) =>200 to disable

    xR{nn}.MOCtau = 1;

    xR{nn}.useAid = 1;

    %Now just to wrap it up ready for processing

    if isMasterNode && ~isdir(xR{nn}.opFolder)

        mkdir(xR{nn}.opFolder);

        xR{nn} = xR{nn}.assignWavPaths('R');

        xR{nn} = xR{nn}.assignFiles;

        xR{nn}.storeSelf;

    end

end

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ** Generate features **

% This is the time consuming, processing intensive portion of the program.

% Nodes that are not the master node are only interested in the opFolder

% member of the jobjects.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% worker(xL.opFolder);

maxConds = nn;

if ~isMasterNode %dont bother wasting master node effort on generating testing features (for now)

    for nn = 1:maxConds

        worker(xR{nn}.opFolder);

    end

end

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Train and test the recogniser - a job for the master node only

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if isMasterNode    

%     while(~all(xL.todoStatus==2))        

%         disp('Waiting on straggler nodes to complete their jobs before HMM is trained . . .')

%         pause(30); %Wait for 30 seconds before looking again

%         xL.lockJobList;

%         xL = xL.loadSelf; %Reload incase changed

%         xL.unlockJobList;

%     end

    y = cHMM(hmmFolder);    

    y.numCoeff = 14*3;

%     y.createSCP(xL.opFolder)

%     y.createMLF(xL.opFolder)

%     y.train(xL.opFolder) %This node can be busy training, even if other jobs are being processed for testing

    % ALLOW MASTER NODE TO MUCK IN WITH GENERATING TESTING FEATURES ONCE

    % HMM HAS BEEN TRAINED

    for nn = 1:maxConds

        worker(xR{nn}.opFolder);

    end    

    xR{end}.lockJobList;

    xR{end} = xR{end}.loadSelf; %Reload changes

    xR{end}.unlockJobList;

    while(~all(xR{end}.todoStatus==2))        

        disp('Waiting on straggler nodes to complete their jobs before HMM is tested . . .')

        pause(30); %Wait for 30 seconds before looking again

        xR{end}.lockJobList;

        xR{end} = xR{end}.loadSelf; %Reload incase changed

        xR{end}.unlockJobList;

    end

    for nn = 1:maxConds

        y.createSCP(xR{nn}.opFolder);

        y.test(xR{nn}.opFolder);

    end

    %Show all of the scores in the command window at the end

    for nn = 1:maxConds

        y.score(xR{nn}.opFolder);

    end

end

% This function wraps up whatever version of MAP I want to call. It is

% implemented partly because I want to avoid messing with jobject too much

% and mostly because I dont want to declare globals in my class.

function [myANprobRateOutput, mydt, myBF] = MAPwrap(stimulus, sampleRate, BFlist, participant, AN_spikesOrProbability, paramChanges)

global ANprobRateOutput  dt BFlist

% disp(20*log10(sqrt(mean(stimulus.^2))/20e-6))

MAP1_14(stimulus, sampleRate, BFlist, participant, AN_spikesOrProbability, paramChanges);

% disp(20*log10(sqrt(mean(stimulus.^2))/20e-6))

myANprobRateOutput   = ANprobRateOutput;

mydt = dt;

myBF = BFlist;

classdef cEssexAid

    %ESSEXAID_WRAPCLASS Wrapper for the EssexAid - Nick Clark July 2011

    %   This class wraps up the EssexAid algorithm function that processes

    %   each block of samples. This wrapper closely emulates the GUI used

    %   in the lab and runs stimuli through the exact same algorithm used

    %   in the lab. It even includes a helper function to generate C code

    %   from the algorithm for use in a real-time framework.

    %% *********************************************************

    %  properties                      _   _

    %                                 | | (_)

    %  _ __  _ __ ___  _ __   ___ _ __| |_ _  ___  ___

    % | '_ \| '__/ _ \| '_ \ / _ \ '__| __| |/ _ \/ __|

    % | |_) | | | (_) | |_) |  __/ |  | |_| |  __/\__ \

    % | .__/|_|  \___/| .__/ \___|_|   \__|_|\___||___/

    % | |             | |

    % |_|             |_|

    %************************************************************

    %% **********************************************************

    % Public properties - can be set by user

    %************************************************************

    properties(Access = public)

        sr         = 48e3;

        numSamples = 1024; %MAX=6912, LAB_USE=48

        stimulusUSER                   

        %------------------------------------------------------------------

        % Params for audiometric freqs 250, 500, 1000, 2000, 4000, 8000 Hz

        %------------------------------------------------------------------

        audiometry_dB= [ 0;    0;    0;    0;    0;   0];   %Pure tone threshold in dB SPL

        mainGain_dB  = [ 0;    0;    0;    0;    0;   0];   %Gain applied at audiometric frequencies

        TC_dBHL      = [40;   40;   40;   40;   40;  40];   %Compression thresholds (in dB HL from 2nd filt)

        TM_dBHL      = [10;   10;   10;   10;   10;  10];   %MOC thresholds (in dB OUTPUT from 2nd filt)

        DRNLc        = [ 0.2;  0.2;  0.2;  0.2;  0.2; 0.2]; %Compression exponent at audiometric frequencies

        %------------------------------------------------------------------

        % Dynamic compression properties

        %------------------------------------------------------------------

        ARtau = 60e-3;       %decay time constant

        ARthreshold_dB = 85; %dB SPL (input signal level) =>200 to disable

        MOCtau = 450e-3;     %Time constant in Seconds

        MOCfactor = 0.5;     %dB attenuation applied to the input per dB exceeding output threshold

        %------------------------------------------------------------------

        % Band filtering properties

        %------------------------------------------------------------------

        bwOct = 1/2; %1/1, 1/2, 1/3, 1/4, 1/5

        filterOrder = 2 %BUTTER=2, GTF=3

        useGTF = false; %If false, revert to butterworth

    end

    %% **********************************************************

    % Read only properties that are not dependent

    %************************************************************

    properties(SetAccess = private)          

        MOCrecord

    end

    %% **********************************************************

    % Constant properties 

    %************************************************************

    properties(Constant = true, Hidden = true)  

        numAudiometricFreqs = 6;

    end

    %% **********************************************************

    % Dependent visable properties - calculated at runtime

    %************************************************************

    properties(Dependent = true, Hidden = false)

        channelBFs %= 250 * 2.^((0:fNmax)'*params.bwOct);        

        numChannels %= numel(channelBFs);        

        aidOPnice %aid output reformatted to be exactly the same dimensions as the input stimulus                               

    end

    %% **********************************************************

    % Dependent invisable properties - calculated at runtime

    %************************************************************

    properties(Dependent = true, Hidden = true)

        TC_dBO_INTERP % Compression threshold in terms of 2nd filter o/p in dB SPL

        TM_dBO_INTERP % MOC threshold in terms of 2nd filter o/p in dB SPL                    

        bwOct_INTERP

        DRNLb_INTERP %=  ( 2e-5 .* 10.^(TCdBO/20)) .^ (1-DRNLc)  ; 

        DRNLc_INTERP               

        mainGain_INTERP %Interp'd and in linear units 

        ARthresholdPa %=  20e-6*10^(ARthreshold_dB/20);% Pa thresh for triggering AR        

        stimulusINTERNAL %input stimulus in correct format for the Aid algo

    end

    %% **********************************************************

    % Protected properties - The user never needs to set

    %************************************************************

    properties(Access = protected)  

        aidOP

        emlc_z

        %--------------------------------------------------------------

        % ENUMERATIONS USED IN THE FRAME PROCESSOR

        %--------------------------------------------------------------

        enumC_ARb  = 0;

        enumC_ARa  = 2;

        enumC_MOCb = 4;

        enumC_MOCa = 6;

        % enumC_BPb1 = 8;

        % enumC_BPa1 = 13;

        % enumC_BPb2 = 18;

        % enumC_BPa2 = 23;

        % enumC_BPb3 = 28;

        % enumC_BPa3 = 33;

        % enumC_BPb4 = 38;

        % enumC_BPa4 = 43;

        enumS_AR  = 0;

        % enumS_MOC1  = 1;

        % enumS_BPin_1_1 = 2;

        % enumS_BPin_2_1 = 6;

        % enumS_BPout_1_1 = 10;

        % enumS_BPout_2_1 = 14;

        %

        % enumS_MOC2 = 18;

        % enumS_BPin_1_2 = 19;

        % enumS_BPin_2_2 = 23;

        % enumS_BPout_1_2 = 27;

        % enumS_BPout_2_2 = 31;

        % ...

    end

    %% **********************************************************

    % methods        _   _               _

    %               | | | |             | |

    % _ __ ___   ___| |_| |__   ___   __| |___

    %| '_ ` _ \ / _ \ __| '_ \ / _ \ / _` / __|

    %| | | | | |  __/ |_| | | | (_) | (_| \__ \

    %|_| |_| |_|\___|\__|_| |_|\___/ \__,_|___/

    %************************************************************

    methods

        %% **********************************************************

        % Constructor

        %************************************************************

        function obj = EssexAid_WrapClass(sr, stimulus)          

            if nargin > 0

                obj.sr = sr;

            end

            if nargin > 1

                obj.stimulusUSER = stimulus;

            else         

                obj.stimulusUSER = obj.pipSequence(obj.sr);

            end

        end

        %% **********************************************************

        % Get method for channelBFs 

        %************************************************************

        function value = get.channelBFs(obj)

            fNmax = 5/obj.bwOct;

            value = 250 * 2.^((0:fNmax)'*obj.bwOct);

        end

        %% **********************************************************

        % Get method for numChannels 

        %************************************************************

        function value = get.numChannels(obj)

            value = numel(obj.channelBFs);

        end

        %% **********************************************************

        % Get method for ARthresholdPa 

        %************************************************************

        function value = get.ARthresholdPa(obj)

            value = 20e-6*10^(obj.ARthreshold_dB/20);% Pa thresh for triggering AR

        end

        %% **********************************************************

        % Get method for TC_dBO_INTERP

        %************************************************************

        function value = get.TC_dBO_INTERP(obj)

            TC_dBO = obj.audiometry_dB - obj.mainGain_dB + obj.TC_dBHL;

            value  = obj.interpPars(TC_dBO, obj.numChannels);

        end

        %% **********************************************************

        % Get method for TM_dBO_INTERP

        %************************************************************

        function value = get.TM_dBO_INTERP(obj)

            TM_dBO = obj.audiometry_dB - obj.mainGain_dB + obj.TM_dBHL;

            value  = obj.interpPars(TM_dBO, obj.numChannels);

        end

        %% **********************************************************

        % Get method for bwOct_INTERP 

        %************************************************************

        function value = get.bwOct_INTERP(obj)

            value = repmat(obj.bwOct, 1, obj.numChannels);

        end

        %% **********************************************************

        % Get method for DRNLb_INTERP 

        %************************************************************

        function value = get.DRNLb_INTERP(obj)

            value = ( 2e-5 .* 10.^(obj.TC_dBO_INTERP/20)) .^ (1-obj.DRNLc_INTERP); 

        end

        %% **********************************************************

        % Get method for DRNLc_INTERP 

        %************************************************************

        function value = get.DRNLc_INTERP(obj)

            value  = obj.interpPars(obj.DRNLc, obj.numChannels);

        end

        %% **********************************************************

        % Get method for mainGain_INTERP 

        %************************************************************

        function value = get.mainGain_INTERP(obj)

            mainGainLin = 10.^(obj.mainGain_dB/20); %lin units

            value  = obj.interpPars(mainGainLin, obj.numChannels);

        end                

        %% ***********************************************************

        % Get method for stimulus

        % -----------------------

        % The hearing aid expects a stereo signal, as the MOC control is

        % linked for left and right channels. It would be more efficient to

        % use a mono version of the aid for simulation in Matlab. However,

        % I always want to use the exact same code for the hardware in the

        % lab and current simulations. This code will make a mono signal

        % stereo if needs be and/or rotate to 2xN array.

        %*************************************************************

        function value = get.stimulusINTERNAL(obj)            

            [nRows, nCols] = size(obj.stimulusUSER);

            % Assume that the stimulus duration is greater than 2 samples.

            % Therefore the number of channels is the min dim.            

            [nChans, I] = min([nRows nCols]);                                                

            if nChans == 2

                if I == 2

                    value = obj.stimulusUSER;

                else

                    value = obj.stimulusUSER';

                end

            elseif nChans == 1 %Just to be explicit

                if I == 2

                    value = [obj.stimulusUSER obj.stimulusUSER];   

                else

                    value = [obj.stimulusUSER; obj.stimulusUSER]';

                end

            end

        end

        %% ***********************************************************

        % Get method for aid output

        % -----------------------

        % This get method is linked to the above internal stimulus method

        % and allows the user to extract the hearing aid output in exactly

        % the same shape and size as the original input stimulus. This is

        % very useful for the speech recognition work and presumably

        % for multithreshold also.

        %*************************************************************

        function value = get.aidOPnice(obj)

            if ~isempty(obj.aidOP)

                [nRows, nCols] = size(obj.stimulusUSER);

                % Assume that the stimulus duration is greater than 2 samples.

                % Therefore the number of channels is the min dim.

                [nChans, I] = min([nRows nCols]);

                %** The aid output will ALWAYS be a 2xN array **

                %The fist job is to remove trailing zeros that may have been

                %introduced by the framing process

                aidOPtruncated = obj.aidOP(:, 1:max([nRows nCols]));

                %The next task is to arrange the op like the ip

                if nChans == 2

                    if I == 1

                        value = aidOPtruncated;

                    else

                        value = aidOPtruncated';

                    end

                elseif nChans == 1 %Just to be explicit

                    if I == 1

                        value = aidOPtruncated(1,:);

                    else

                        value = aidOPtruncated(1,:)';

                    end

                end

            else % ---- of if isempty statement

                value = [];

            end

        end

        %% ***********************************************************

        % *** Set methods ***

        % -----------------------

        % This is a bunch of unexciting error hunting functions. They also

        % flush the aid output if any parameters change. Therefore,

        % processStim will have to be called explicity by the user once

        % again. 

        %*************************************************************

        function obj = set.stimulusUSER(obj,value)

            [nRows, nCols] = size(value);

            % Assume that the stimulus duration is greater than 2 samples.

            % Therefore the number of channels is the min dim.            

            nChans = min([nRows nCols]);            

            assert(nChans<3 && nChans, 'Number of stimulus channels must be 1 or 2')

            obj = obj.flushAidData; %flush any previous hearing aid data if the input stimulus changes

            obj.stimulusUSER = value;

        end                     

        function obj = set.sr(obj,value)

            assert(value>=20e3 && value<=192e3, 'sr must be between 20 and 192 kHz')            

            obj = obj.flushAidData;

            obj.sr = value;

        end

        function obj = set.numSamples(obj,value)

            assert(value>=48 && value<=6912, 'must be between 48 and 6912 samples')            

            obj = obj.flushAidData;

            obj.numSamples = value;

        end

        function obj = set.audiometry_dB(obj,value)

            [nRows,nCols] = size(value);

            assert(nRows==obj.numAudiometricFreqs && nCols==1, 'must be 6x1 column vector') %#ok<MCSUP>

            obj = obj.flushAidData;

            obj.audiometry_dB = value;

        end

        function obj = set.mainGain_dB(obj,value)

            [nRows,nCols] = size(value);

            assert(nRows==obj.numAudiometricFreqs && nCols==1, 'must be 6x1 column vector') %#ok<MCSUP>

            obj = obj.flushAidData;

            obj.mainGain_dB = value;

        end

        function obj = set.TC_dBHL(obj,value)

            [nRows,nCols] = size(value);

            assert(nRows==obj.numAudiometricFreqs && nCols==1, 'must be 6x1 column vector') %#ok<MCSUP>

            obj = obj.flushAidData;

            obj.TC_dBHL = value;

        end

        function obj = set.TM_dBHL(obj,value)

            [nRows,nCols] = size(value);

            assert(nRows==obj.numAudiometricFreqs && nCols==1, 'must be 6x1 column vector') %#ok<MCSUP>

            obj = obj.flushAidData;

            obj.TM_dBHL = value;

        end

        function obj = set.DRNLc(obj,value)

            [nRows,nCols] = size(value);

            assert(nRows==obj.numAudiometricFreqs && nCols==1, 'must be 6x1 column vector') %#ok<MCSUP>

            assert(all(value)>=0 && all(value)<=1, 'all DRNLc values must be between 0 and 1')

            obj = obj.flushAidData;

            obj.DRNLc = value;

        end

        function obj = set.ARtau(obj,value)

            assert(value>=1e-3 && value<=1, 'must be between 1e-3 and 1s')            

            obj = obj.flushAidData;

            obj.ARtau = value;

        end

        function obj = set.ARthreshold_dB(obj,value)

            assert(value>0, 'set AR to a high value to disable it')            

            obj = obj.flushAidData;

            obj.ARthreshold_dB = value;

        end

        function obj = set.MOCtau(obj,value)

            assert(value>=1e-3 && value<=2, 'must be between 1e-3 and 2s')            

            obj = obj.flushAidData;

            obj.MOCtau = value;

        end

        function obj = set.MOCfactor(obj,value)

            assert(value>=0 && value<=1, 'must be between 0 and 1')

            obj = obj.flushAidData;

            obj.MOCfactor = value;

        end

        function obj = set.bwOct(obj,value)

            assert(value==1/1 || value==1/2 || value==1/3 || value==1/4 || value==1/5, 'must be one of 1./(1:5)')

            obj = obj.flushAidData;

            obj.bwOct = value;

        end

        function obj = set.filterOrder(obj,value)

            assert(value>0 && value<5, 'must be one of 1:4')

            obj = obj.flushAidData;

            obj.filterOrder = value;

        end

        function obj = set.useGTF(obj,value)

            obj = obj.flushAidData;

            obj.useGTF = value;

        end                    

        %% **********************************************************

        % flushAidData        

        % This second function is a workaround allowing a set method to

        % change another property value.

        %************************************************************

        function obj = flushAidData(obj) 

            obj.aidOP = [];

            obj.MOCrecord = [];

        end 

        %% **********************************************************

        % OVERLOADED plot method 

        %************************************************************

        function plot(obj)

            clf

            sig2dBSPL = @(sig)20*log10(abs(sig/20e-6)+(1/(2^32)));

            dt = 1/obj.sr;

            tAxis = dt:dt:dt*size(obj.stimulusINTERNAL,1);

            subplot(2,1,1)

            plot(tAxis(1:length(obj.stimulusUSER)), sig2dBSPL(obj.stimulusUSER), 'k')

            if ~isempty(obj.aidOPnice)

                hold on

                plot(tAxis(1:length(obj.stimulusUSER)), sig2dBSPL(obj.aidOPnice), 'r')

            end                                   

            ylim([0 100])

            xlim([0 tAxis(length(obj.stimulusUSER))])

            title('Level response')

            xlabel('Time in seconds')

            ylabel('Level in dB SPL')

            subplot(2,1,2) 

            if ~isempty(obj.MOCrecord)

                imagesc(tAxis, 1:obj.numChannels, flipud(-20*log10(obj.MOCrecord)))

                colorbar                               

            end         

            title('MOC attenuation')

            xlabel('Time in seconds')

            ylabel('Band frequency in Hz')

            numSpacers = 1 + (obj.numChannels-numel(obj.DRNLc)) / (numel(obj.DRNLc)-1);

            set(gca, 'YTick', 1:numSpacers:obj.numChannels);

            set(gca, 'YTickLabel', num2str(flipud([250; 500; 1000; 2000; 4000; 8000])));

        end% ------ OVERLOADED plot method

        %% **********************************************************

        % OVERLOADED soundsc method 

        %************************************************************

        function soundsc(obj)

            soundsc(obj.aidOPnice, obj.sr)

        end                 

        %% **********************************************************

        % processStim

        %************************************************************

        function obj = processStim(obj)

            %--------------------------------------------------------------

            % EMULATION OF THE GUI PARAMETER CONVERSIONS

            %--------------------------------------------------------------

            biggestNumSamples = obj.numSamples; 

            filterStatesL = (zeros(3000,1));

            filterStatesR = filterStatesL;

            filterCoeffs = (zeros(5000,1));

            %filter coefficients

            ARcutOff=1/(2*pi*obj.ARtau);

            [b,a] = butter(1,ARcutOff/(obj.sr/2));

            filterCoeffs(obj.enumC_ARb+1:obj.enumC_ARb+2) = b;

            filterCoeffs(obj.enumC_ARa+1:obj.enumC_ARa+2) = a;

            MOCcutOff=1/(2*pi*obj.MOCtau);

            [bMOC,aMOC] = butter(1,MOCcutOff/(obj.sr/2));

            filterCoeffs(obj.enumC_MOCb+1:obj.enumC_MOCb+2) = bMOC;

            filterCoeffs(obj.enumC_MOCa+1:obj.enumC_MOCa+2) = aMOC;

            for filterCount = 1:obj.numChannels

                %-----------------------------------

                % nonlinear path - filter bws

                %-----------------------------------

                lowerCutOff=obj.channelBFs(filterCount)*2^(-obj.bwOct_INTERP(filterCount)/2);

                upperCutOff=obj.channelBFs(filterCount)*2^( obj.bwOct_INTERP(filterCount)/2);

                if obj.useGTF

                    bwHz = upperCutOff - lowerCutOff;

                    [b_DRNL,a_DRNL] = obj.gammatone(bwHz, obj.channelBFs(filterCount), 1/obj.sr);

                    filterCoeffs(10*(filterCount-1)+9 :10*(filterCount-1)+10) = b_DRNL;

                    filterCoeffs(10*(filterCount-1)+14:10*(filterCount-1)+16) = a_DRNL;

                else                                             

                    [b_DRNL,a_DRNL] = butter(2,[lowerCutOff upperCutOff]/(obj.sr/2));

                    filterCoeffs(10*(filterCount-1)+9 :10*(filterCount-1)+13) = b_DRNL;

                    filterCoeffs(10*(filterCount-1)+14:10*(filterCount-1)+18) = a_DRNL;

                end

            end

            %--------------------------------------------------------------

            % EMULATION OF THE IO CALLBACK THREAD

            %--------------------------------------------------------------

            frameBufferL = buffer(obj.stimulusINTERNAL(:,1), obj.numSamples);

            frameBufferR = buffer(obj.stimulusINTERNAL(:,2), obj.numSamples);

            nFrames = size(frameBufferL,2);

            pad = zeros(1,biggestNumSamples-obj.numSamples);

            ARampL=ones(1,biggestNumSamples);

            ARampR = ARampL;

            MOCcontrol = ones(obj.numChannels, biggestNumSamples);

            peakIPL = zeros(5,1);

            peakOPL = peakIPL;

            rmsIPL  = peakIPL;

            rmsOPL  = peakIPL;

            peakIPR = peakIPL;

            peakOPR = peakIPL;

            rmsIPR  = peakIPL;

            rmsOPR  = peakIPL;

            MOCend = zeros(obj.numChannels,1);

            op = [];

            moc= [];

            for nn = 1:nFrames

                frameBufferPadL = [frameBufferL(:,nn)' pad];

                frameBufferPadR = [frameBufferR(:,nn)' pad];

                [ outBufferL, outBufferR, filterStatesL, filterStatesR,  ARampL, ARampR, MOCend, peakIPL, peakOPL, rmsIPL, rmsOPL, peakIPR, peakOPR, rmsIPR, rmsOPR, MOCcontrol ] =...

                    EssexAidProcessVFrameSwitchable( ...

                    frameBufferPadL,...

                    frameBufferPadR,...

                    filterStatesL,...

                    filterStatesR,...

                    filterCoeffs,...

                    obj.numChannels,...

                    obj.numSamples,...

                    ARampL,...

                    ARampR,...

                    obj.ARthresholdPa,...

                    obj.filterOrder,...

                    obj.DRNLb_INTERP,...

                    obj.DRNLc_INTERP,...

                    obj.TM_dBO_INTERP,...

                    obj.MOCfactor,...

                    peakIPL,...

                    peakOPL,...

                    rmsIPL,...

                    rmsOPL,...

                    peakIPR,...

                    peakOPR,...

                    rmsIPR,...

                    rmsOPR,...

                    MOCend,...

                    MOCcontrol,...

                    obj.mainGain_INTERP,...

                    obj.useGTF);

                outBuffer = ( [outBufferL(:, 1:obj.numSamples); outBufferR(:, 1:obj.numSamples)] );                

                op = [op outBuffer]; %#ok<AGROW>   

                moc= [moc MOCcontrol]; %#ok<AGROW>

            end %End of frame processing emulation loop

            obj.aidOP = op;

            obj.MOCrecord=moc;

        end %End of process stim method          

    end %End of methods block

    %%  *********************************************************

    %      _        _   _                       _   _               _

    %     | |      | | (_)                     | | | |             | |

    %  ___| |_ __ _| |_ _  ___   _ __ ___   ___| |_| |__   ___   __| |___

    % / __| __/ _` | __| |/ __| | '_ ` _ \ / _ \ __| '_ \ / _ \ / _` / __|

    % \__ \ || (_| | |_| | (__  | | | | | |  __/ |_| | | | (_) | (_| \__ \

    % |___/\__\__,_|\__|_|\___| |_| |_| |_|\___|\__|_| |_|\___/ \__,_|___/

    %************************************************************

    methods(Static)        

        %% ********************************************************

        % pipOut - sequence of tone pips at various levels

        %**********************************************************        

        function pipOut = pipSequence(sampleRate, freq, dBlevs, pulseDur, silDur)

            if nargin < 5

                silDur = 0.3;

            end

            if nargin < 4

                pulseDur = 0.1;

            end

            if nargin < 3

                dBlevs = 20:20:100;

            end

            if nargin < 2

                freq = 500;

            end

            if nargin < 1

                sampleRate = 48e3;

            end

            dt = 1/sampleRate;

            tAxis = dt:dt:pulseDur;

            sPulse = sin(2*pi*freq*tAxis);

            sPulse = sPulse./sqrt(mean(sPulse.^2));

            rms2dBspl = @(dBspl)20e-6*10^(dBspl/20); %sneaky short-hand function by (ab)using function handles

            zPad = zeros(1,ceil(sampleRate*silDur));

            pipOut = [];

            for nn = 1:numel(dBlevs)                

                pipOut = [ pipOut sPulse*rms2dBspl(dBlevs(nn))  zPad]; %#ok<AGROW>

            end

        end% ------ OF pipSequence

        %% ********************************************************

        % interpPars - Linear interpolation of given parameter to mimic GUI

        % fitting functionality.

        %**********************************************************        

        function fullArray = interpPars(shortArray, numBands)

            nGUIbands = numel(shortArray);

            if numBands == nGUIbands

                fullArray = shortArray;

            else

                numSpacers = (numBands-nGUIbands) / (nGUIbands-1);

                fullArray = shortArray(1);

                for nn = 2:nGUIbands

                    fullArray = [fullArray,...

                        repmat(mean([shortArray(nn) shortArray(nn-1)]),1,numSpacers),...

                        shortArray(nn)]; %#ok<AGROW>

                end                    

            end                        

        end% ----- OF interpPars

        %% ********************************************************

        % gammatone - get filter coefficients

        %********************************************************** 

        function [b,a] = gammatone(bw, cf, dt)

            phi = 2 * pi * bw * dt;

            theta = 2 * pi * cf * dt;

            cos_theta = cos(theta);

            sin_theta = sin(theta);

            alpha = -exp(-phi) * cos_theta;

            b0 = 1.0;

            b1 = 2 * alpha;

            b2 = exp(-2 * phi);

            z1 = (1 + alpha * cos_theta) - (alpha * sin_theta) * 1i;

            z2 = (1 + b1 * cos_theta) - (b1 * sin_theta) * 1i;

            z3 = (b2 * cos(2 * theta)) - (b2 * sin(2 * theta)) * 1i;

            tf = (z2 + z3) / z1;

            a0 = abs(tf);

            a1 = alpha * a0;

            a = [b0, b1, b2];

            b = [a0, a1];

        end% ------ OF gammatone

    end% ------ OF static methods

end %End of classdef

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%   This program is free software; you can redistribute it and/or modify

%   it under the terms of the GNU General Public License as published by

%   the Free Software Foundation; either version 2 of the License, or

%   (at your option) any later version.

%

%   This program is distributed in the hope that it will be useful,

%   but WITHOUT ANY WARRANTY; without even the implied warranty of

%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

%   GNU General Public License for more details.

%

%   You can obtain a copy of the GNU General Public License from

%   http://www.gnu.org/copyleft/gpl.html or by writing to

%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

classdef cHMM

    %HMMCLASS Handles all of the HTK related gubbins

    %   Please see the documentation located in a separate file for further

    %   information.     

    %%  *********************************************************

    %  properties                      _   _

    %                                 | | (_)

    %  _ __  _ __ ___  _ __   ___ _ __| |_ _  ___  ___

    % | '_ \| '__/ _ \| '_ \ / _ \ '__| __| |/ _ \/ __|

    % | |_) | | | (_) | |_) |  __/ |  | |_| |  __/\__ \

    % | .__/|_|  \___/| .__/ \___|_|   \__|_|\___||___/

    % | |             | |

    % |_|             |_|

    %************************************************************

    %% **********************************************************

    % Public properties - can be set by user

    %************************************************************

    properties(Access = public)

        hmmFolder

        paramType               = 'USER_D_A'; %DELTAS and ACCELERATIONS

        numCoeff                = 27; %9*3 THIS IS FOR PROBABILITY MODEL (not high spont+low spont which would be 18*3=54)

        HERestDataPath      = fullfile(pwd, 'def', 'HERest_digit');

        binPath             = fullfile(pwd, 'def', 'bin');

        configFile          = fullfile(pwd, 'def', 'config_STANDARD');

        trainWordListFile   = fullfile(pwd, 'def', 'Grammar_digit', 'words3');

        testWordListFile    = fullfile(pwd, 'def', 'Grammar_digit', 'wordsNoSevenZero');

        wordNetFile         = fullfile(pwd, 'def', 'Grammar_digit', 'wdnetNoSP.slf');

        dictFile            = fullfile(pwd, 'def', 'Grammar_digit', 'noSevenZeroDict');

    end

    %% **********************************************************

    % Dependent - never set by user. Only calculated when needed

    %************************************************************

    properties(Dependent = true)

        protoFile              % = fullfile(pwd, 'def', 'proto_RobANonly_9'); %probability only

    end

    %%  *********************************************************

    % methods        _   _               _

    %               | | | |             | |

    % _ __ ___   ___| |_| |__   ___   __| |___

    %| '_ ` _ \ / _ \ __| '_ \ / _ \ / _` / __|

    %| | | | | |  __/ |_| | | | (_) | (_| \__ \

    %|_| |_| |_|\___|\__|_| |_|\___/ \__,_|___/

    %************************************************************

    methods

        %% **********************************************************

        % Constructor

        %************************************************************

        function obj = cHMM(hmmFolder)

            if nargin > 0

                obj.hmmFolder = hmmFolder;

            end

        end % ------ OF CONSTRUCTOR

        %% **********************************************************

        % genProto - generate task specific prototype

        %************************************************************

        function genProto(obj)

            % models_1mixsil.exe - takes input (hmmdef) and copies it making a

            % basis of one, two three etc. etc.

            obj.models_1mixsilMat(fullfile(obj.hmmFolder,'hmm0','hmmdef'), fullfile(obj.hmmFolder,'hmm0','models'));            

        end % ------ OF GENPROTO

        %% **********************************************************

        % istrained

        %************************************************************

        function boolans = istrained(obj)

            boolans = numel(dir(fullfile(obj.hmmFolder,'hmm36','models')));

        end

        %% **********************************************************

        % Train

        %************************************************************

        function train(obj, trainFeatureFolder)

            % Most of the following code block can be replaced by doing a

            % find and replace across the  code below. I didn't

            % want to mess with the working code too much, so I just

            % copied the object properties needed into the variable names

            % expected.

            ED_CMDFILE1 = fullfile(obj.HERestDataPath, 'sil1.hed');

            ED_CMDFILE2 = fullfile(obj.HERestDataPath, 'mix2_16.hed');

            ED_CMDFILE3 = fullfile(obj.HERestDataPath, 'mix3_16.hed');

            ED_CMDFILE4 = fullfile(obj.HERestDataPath, 'mix5_16.hed');

            ED_CMDFILE5 = fullfile(obj.HERestDataPath, 'mix7_16.hed');

            NUM_COEF        = obj.numCoeff;

            PAR_TYPE        = obj.paramType;

            LIST_FILE       = fullfile(obj.hmmFolder, 'tmp.list');

            word_list       = obj.trainWordListFile;%fullfile(obj.grammarPath, 'words');

            word_listSP     = word_list; % for use in hmm4 onwards - UGLY HACK NOW SP ABANDONED

            proto           = obj.protoFile;%Does not exist on disk just yet probably - see a few lines down

            config          = obj.configFile;

            train_list      = fullfile(trainFeatureFolder, 'list.scp');

            labels          = fullfile(trainFeatureFolder, 'labels.mlf');

            labelssp        = labels; % for use in hmm4 onwards - UGLY HACK NOW SP ABANDONED

            hmm_dir         = obj.hmmFolder;

            FEAT_ROOT       = trainFeatureFolder;

            % Now for the actual HMM training code

            mkdir(hmm_dir)

            for I = 0:36

                h = fullfile(hmm_dir,['hmm' num2str(I)]);

                mkdir(h);

            end

            obj.makeProtoHmm(proto, obj.paramType, obj.numCoeff, 18);

            fid = fopen(train_list,'r');

            disp(train_list)

            S = textscan(fid,'%s','Delimiter','\n');

            fclose(fid);

            fid = fopen(LIST_FILE,'w');

            for I = 1:size(S{1},1)

                str  = fullfile(FEAT_ROOT,S{1}{I});

                fprintf(fid,'%s\n',str);

            end;

            fclose(fid);

            % HCompV just gets the vfloor stuff out so we can begin approximating

            cmd = ['"HCompV" -T 2 -D -C "' config '" -o hmmdef -f 0.01 -m -S "' LIST_FILE '" -M "' hmm_dir filesep 'hmm0" "' proto '"'];

            system(cmd);                        

            %cmd = [BINDIR filesep 'macro' binExt ' ' num2str(NUM_COEF) ' ' PAR_TYPE ' ' '"' hmm_dir filesep 'hmm0' filesep 'vFloors' '" "' hmm_dir filesep 'hmm0' filesep 'macros' '"'];

            %system(cmd);            

            obj.macroMat(NUM_COEF,PAR_TYPE, fullfile(hmm_dir, 'hmm0', 'vFloors'), fullfile(hmm_dir, 'hmm0', 'macros'));

            %MAKE THE INITIAL MODEL PROTOTYPE

            genProto(obj);

            disp('Seed HMM successfully Produced.....');

            %Training

            for I = 1:3

                disp(I)

                j = I-1;

                cmd = ['HERest -D -C ' config ' -I ' labels ' -t 250.0 150.0 1000.0 -S ' LIST_FILE  ' -H ' hmm_dir filesep 'hmm' num2str(j)  filesep 'macros -H ' hmm_dir filesep 'hmm' num2str(j) filesep 'models -M ' hmm_dir filesep 'hmm' num2str(I) ' ' word_list];

                %                 disp(cmd)

                system(cmd);

            end

            disp('3 iterations complete');

            rmdir ([hmm_dir filesep 'hmm4'],'s')

            copyfile ([hmm_dir filesep 'hmm3'], [hmm_dir filesep 'hmm4'])                                                

            % The following command takes state 3 from the silence model

            % and appends it to the end of the model as state 2 of the

            % short pause model.

            % Original:

            % cmd = [BINDIR filesep 'spmodel_gen' binExt ' ' hmm_dir filesep 'hmm3' filesep 'models ' hmm_dir filesep 'hmm4' filesep 'models'];

            % system(cmd);

            % New:

            obj.spmodel_genMat(fullfile(hmm_dir,'hmm3','models'), fullfile(hmm_dir,'hmm4','models'));                                     

            cmd = ['HHEd  -T 2 -H ' hmm_dir filesep 'hmm4' filesep 'macros -H ' hmm_dir filesep 'hmm4' filesep 'models -M ' hmm_dir filesep 'hmm5 ' ED_CMDFILE1 ' ' word_listSP ];

            system(cmd);

            disp ('SP model fixed')

            % after the spmodel_gen command - the word_list is changed to

            % word_listSP. The sp model is just ignored currently

            for I = 6:8

                disp(I)

                j = I-1;

                cmd = ['HERest -C ' config ' -I ' labelssp ' -S ' LIST_FILE  ' -H ' hmm_dir filesep 'hmm' num2str(j)  filesep 'macros -H ' hmm_dir filesep 'hmm' num2str(j) filesep 'models -M ' hmm_dir filesep 'hmm' num2str(I) ' ' word_listSP];

                system(cmd);

            end

            disp('6 iterations complete');

            cmd = ['HHEd  -T 2 -H ' hmm_dir filesep 'hmm8' filesep 'macros -H ' hmm_dir filesep 'hmm8' filesep 'models -M ' hmm_dir filesep 'hmm9 ' ED_CMDFILE2 ' ' word_listSP ];

            system(cmd);

            disp ('2 gaussians per mixture')

            for I = 10:12

                disp(I)

                j = I-1;

                cmd = ['HERest -C ' config ' -I ' labelssp ' -S ' LIST_FILE  ' -H ' hmm_dir filesep 'hmm' num2str(j)  filesep 'macros -H ' hmm_dir filesep 'hmm' num2str(j) filesep 'models -M ' hmm_dir filesep 'hmm' num2str(I) ' ' word_listSP];

                system(cmd);

            end

            disp ('9 iterations completed')

            cmd = ['HHEd  -T 2 -H ' hmm_dir filesep 'hmm12'  filesep 'macros -H ' hmm_dir filesep 'hmm12' filesep 'models -M ' hmm_dir filesep 'hmm13 ' ED_CMDFILE3 ' ' word_listSP ];

            system(cmd);

            disp ('3 gaussians per mixture')

            for I = 14:20

                disp(I)

                j = I-1;

                cmd = ['HERest -C ' config ' -I ' labelssp ' -S ' LIST_FILE  ' -H ' hmm_dir filesep 'hmm' num2str(j)  filesep 'macros -H ' hmm_dir filesep 'hmm' num2str(j) filesep 'models -M ' hmm_dir filesep 'hmm' num2str(I) ' ' word_listSP];

                system(cmd);

            end

            disp ('16 iterations completed')

            cmd = ['HHEd  -T 2 -H ' hmm_dir filesep 'hmm20'  filesep 'macros -H ' hmm_dir filesep 'hmm20' filesep 'models -M ' hmm_dir filesep 'hmm21 ' ED_CMDFILE4 ' ' word_listSP ];

            system(cmd);

            disp ('5 gaussians per mixture')

            for I = 22:28

                disp(I)

                j = I-1;

                cmd = ['HERest -C ' config ' -I ' labelssp ' -S ' LIST_FILE  ' -H ' hmm_dir filesep 'hmm' num2str(j)  filesep 'macros -H ' hmm_dir filesep 'hmm' num2str(j) filesep 'models -M ' hmm_dir filesep 'hmm' num2str(I) ' ' word_listSP];

                system(cmd);

            end

            disp ('23 iterations completed')

            cmd = ['HHEd  -T 2 -H ' hmm_dir filesep 'hmm28'  filesep 'macros -H ' hmm_dir filesep 'hmm28' filesep 'models -M ' hmm_dir filesep 'hmm29 ' ED_CMDFILE5 ' ' word_listSP ];

            system(cmd);

            disp ('7 gaussians per mixture')

            for I = 30:36

                disp(I)

                j = I-1;

                cmd = ['HERest -C ' config ' -I ' labelssp ' -S ' LIST_FILE  ' -H ' hmm_dir filesep 'hmm' num2str(j)  filesep 'macros -H ' hmm_dir filesep 'hmm' num2str(j) filesep 'models -M ' hmm_dir filesep 'hmm' num2str(I) ' ' word_listSP];

                system(cmd);

            end

            disp ('30 iterations completed')

        end % ------ OF TRAIN

        %% **********************************************************

        % Test

        %************************************************************

        function test(obj, testFeatureFolder)

            flags = '-p 0.0 -s 0.0';

            test_word_list = obj.testWordListFile;

            net = obj.wordNetFile;

            dict = obj.dictFile;

            LIST_FILE       = fullfile(testFeatureFolder, 'tmp.list');

            config          = obj.configFile;

            disp ('Now testing with 7 mixture HMMs')

            mod_file = fullfile(obj.hmmFolder, 'hmm36', 'models');

            mac_file = fullfile(obj.hmmFolder, 'hmm36', 'macros');

            RESULTS_DIR = testFeatureFolder;

            TEST_FEAT_ROOT  = testFeatureFolder;

            N1list =  fullfile(testFeatureFolder, 'list.scp');

            ftest = fopen(N1list,'r');

            S = textscan(ftest,'%s','Delimiter','\n');

            fclose(ftest);

            flist = fopen(LIST_FILE,'w');

            for I = 1:size(S{1},1)

                str  = fullfile(TEST_FEAT_ROOT,S{1}{I});

                fprintf(flist,'%s\n',str);

            end;

            fclose(flist);

            cmd = ['HVite -D -H ' mac_file ' -H ' mod_file ' -S ' LIST_FILE ' -C ' config ' -w ' net ' -l ''*'' -i ' RESULTS_DIR filesep 'result.mlf ' flags ' ' dict ' ' test_word_list];

            system(cmd);

        end % ------ OF TEST

        %% **********************************************************

        % Get methods determining feature vector related gubbins

        %************************************************************

        function value = get.protoFile(obj)

            value = fullfile(obj.hmmFolder, 'proto_AutoGen');

        end

    end % ------ OF METHODS

    %%  *********************************************************

    %      _        _   _                       _   _               _

    %     | |      | | (_)                     | | | |             | |

    %  ___| |_ __ _| |_ _  ___   _ __ ___   ___| |_| |__   ___   __| |___

    % / __| __/ _` | __| |/ __| | '_ ` _ \ / _ \ __| '_ \ / _ \ / _` / __|

    % \__ \ || (_| | |_| | (__  | | | | | |  __/ |_| | | | (_) | (_| \__ \