function [upper, lower ] = specconf(nrecords)
%
%    The function [upper, lower ] = specconf(nrecords)
%    estimates the 95% upper and lower confidence intervals 
%    if for the spectrum S there was ensemble averaging of 
%    nrecords spectra.
%
%  ---------------------------------------------------------------
%    Written by Vassilis Zervakis, while at
%               AINCO-Interocean, Madrid, Spain, in June 1995.
%  ---------------------------------------------------------------
%
% --- the matrix chi contains the chi squared values for ------
% --- 0.05 and 0.95 confidence intervals (per column) for -----
% --- the averaging number n provided in the 3rd column ------
% -----Reference Bendat & Piercol -----------------------------
chi = [  0.0039   3.84   1;...
         0.103    5.99   2;...
         0.352    7.81   3;...
         0.711    9.49   4;...
         1.15    11.07   5;...
         1.64    12.59   6;...
         2.17    14.07   7;...
         2.73    15.51   8;...
         3.33    16.92   9;...
         3.94    18.31  10;...
         4.57    19.68  11;...
         5.23    21.03  12;...
         5.89    22.36  13;...
         6.57    23.68  14;...
         7.26    25.00  15;...
         7.96    26.30  16;...
         8.67    27.59  17;...
         9.39    28.87  18;...
        10.12    30.14  19;...
        10.85    31.41  20;...
        11.59    32.67  21;...
        12.34    33.92  22;...
        13.09    35.17  23;...
        13.85    36.42  24;...
        14.61    37.65  25;...
        15.38    38.88  26;...
        16.15    40.11  27;...
        16.93    41.34  28;...
        17.71    42.56  29;...
        18.49    43.77  30;...
        26.51    55.76  40;...
        43.19    79.08  60;...
        95.70   146.57 120];
if nrecords <= 30,
   upper = nrecords / chi(nrecords,2);
   lower = nrecords / chi(nrecords,1);
end
if (30 < nrecords) & (nrecords <= 40),
   upperslope = (chi(31,2)-chi(30,2)) / 10.;
   lowerslope = (chi(31,1)-chi(30,1)) / 10.;
   chiup = upperslope * (nrecords-30) + chi(30,2);
   chidn = lowerslope * (nrecords-30) + chi(30,1);
   upper = nrecords / chiup;
   lower = nrecords / chidn;
end
if (40 < nrecords) & (nrecords <= 60),
   upperslope = (chi(32,2)-chi(31,2)) / 20.;
   lowerslope = (chi(32,1)-chi(31,1)) / 20.;
   chiup = upperslope * (nrecords-40) + chi(31,2);
   chidn = lowerslope * (nrecords-40) + chi(31,1);
   upper = nrecords / chiup;
   lower = nrecords / chidn;
end
if (60 < nrecords) & (nrecords <= 120),
   upperslope = (chi(33,2)-chi(32,2)) / 60.;
   lowerslope = (chi(33,1)-chi(32,1)) / 60.;
   chiup = upperslope * (nrecords-60) + chi(32,2);
   chidn = lowerslope * (nrecords-60) + chi(32,1);
   upper = nrecords / chiup;
   lower = nrecords / chidn;
end
if nrecords > 120,
   n = nrecords;
   chiup = n*(1 - 2/(9*n) + 0.95*sqrt(2/(9*n)))^(1/3);
   chidn = n*(1 - 2/(9*n) + 0.05*sqrt(2/(9*n)))^(1/3);
   upper = nrecords / chiup;
   lower = nrecords / chidn;
end
