% 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);