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