Code covered by the BSD License

# Continuous Sound and Vibration Analysis

### Edward Zechmann (view profile)

09 Sep 2008 (Updated )

This program analyzes sound and vibrations data using metrics for continuous noise and vibrations.

[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
% % 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
% %
% % **********************************************************************
% %
% %
% % 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
% %                                 Removed pure rejection regime warnings
% %                                 since the catch is reliable.
% %
% % modified   4 March      2009    Updated Comments.
% %
% %
% % **********************************************************************
% %
% % Feel free to modify this code.
% %
% %

% 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
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

```