
function alpha=get_alpha(image,f1,f2)

% This function calculates the slope of the average amplitude
% spectrum across all orientations, and plots the output on a loglog scale.
% 
% Input:
% image: name of a NxM matrix representing an image (any image can be read 
% into this matrix using matlab's function "imread". The image
% can be rescaled to gray scale with matlab's "rgb2gray" function).
% 
% f1, f2 - frequency band between which fitting of the amplitude spectrum will 
% be done. The frequencies should be given in cycles per image.
% 
% Olga Dyakova, Uppsala University, Motion Vision Group, 2014
% ref. Bruno Olshausen function "rotavg", available online 
% at http://redwood.berkeley.edu/bruno/ VS212B/ lab2/rotavg.m
% and Dr Cyril Pernet, Fourier Analysis for neuroscientists A practical 
% guide using Matlab: http://www.sbirc.ed.ac.uk/cyril/download/DTP_Fourier_analyses.pdf


fft_image = fftshift(fft2(image));%perform Fourier transform and 
%move the zero-frequency component to the center
abs_im=abs(fft_image);
[N, M]=size(abs_im);
% find the rotational average
[X, Y]=meshgrid(-M/2:M/2-1,-N/2:N/2-1);
[~, rho]=cart2pol(X,Y); %get the polar coordinates
rho=round(rho);
abs_av=nan(1,N/2);
for r=0:N/2-1
  idx= rho==r;
  abs_av(r+1)=mean(abs_im(idx));%get the rotationally averaged ampl. spectrum
end

freq=0:N/2-1;
figure
loglog(freq,abs_av,'k') % plot the average amplitude spectrum on loglog scale
xlabel('spatial frequency')
ylabel('average amplitude spectrum')
xx=log(freq(f1:f2));         % cut low and high frequencies
yy=log(abs_av(freq(f1:f2))); % get rotationally averaged ampl. spectrum 
                             % data points for fitting

p=polyfit(xx,yy,1); % perform fitting
alpha=(-1)*p(1); % get the alpha
display(alpha);