function avgPeaks = harmonicTest(n, trials)
   % Evaluate whether average number of peaks is Harmonic number
   % USAGE:  harmonicTest(n, trails)
   %         n is the length of a random vector
   %         trials is number of independent tests
   %            (trials is 1000 by default)
   %         avgPeaks is the overall average per trial
  
   if (nargin < 2)
       trials = 1000;    % DEFAULT choice
   end
   
   p = 0;   % number of peaks
   for t = 1:trials
       v = rand(1,n);
       p = p + peaks(v);
   end
   avgPeaks = p / trials;