%   This script will simulate a 1000 Hz tone (for 16 ms),
%   add white Gaussian noise, take its FT, and test to see if
%   1000 Hz has the greatest magnitude.  If so, this is considered
%   a detection.  The script compute 30 probabilities of detection
%   using 1000 simulations per run.  The mean and the standard
%   deviation of the probability of detection is computed in
%   the end

runs = 30;          %  Number of runs
runlen = 1000;      %  Number of simulations per run
fs = 4e3;           %  Sampling Frequency
ftone = 1e3;        %  Tone Singal Frequency
t = [0:63]/fs;      %  Time axis (0 to 16 ms).
sig = sin(2*pi*ftone*t);   %  Signal (doesn't change)

%  Big  Loop to compute a new 'pd' (probability of detection)

for n=1:runs,
  d=0;             %  detection counter

%  Inner Loop to run each simulation

  for k=1:runlen
    nsig = sig + .9988*randn(1,64);     %  Signal plus noise
    spec = abs(fft(nsig));              %  FT of signal plus noise
    if(spec(17) == max(spec(1:32)));    %  The 17-th FT component
                                        %   corresponds to 1000 Hz
                                        %   If it is the max spectral
                                        %   component, then it counts as
                                        %   a detection.
      d=d+1;
    end
  end
  pd(n) = d/runlen;          %  Compute fraction of times it exceeded threshold       
end
mpd = mean(pd);            %  Compute mean Probability of Detection
spd = std(pd);             %  Compute standard Deviation of PD