function Benchmark_NN(input, output, hiddenLayerSize, step)

    for i=1:size(input,2)
        input(:,i) = mat2gray(input(:,i));
    end

    function [accuracyTrain, accuracyVal, accuracyTest, accuracyAll] = BackPropagationNeuralNetwork()
        x = input';
        t = output';

        % Choose a Training Function
        % For a list of all training functions type: help nntrain
        % 'trainlm' is usually fastest. (Levenberg-Marquardt)
        % 'trainbr' takes longer but may be better for challenging problems.
        % 'trainscg' uses less memory. Suitable in low memory situations.
        trainFcn = 'trainscg';

        % Create a Pattern Recognition Network
        net = patternnet(hiddenLayerSize, trainFcn);
        net.trainParam.showWindow = 0;

        % Setup Division of Data for Training, Validation, Testing
        net.divideParam.trainRatio = 70/100;
        net.divideParam.valRatio = 15/100;
        net.divideParam.testRatio = 15/100;

        % Train the Network
        [net, tr] = train(net,x,t);

        % Training Confusion Plot Variables
        yTrn = net(x(:,tr.trainInd));
        tTrn = t(:,tr.trainInd);
        [c,~] = confusion(tTrn,yTrn);
        accuracyTrain = 100*(1-c);      

        % Validation Confusion Plot Variables
        yVal = net(x(:,tr.valInd));
        tVal = t(:,tr.valInd);
        [c,~] = confusion(tVal,yVal);
        accuracyVal = 100*(1-c);     

        % Test Confusion Plot Variables
        yTst = net(x(:,tr.testInd));
        tTst = t(:,tr.testInd);
        [c,~] = confusion(tTst,yTst);
        accuracyTest = 100*(1-c);    

        % Overall Confusion Plot Variables
        yAll = net(x);
        tAll = t;
        [c,~] = confusion(tAll,yAll);
        accuracyAll = 100*(1-c);

        % Plot Confusion
        %plotconfusion(tTrn, yTrn, 'Training', tVal, yVal, 'Validation', tTst, yTst, 'Test', tAll, yAll, 'Overall')

        % View the Network
        %view(net)

        % Plots
        % Uncomment these lines to enable various plots.
        %figure, plotperform(tr)
        %figure, plottrainstate(tr)
        %figure, ploterrhist(e)
        %figure, plotconfusion(t,y)
        %figure, plotroc(t,y)
    end

    acc1 = zeros(step,1);
    acc2 = zeros(step,1);    
    acc3 = zeros(step,1);    
    acc4 = zeros(step,1);    

    for i=1 : step
        %disp(i);
        [accuracyTrain, accuracyVal, accuracyTest, accuracyAll] = BackPropagationNeuralNetwork();
        acc1(i) = accuracyTrain;
        acc2(i) = accuracyVal;
        acc3(i) = accuracyTest;
        acc4(i) = accuracyAll;
    end
    
    fprintf('%2.2f±%2.2f\t%2.2f±%2.2f\t%2.2f±%2.2f\t%2.2f±%2.2f\n', ...
        mean(acc1), std(acc1), mean(acc2), std(acc2), mean(acc3), std(acc3), ...
        mean(acc4), std(acc4));
end
