fitOneOverF

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 Fit a white noise plus 1/f power spectrum form
 to the given timestream.
 The fitting is moderately robust

 [alpha,fknee,white,f,measured,fitted] = fitOneOverF(fs,timestream)

 fs - sample frequency
 timestream - data stream

 alpha = 1/f^alpha index
 fknee is 1/f knee frequency
 white is standard deviation of white noise part of spectrum
 f is frequency axis
 measured is meaured power spectrum
 fitted is fitted power spectrum


 NB: This will fail if the knee frequency is higher than the sampling
 frequency, but you're pretty much screwed in that case anyway

0001 function [alpha,fknee,white,f,measured,fitted,exitflag] = fitOneOverF(fs,timestream)
0002     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0003     %
0004     % Fit a white noise plus 1/f power spectrum form
0005     % to the given timestream.
0006     % The fitting is moderately robust
0007     %
0008     % [alpha,fknee,white,f,measured,fitted] = fitOneOverF(fs,timestream)
0009     %
0010     % fs - sample frequency
0011     % timestream - data stream
0012     %
0013     % alpha = 1/f^alpha index
0014     % fknee is 1/f knee frequency
0015     % white is standard deviation of white noise part of spectrum
0016     % f is frequency axis
0017     % measured is meaured power spectrum
0018     % fitted is fitted power spectrum
0019     %
0020     %
0021     % NB: This will fail if the knee frequency is higher than the sampling
0022     % frequency, but you're pretty much screwed in that case anyway
0023     
0024     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0025     %We fit to the log of the power spectrum.
0026     %First we normalize the power spectrum to mean 0, std deviation 1
0027     stdev= std(timestream);
0028     mu = mean(timestream);
0029     d = (timestream-mu)/stdev;
0030     
0031     %Next we compute the power spectrum to which we will fit.
0032     %Its better to fit in log space.
0033     [f,fourier] = ffto(d,fs);
0034     fourier= fourier/sqrt(length(fourier));
0035     fourier=fourier(f>0);
0036     f=f(f>0);
0037     measured=abs(fourier).^2;
0038     logP = log(measured);
0039     measured=measured*(stdev^2);
0040     
0041     %The actual fitting code.  Seems fairly robust.
0042     %Fit using the nested function residual, below.
0043     %The initial state has white noise one and reasonable alpha and fknee
0044     %given the data.
0045     x0 = [1., 0.05*(fs/100), 1.];
0046     parameter_lower_bounds=[0.5 f(1) 0];
0047     parameter_upper_bounds=[10.0 f(end) 100];
0048     
0049     options = optimset('MaxFunEvals',3000,'MaxIter',2000,'Display','off','Algorithm','interior-point');
0050 %    [x,fval, exitflag] = fminsearch(@(y)residual(y),x0,options);
0051     [x,fval,exitflag] = fmincon(@(y)residual(y),x0,[],[],[],[],parameter_lower_bounds,parameter_upper_bounds,[],options);
0052     x0_white = [1.];
0053     options = optimset('MaxFunEvals',3000,'MaxIter',2000,'Display','off','Algorithm','interior-point');
0054     
0055     [x_white,fval_white, exitflag_white] = fminsearch(@(y)residualWhite(y),x0_white,options);
0056 
0057 
0058     if (fval_white<fval)
0059         exitflag=exitflag_white;
0060         fitted = ones(size(f))*(x_white(1)*stdev)^2;
0061         alpha=1;
0062         fknee=0;
0063         white=x_white*stdev;
0064     else
0065         alpha = x(1);
0066         fknee = x(2);
0067         x(3)=x(3)*stdev;
0068         white = x(3);
0069         fitted = oneOverFSpectrum(f,x);
0070 
0071     end
0072 
0073     
0074     
0075     %Extract the parameters and denormalize.
0076 
0077     %It might be nice to return this too.
0078   
0079     %Because we are fitting to the logarithm, and to the power spectrum
0080     %at that, we have to include the magic Euler-Mascharoni number.
0081     %Read Dunkeley et al.
0082     %astro-ph/0405462
0083     function r = residual(y)
0084        logFitSpectrum=log(oneOverFSpectrum(f,y));       
0085        difference = (logFitSpectrum - logP) - 0.5772156649015328606065;
0086        r=sum(difference .* difference);
0087        if (y(2)<0 || y(2)>fs) 
0088            r=r+1e30;
0089        end
0090        if (y(3)<0)
0091         r=r+1e30;
0092        end
0093     end
0094 
0095     function r = residualWhite(y)
0096        logFitSpectrum=log(y(1));       
0097        difference = (logFitSpectrum - logP) - 0.5772156649015328606065;
0098        r=sum(difference .* difference);
0099     end
0100 
0101 
0102 
0103 end
0104 
0105 
0106 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0107 function Ys = oneOverFSpectrum(f,params) 
0108 % Functional form of the spectrum to fit to the data
0109 alpha = params(1);
0110 fknee = params(2);
0111 white = params(3);
0112 
0113     Ys = (1.0 +(fknee./f).^alpha)*(white^2);
0114 end
0115 
0116