time   calls  line
                  1 function [alpha,fknee,white,f,measured,fitted,exitflag] = fitOneOverF(fs,timestream)
                  2     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                  3     %
                  4     % Fit a white noise plus 1/f power spectrum form
                  5     % to the given timestream.
                  6     % The fitting is moderately robust
                  7     %
                  8     % [alpha,fknee,white,f,measured,fitted] = fitOneOverF(fs,timestream)
                  9     %
                 10     % fs - sample frequency
                 11     % timestream - data stream
                 12     %
                 13     % alpha = 1/f^alpha index
                 14     % fknee is 1/f knee frequency
                 15     % white is standard deviation of white noise part of spectrum
                 16     % f is frequency axis
                 17     % measured is meaured power spectrum
                 18     % fitted is fitted power spectrum
                 19     %
                 20     %
                 21     % NB: This will fail if the knee frequency is higher than the sampling
                 22     % frequency, but you're pretty much screwed in that case anyway
                 23     
                 24     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                 25     %We fit to the log of the power spectrum.
                 26     %First we normalize the power spectrum to mean 0, std deviation 1
  0.16     860   27     stdev= std(timestream); 
  0.03     860   28     mu = mean(timestream); 
           860   29     d = (timestream-mu)/stdev; 
                 30     
                 31     %Next we compute the power spectrum to which we will fit.
                 32     %Its better to fit in log space.
  0.27     860   33     [f,fourier] = ffto(d,fs); 
  0.01     860   34     fourier= fourier/sqrt(length(fourier)); 
           860   35     fourier=fourier(f>0); 
           860   36     f=f(f>0); 
           860   37     measured=abs(fourier).^2; 
           860   38     logP = log(measured); 
  0.01     860   39     measured=measured*(stdev^2); 
                 40     
                 41     %The actual fitting code.  Seems fairly robust.
                 42     %Fit using the nested function residual, below.
                 43     %The initial state has white noise one and reasonable alpha and fknee
                 44     %given the data.
  0.01     860   45     x0 = [1., 0.05*(fs/100), 1.]; 
           860   46     parameter_lower_bounds=[0.5 f(1) 0]; 
           860   47     parameter_upper_bounds=[10.0 f(end) 100]; 
                 48     
  2.72     860   49 	options = optimset('MaxFunEvals',3000,'MaxIter',2000,'Display','off','Algorithm','interior-point'); 
                 50 %	[x,fval, exitflag] = fminsearch(@(y)residual(y),x0,options);
 211.44     860   51     [x,fval,exitflag] = fmincon(@(y)residual(y),x0,[],[],[],[],parameter_lower_bounds,parameter_upper_bounds,[],options); 
           860   52     x0_white = [1.]; 
  2.62     860   53 	options = optimset('MaxFunEvals',3000,'MaxIter',2000,'Display','off','Algorithm','interior-point'); 
                 54     
  3.36     860   55 	[x_white,fval_white, exitflag_white] = fminsearch(@(y)residualWhite(y),x0_white,options); 
                 56 
                 57 
           860   58     if (fval_white<fval) 
            30   59         exitflag=exitflag_white; 
            30   60         fitted = ones(size(f))*(x_white(1)*stdev)^2; 
            30   61         alpha=1; 
            30   62         fknee=0; 
            30   63         white=x_white*stdev; 
           830   64     else 
           830   65         alpha = x(1); 
  0.01     830   66         fknee = x(2); 
           830   67         x(3)=x(3)*stdev; 
           830   68         white = x(3); 
  0.08     830   69         fitted = oneOverFSpectrum(f,x); 
                 70 
           830   71     end 
                 72 
                 73     
                 74     
                 75     %Extract the parameters and denormalize.
                 76 
                 77     %It might be nice to return this too.
                 78   
                 79     %Because we are fitting to the logarithm, and to the power spectrum
                 80     %at that, we have to include the magic Euler-Mascharoni number.
                 81     %Read Dunkeley et al.
                 82     %astro-ph/0405462 
                 83     function r = residual(y)
                 84        logFitSpectrum=log(oneOverFSpectrum(f,y));       
                 85        difference = (logFitSpectrum - logP) - 0.5772156649015328606065;
                 86        r=sum(difference .* difference);
                 87        if (y(2)<0 || y(2)>fs) 
                 88            r=r+1e30;
                 89        end
                 90        if (y(3)<0)
                 91         r=r+1e30;
                 92        end
                 93     end
                 94 
                 95     function r = residualWhite(y)
                 96        logFitSpectrum=log(y(1));       
                 97        difference = (logFitSpectrum - logP) - 0.5772156649015328606065;
                 98        r=sum(difference .* difference);
                 99     end
                100 
                101 
                102 
           860  103 end 

Other subfunctions in this file are not included in this listing.