fit1overf

 [a,A,B,fknee] = fit1overf(V,fs,popt)
 This function fits a spectrum of the form:
 y = A*f^(-a) + B
 to the spectrum of the samples in V. fs is the sample frequency in Hz.

 If popt == 1 then plot the fitted spectrum

 The function returns:
   a: the power of f, nominally 1.
   A: the scale factor for 1/f part of spectrum
   B: the white noise average power level
   fknee: the knee frequency in Hz of the spectrum
       fknee = exp(1/a*log(A/(sqrt(2)*B)))

 ogk

0001 function [a,A,B,fknee] = fit1overf(V,fs,popt)
0002 % [a,A,B,fknee] = fit1overf(V,fs,popt)
0003 % This function fits a spectrum of the form:
0004 % y = A*f^(-a) + B
0005 % to the spectrum of the samples in V. fs is the sample frequency in Hz.
0006 %
0007 % If popt == 1 then plot the fitted spectrum
0008 %
0009 % The function returns:
0010 %   a: the power of f, nominally 1.
0011 %   A: the scale factor for 1/f part of spectrum
0012 %   B: the white noise average power level
0013 %   fknee: the knee frequency in Hz of the spectrum
0014 %       fknee = exp(1/a*log(A/(sqrt(2)*B)))
0015 %
0016 % ogk
0017 %
0018 
0019 % Step 1: Get the spectrum of the data
0020 [f,Y] = ffto(V,fs);
0021 Y = Y(f>0)/sqrt(length(f));
0022 f = f(f>0);
0023 c = {f abs(Y)};
0024 
0025 % Starting point for minimization:
0026 x0 = [1;0.05;std(V-mean(V))/2];
0027 %x0 = [1;30;140];
0028 options = optimset('MaxFunEvals',1000,'MaxIter',1000);
0029 [x,fval] = fminsearch(@(x) ssqErr(x,c),x0,options);
0030 
0031 a = x(1);
0032 A = x(2);
0033 B = x(3);
0034 fknee = exp(1/a*log(A/((sqrt(2)-1)*B)));
0035 
0036 if popt == 1
0037     
0038     %figure
0039    loglog(f,abs(Y),'k-',f,A*f.^(-a)+B,'r-',...
0040    [fknee fknee],[min(abs(Y)) max(abs(Y))],'g-.',...
0041    [min(f) max(f)],[sqrt(2)*B sqrt(2)*B],'g-.')
0042     ylim([min(abs(Y)) max(abs(Y))])
0043     xlim([min(f) max(f)])
0044     xlabel('Frequency [Hz]')
0045     ylabel(sprintf('Power spectrum: f_{knee} = %4.3f Hz',fknee))
0046     legend('Data','Fit')
0047     
0048 end
0049 
0050 end
0051 
0052 function err = ssqErr(x,c)
0053 % The error function we want to minimize
0054 Y = c{2};
0055 f = c{1};
0056 
0057 Yfit = spec(f,x);
0058 err = sum((Yfit-Y).^2);
0059 % Minimize in log space:
0060 % err = sum((log10(Yfit)-log10(Y)).^2);
0061 end
0062 
0063 function Ys = spec(f,x)
0064 % Functional form of the spectrum to fit to the data
0065 a = x(1);
0066 A = x(2);
0067 B = x(3);
0068 Ys = A.*f.^(-a)+B;
0069 end