oneftest.m

%  This M-file creates a filter that, when excited by white noise, will
%  generate a time series that has Spectral Density that is proportional
%  to 1/f, where f is the frequency.
%
%  Start by setting the poles to be an octave apart and the
%  zeros symmetric with respect to the poles (half way between the poles).

numoctaves = 45;      % The number of octaves that the Spectrum has 1/f.
w = logspace(-1,12);  % For the psd plot.
alpha = .5;           % For setting the zero location.
omegalow = 2*pi*0.1;  % The lowest frequency is 0.1 Hz.

n      = 2.^(0:numoctaves-1); % Get the spacing of the poles.
p      = omegalow*n;          % Poles spaced an octave apart (2^n*omegalow)
z      = alpha*p;             % Zeros half way between poles.

A = diag(-p);                 % Create the A matrix for the system
for i = 1:numoctaves-1        % The lower triangular part of the A matrix
    A = A+diag(p(i:end-1)-z(i:end-1),-i);
end

B = (z-p)';                   % Create the B vector
C = ones(size(B))';           % Create the C vector
d = 1;

sys_1 = ss(A,B,C,d);          % Use the Control System Toolbox to get lti
bode(sys_1,w)                 %   object sys.  This is used to get the 
grid                          %   Bode plot.