% The function oneoverf_filt.m generates a unit-area discrete-time
% filter with a  1/(f^alpha) magnitude-squared frequency response
% over a prescribed bandwidth using the Keshner/Saletti method.
% In particular, a continuous-time filter with a magnitude-squared 
% frequency response of 1/(f^alpha) is designed, and then it is 
% converted to discrete-time using the bilinear transformation.  
%
% Note that if the product Omega*Ts is too small, there will be
% numerical issues when creating a and b.  Also note that dmin = -4 
% and Ts = 0.0625 provides sound results.  These parameters should
% be fine for short-term data collection from simulate.m
%
% Function arguments:
%	th - current parameters values
%	Sgran - desired granularity
%	dmin - minimum frequency of desired 1/f behavior in decades
%	dmax - maximum frequency of desired 1/f behavior in decades
%
% Function output:
%	h - unit-area, discrete-time filter with 1/(f^alpha)
%	    magnitude-squared frequency response over a
%	    prescribed bandwidth
%

function h = oneoverf_filt(th,Sgran,dmin,dmax)

% Assigning variables.
alpha = th(59);
dnum = dmax-dmin;

% Defining continous-time poles and zeros in units of rad/sec.
i = 0:dnum-1;
Omegap = -(10^dmin)*(10^((1-alpha/2)/2))*(10.^i)*2*pi;
Omegaz = -(10^dmin)*(10^((1-alpha/2)/2))*(10^(alpha/2))*(10.^i)*2*pi;

falloff = 0;
Omegap = [Omegap Omegap(dnum)*10*ones(1,falloff)];

% Using the BLT to transform to discrete-time poles and zeros.
omegap = ((1+Omegap*(Sgran/2))./(1-Omegap*(Sgran/2)));
omegaz = ((1+Omegaz*(Sgran/2))./(1-Omegaz*(Sgran/2)));

% Calculating ARMA coefficients from the discrete-time poles and zeros.
a = 1;
b = 1;

for i = 1:dnum

	b = conv(b,[1 -omegaz(i)]);

end

b = [zeros(1,falloff) b]*sqrt((2*pi)^alpha);

for i = 1:dnum+falloff

	a = conv(a,[1 -omegap(i)]);

end

% Calculating moving average filter from ARMA coefficients.
% The length of the filter is determined such that the last
% sample of the created impulse response is less than 2.5
% percent of the first sample. 
ratio = 0.1;
count = 0;
while (ratio > 0.025)

      count = count+20;
      h = filter(b,a,[1; zeros(count,1)]);
      ratio = h(end)/h(1);

end

% Normalizing function output to unity area.
h = h/sum(h);