Code covered by the BSD License  

Highlights from
Nth_Oct_Hand_Arm_&_AC_Filter_Tool_Box

Nth_Oct_Hand_Arm_&_AC_Filter_Tool_Box

by

 

22 Dec 2008 (Updated )

Features Nth octave band, Hand Arm, and A and C weighting filters

[ndraw, count, count2, error]=rand_int(a, b, n, sorted, no_duplicates, make_plot)
function [ndraw, count, count2, error]=rand_int(a, b, n, sorted, no_duplicates, make_plot)
% % rand_int: Quckly generates n random integers from a to b integer 
% % 
% % Syntax;
% % 
% % [ndraw, count, count2, error]=rand_int(a, b, n, sorted, no_duplicates, make_plot);
% % 
% % **********************************************************************
% % 
% % Description
% % 
% % This program quickly and randomly selects n integers from a to b.
% % The uniform distribution is used so that all integers are more equally
% % probable.  No claims are made about the actual probability
% % distributions.
% %
% % **********************************************************************
% % 
% % Input Arguments      
% %
% % a and b are the lower and upper values of the range of integers.
% % Either a or b can be the higher and lower values of the range.
% % The defaults are a=1 and b=10.
% % 
% % n is the number of integers to select and can be a constant
% % or a two element vector.
% % The defaults are a=1 and b=10.
% % 
% % If n is a single number, then ndraw is a column vector.
% % The defaults are n=5.
% % 
% % If n is a two element vector then n specifies the size of the output
% % matrix.  Example: n=[mm nn]; size(ndraw)=[mm nn];
% %
% % The user specifies whether the output is sorted 1 or randomly ordered 0.
% % The default is not sorted.  If sorted the output is in ascending order.
% %
% % The user specifies to remove duplicate integers 1 or allow duplicate
% % integers 0.  The default is to remove duplicates 1.  If duplicates are
% % removed, then each integer can occur only once.
% % 
% % If make_plot equals 1, then a plot of the random integers is made. 
% % If make_plot does not equal 1, then no plots are made.
% % The default is to not make plots.
% % The x-axis is the row or column indices whichever has more indices.  
% % The y-axis is the value of the integers.
% %  
% % This option was created to visually check the uniformity of the output
% % distribution. 
% %  
% % **********************************************************************
% % 
% % Output Arguments      
% % 
% % ndraw   is the matrix of random integers.
% % 
% %         If n is a single number then ndraw is a column vector.
% % 
% %         if n is a two element vector then ndraw is a matrix with
% %         size equal to n.  Example n=[mm nn]; size(ndraw)=[mm nn]
% % 
% % count   is the number of iterations of selecting integers.
% % 
% % count2  is the number of iterations of deleting extra random numbers
% % 
% % error   1 there is at least 1 error
% %         0 there are no errors
% % 
% % **********************************************************************
% % 
% % Comments on Speed      
% %
% % This program quickly selects n integers even when
% % b-a is very large (See the Large Integer Range Examples).
% %
% % When duplicates are allowed the program just selects the required
% % number of integers.  It is very simple case.
% %
% % When duplicates are removed, three regimes are used to optimize speed.
% %
% % The code runs the fastest in regime 2 pure rejection code contributed
% % by John D'Errico (woodchips@rochester.rr.com).  The pure rejection 
% % code is forced to quit after 100 unsuccessful iterations, then A catch
% % is used to select the required number of integers.
% %
% % The code runs the slowest when min((b-a+1)-n, n)/(b-a) > 0.43
% % and utilizes a contributed alternative approach in this range
% % from Jos x@y.z.
% %
% % For all other cases, the program selects extra integers, then keeps the
% % unique integers and randomly deselects extra integers.
% % There is a selection while loop and a deselection while loop.
% %
% % See example for syntax and definitions
% % of the input and output arguments
% %
% % **********************************************************************
%
% Example='';       % Large integer Range Example with very small n
%
% a=1;              % a is the lowest integer that can be selected.
%                   % default value is 1
%
% b=10^15;          % b is the highest integer that can be selected.
%                   % default value is 10
%
% n=5;              % n is the number of integers to select.
%                   % default value is 5
%                   %
%                   % n can be a constant or a two element vector
%                   %
%                   % if n is a constant then the output is a column
%                   % vector with n rows
%                   %
%                   % if n is a two element matrix [mm nn]
%                   % then the output is a matrix of size [mm nn]
%
% sorted=1;         % 1 output is sorted (ascending order)
%                   % otherwise not sorted
%                   % default value is 0 (not sorted)
%
% no_duplicates=1;  % 1 there are no duplicate integers in the output
%                   % otherwise allow duplicate integers
%                   % default value is 1 (duplicates removed)
% 
% make_plot=1;      % 1 make plots of the integers
%                   % otherwise do not make any plots
% 
% [ndraw, count]=rand_int(a, b, n, sorted, no_duplicates, make_plot);
% 
% % **********************************************************************
% 
% Example='';       % Very Large integer Range Example with small n
% 
% a=-10^15;
% b=10^15;
% n=1000;
% sorted=1;         %(ascending order)
% no_duplicates=1;
% [ndraw, count, count2]=rand_int(a,b,n,sorted, no_duplicates, make_plot);
% 
% % **********************************************************************
% 
% Example='';       % Large integer Range Example with large n
% 
% a=1;
% b=10^7;
% n=[2.5*10^6 2];  % note that n is a two element matrix and very large
% sorted=0;
% no_duplicates=0;
% [ndraw, count, count2]=rand_int(a,b,n, sorted, no_duplicates, make_plot);
% 
% Examples of the 3 regimes
% 
% Regime 1: Build entire matrix then remove extras
% [ndraw, count, count2, error]=rand_int(-10000, 10000, [ 1000 10 ], 0, 0, 1);
% 
% Regime 2: Pure rejection code (100 percent purge regime).
% [ndraw, count, count2, error]=rand_int(-10000, 10000, [ 10 1 ], 0, 0, 1);
% 
% Regime 3: Iteratively select integers then remove extras if necessary
% [ndraw, count, count2, error]=rand_int(-10000, 10000, [ 100 10 ], 0, 0, 1);
% 
% % **********************************************************************
% 
% Example=''; % Exceeding count limit.  The Pure Rejection Regime Failure
%
% % If the following code is run enough times the pure rejection regime
% % will fail and the catch will select the additional integers.
% 
% a=1;
% b=10000;
% n=250;  % note that n is a two element matrix and very large
% sorted=1;
% no_duplicates=1;
% [ndraw, count, count2]=rand_int(a,b,n, sorted, no_duplicates);
% 
% % The pure rejection regime is forced into failure if it fails to select 
% % enough integers after 100 iterations.  A catch is put into place 
% % to select  all the necessary integers.
% % 
% % 
% % **********************************************************************
% %
% %
% %
% % This program was written by Edward L. Zechmann
% %     date  23 January 2008
% %
% % modified  23 January    2008    Updated comments.
% %                                 Increased speed for a large range of
% %                                 integers.
% %
% % modified  24 January    2008    Modified the code to be faster
% %                                 Added contributed alternative approach
% %                                 from Jos x@y.z.
% %                                 Created a second regime to quicken
% %                                 integer selection when number of selected
% %                                 integers is nearly one-half of the range
% %                                 Added some large integer range examples
% %
% % modified  25 January    2008    Added sorted option for ordering the 
% %                                 output array.  Now output can be sorted 
% %                                 or random.
% %
% %                                 Added an option for n to be a two element
% %                                 matrix stipulating the output matrix size
% %                                 so n=[mm nn];
% %                                 mm is number of rows
% %                                 nn is number of columns
% %
% % modified  26 January    2008    Updated Comments
% %
% %                                 Fixed a bug so that output is within the
% %                                 selected range.
% %
% %                                 Added contributed code from
% %                                 John D'Errico (woodchips@rochester.rr.com)
% %                                 Created a third regime to quickly select
% %                                 a very small number of integers.
% %
% %                                 Added warnings and error output when input
% %                                 does not make sense or is out of range
% %
% % modified  30 January    2008    Rearranged the three regimes so that
% %                                 third regime can catch failures from the
% %                                 100% purge regime when it exceeds
% %                                 the max_count limit.
% %
% % modified  31 January    2008    Fixed bug.  The ouptput is reshaped
% %                                 to the specified size or the best size
% %                                 for ndraw
% %
% %                                 The comments were revised and updated
% %
% % modified  20 February   2008    Added option to make a plot of the integers
% %                                 Added additional examples.
% %                                 Removed pure rejection regime warnings
% %                                 since the catch is reliable.  
% %
% % modified   4 March      2009    Updated Comments.
% %
% %
% % **********************************************************************
% %
% % Feel free to modify this code.
% % 
% % See Also: MYRANDINT
% % 

% set the error flag to no error
% flag = 0 means there is typical output
% flag = 1 mean the output must be truncated because of bad input
flag=0;

if nargin < 1 || (isempty(a) || ~isnumeric(a))
    a=1;
end

if nargin < 2 || (isempty(b) || ~isnumeric(b))
    b=10;
end

if nargin < 3 || (isempty(n) || ~isnumeric(n))
    n=5;
end

if nargin < 4  || (isempty(sorted) || ~isnumeric(sorted))
    sorted=0;
end

if nargin < 5  || (isempty(no_duplicates) || ~isnumeric(no_duplicates))
    no_duplicates=1;
end

if nargin < 6  || (isempty(make_plot) || ~isnumeric(make_plot))
    make_plot=0;
end

if ~isequal(sorted, 1)
    sorted=0;
end

if ~isequal(no_duplicates, 1)
    no_duplicates=0;
end



% Initialize output
ndraw=[];
count=0;
count2=0;
error=0;
mm=0;
nn=0;

if isequal(flag, 1)

else

    % set the error flag to no error
    % flag = 0 means there is typical output
    % flag = 1 mean the output must be truncated because of bad input
    flag=0;
    % set the other flags to null
    flag2=0;
    flag3=0;

    % make sure certain inputs are constant integers
    %
    % n can be a 2 element vector
    % round a up to nearest integer
    a=ceil(a(1));

    % round b down to nearest integer
    b=floor(b(1));

    % make sure n has only 2 dimensions
    n=round(n(:, :, 1));
    sorted=round(sorted(1));
    no_duplicates=round(no_duplicates(1));



    % If n is a constant the output is a column vector
    %
    % if n is a two element vector n=[mm nn];
    % Then the output matrix has size [mm nn];

    if length(n) > 1
        mm=n(1);
        nn=n(2);
        n=prod(n);
    else
        mm=n;
        nn=1;
    end

    % Make sure b > a
    % swap values if necessary
    if a > b
        c=a;
        a=b;
        b=c;
    end

    % Make sure b-a >= 1 and that n >= 1.
    % Make sure there are at least n integers from b to a.
    % There are b-a+1 integers from b to a inclusive.
    % if conditions are false return reasonable output.
    if logical((b-a+1) < 1)
        flag=1;
        ndraw=[];
        warning('Input error!  Empty range of integers.  Returning empty output array');
    end

    % If specified to select less than 1 integer return an empty matrix.
    if logical(n < 1)
        flag=1;
        ndraw=[];
        warning('Input error!  Select less than 1 integer.  Returning empty output array');
    end

    % If specified to select more integers than available then return all
    % availabe integers and print a warning.
    if logical(n > (b-a+1)) && isequal(no_duplicates, 1)
        flag=1;
        ndraw=(a:b)';
        warning('Input error!  Selecting more integers than in the range and choosing no duplicates.  Returning a:b as the output array');
    end

    % if flag is set to 1 then set the error
    % and truncate the output to reasonable output.
    if isequal(flag, 1)
        error=1;
        n2=min([length(ndraw), n]);
        if n2 >= 1
            ndraw=ndraw(1:n2, 1);
        else
            ndraw=[];
        end
    end

    n3=n;
    if n >= (b-a+1)/2
        %reverse algorithm and deselect data points
        flag2=1;
        n=(b-a+1)-n;
    end

    % define the threshold limits for the different regimes
    thresh_regime1=0.43;    % generates entire matrix

    % The regime 2 threshold depends on whether selecting almost all the
    % data points or very few data points
    % not sure what teh optimum threshold is
    % this does not crash
    thresh_regime2l=min([0.07, 7/n3]);  % The probability of a duplicate
    % grows with range and number
    % of selected integers
    thresh_regime2h=min([0.05, 5/((b-a+1)-n3)]);  %

    if min((b-a+1)-n3, n3)/(b-a) > thresh_regime1
        % Regime 1
        % n is nearly equal to half the number of data points
        flag3=1;
    elseif logical(n3 <= (b-a+1)*thresh_regime2l) || logical(n3 >= (b-a+1)*(1-thresh_regime2h))
        % Regime 2
        % n is a small fraction of the range
        flag3=2;
    else
        % Regime 3 catch any failures from the first two regimes
        % Regime 2 will fail. regime 1 should b never fail.
        % Select random integers then remove the extra integers
        flag3=3;
    end

    % maximum number of times the selection while loop will run
    max_count=max(ceil(n3/(n3-(b-a))), 100);

    % nunique is the number of unique integers selected
    nunique=0;

    % count is the number of iterations of the selection while loop
    count=0;

    % count2 is the number of iterations of the deselection while loop
    count2=0;

    num_points=b-a+1;

    if isequal(flag, 1)

    else
        if isequal(no_duplicates, 1)

            % For the case of removing all duplicates
            fail=0;

            if  isequal(flag3, 1)
                % In this regime the entire matrix of values if filled, 
                % then the extras are discarded.
                % This is a faster approach for the case
                % min((b-a+1)-n,n)/(b-a) > thresh_regime1
                % Code was contributed by Jos x@y.z
                % code was modified by sorting the output
                % when sorting is specified
                r=randperm(num_points);
                if isequal(sorted, 1)
                    ndraw=sort(a-1+r(1:n))';
                else
                    ndraw=a-1+r(1:n)';
                end

            elseif isequal(flag3, 2)
                % This is the 100 percent purge regime.
                % There is a possibility of returning an empty array.
                %
                % The program runs faster by rejecting all integers rather than
                % sorting through the input and keeping the good integers
                %
                % This is faster approach for the case
                % (b-a+1)/n > 1-thresh_regime2 || n/(b-a) < thresh_regime2
                %
                % The judicious cuttoff of ten percent minimizes the chances of
                % an error and returning an empty array.

                flag5 = 1;
                count=0;
                while flag5 && logical(count < max_count)
                    count=count+1;
                    ndraw = unique(round((a - 1/2) + (b-a+1)*rand(n,1)))';
                    flag5 = (length(ndraw)~=n);
                end

            end

            % check to find out if enough integers were selected
            if logical(length(ndraw)~=n) && ~isequal(flag3, 3)
                ndraw=[];
                error=1;
                fail=1;
                % This warning is caught by the if statement below
                % warning('Failed to select enough integers');
            end

            % if not enough integers were selected
            if (~isequal(flag3, 1) && ~isequal(flag3, 2)) || fail

                ndraw=[];
                count=0;

                while (logical(nunique < n) || logical(count < 1)) && logical(count < max_count)

                    count=count+1;

                    % To make all combinations more equally likely, stretch the random
                    % distribution to a-0.5 and b+0.5.
                    % Round the numbers to the nearest integers.
                    % Make sure the selected integers are unique.
                    % Only one selection per integer.
                    % select an extra thirty percent of points because there will
                    % duplicates that are discarded.
                    num_draws=max(ceil((1.1+n/num_points)*(n-nunique)));
                    ndraw1=round((a+b)/2+((b+0.5)-(a-0.5))*(-0.5+rand(num_draws, 1)));

                    % make sure the numbers fit in the proper range
                    ndraw1(ndraw1 > b)=b;
                    ndraw1(ndraw1 < a)=a;

                    ndraw=unique([ndraw' ndraw1'])';

                    % Count the number of unique integers selected
                    nunique=length(ndraw);
                    count2=0;

                    if nunique > n
                        ndraw_del=[];
                        num_draws=ceil(1.3*(nunique-n));
                        aa=1;
                        bb=nunique;
                        nunique_del=0;

                        while nunique_del < (nunique-n) && logical(count2 < 100)
                            count2=count2+1;
                            ndraw_del=unique([ndraw_del' round((aa+bb)/2+((bb+0.5)-(aa-0.5))*(-0.5+rand(num_draws, 1)))']');
                            nunique_del=length(ndraw_del);
                        end
                        p = randperm(nunique_del);
                        delete_ix=ndraw_del(p);

                        save_ix=setdiff(aa:bb, delete_ix(1:(nunique-n)));
                        ndraw=ndraw(save_ix);

                    end

                end

            end


            % If the number of integers to be selected is greater than half the
            % size of the draw pool of then use a deselection method.
            if isequal(flag2, 1)
                if length(ndraw) < 2
                    ndraw=setdiff((a:b)', ndraw');
                else
                    ndraw=setdiff((a:b)', ndraw')';
                end
            end

            if ~isequal(sorted, 1)
                draw_pool=randperm(length(ndraw))';
                ndraw=ndraw(draw_pool);
            end

            if fail
                if length(ndraw)~=n3
                    warning('Returning empty array');
                    ndraw=[];
                    error=1;
                else
                    % The error has been succesfully corrected 
                    %warning('Error Canceled.  Correct Number of Integers Selected.  Returned full array');
                    error=0;
                end
            end

        else

            % For the case of allowing duplicates

            if isequal(flag2, 1)
                n=(b-a+1)-n;
            end

            ndraw=round((a+b)/2+((b+0.5)-(a-0.5))*(-0.5+rand(n, 1)));

            % make sure the numbers fit in the proper range
            ndraw(ndraw > b)=b;
            ndraw(ndraw < a)=a;

            if isequal(sorted, 1)
                % sort output
                ndraw=sort(ndraw);
            end

        end

    end
end

%
% ************************************************************************
%
% Reshape the output array to the desired dimensions
% If there are not enough data points, then
% determine the best size for the output.
%
if logical(nn >= 1) && logical(mm >= 1)
    if logical( numel(ndraw) >= mm*nn)
        % reshape the data to the specified size
        ndraw=reshape(ndraw, mm, nn);
    else
        % Not enough data points to
        % reshape the data to the specified size.
        %
        % Get the number of data points available
        % and reshape the output to a reasonable size.
        np=numel(ndraw);
        if np < 1
            ndraw=[];
        else

            % Find the smaller dimension size
            [md min_dim]=min([mm, nn]);

            % determine the best way to reshape the output matrix
            % try to preserve the value of the smaller dimension
            if isequal(min_dim, 1)
                if mm <= np
                    nn=floor(np/mm);
                else
                    mm=np;
                    nn=1;
                end
            else

                if nn <= np
                    mm=floor(np/nn);
                else
                    nn=np;
                    mm=1;
                end
            end

            ndraw=reshape(ndraw(1:(mm*nn)), mm, nn);
            
        end

    end

end

% Plot the output if desired
% Good for inspecting the randomness of the output.
if isequal(make_plot, 1) && ~isempty(ndraw)
    
    figure(1);
    delete(1);
    figure(1);
    [m1, n1]=size(ndraw);
    
    if m1 > n1
        ndraw2=ndraw';
        axis_dir='Row';
    else
        ndraw2=ndraw;
        axis_dir='Column';
    end
    
    [m1, n1]=size(ndraw2);
    
    plot(1:n1, ndraw2, 'linestyle', 'none', 'Marker', '.', 'MarkerEdgeColor', 'k');
    ylim([a b]);
    
    xlim([0.5 n1+0.5]);
    
    regime_desc={'Build Entire Array, Then Discard Extras', 'Pure Purge', 'Select 10% Extra ,Then Discard Extras'};
    
    title({'Plot of Random Integers', ['Regime ', num2str(flag3), ' ', regime_desc{flag3}]}, 'Fontsize', 18);
    xlabel([axis_dir,' Indices of Integers'], 'Fontsize', 16);
    ylabel('Value of Integers', 'Fontsize', 16);
    
end

Contact us