# arma2allan_v000

```%assumption: number of zeros must be strictly less than number of poles.
%(a0)x(n+m)+(a1)x(n+m-1)...+(am)x(n)=(b0)w(n+m)+(b1)w(n+m-1)...+(bk)a(n+m-k
%)
function allan_th=arma2allan(m_A, m_B, fs,t_axis, allan)

[nmod na]=size(m_A);
nb=size(m_B,2);
allan_th=zeros(nmod,length(t_axis));

%discrete indeces
n_axis=round(t_axis*fs);

%normalize the model so that m_A(1)=1
m_A=m_A./(m_A(:,1)*ones(1, size(m_A,2)));
m_B=m_B./(m_A(:,1)*ones(1, size(m_B,2)));

for in=1:nmod
%stip the excessive zeros for each model
A=m_A(in,:);
B=m_B(in,:);
sr_a=find(fliplr(A)~=0,1);
if (~isempty(sr_a))
A=A(1,1:end-sr_a+1);
end
sr_a=find(fliplr(B)~=0,1);
if (~isempty(sr_a))
B=B(1,1:end-sr_a+1);
end

if (size(A,2)==1) %WN or MA
if (size(B,2)==1) %WN
allan_th(in,:)=m_B(in,1)./((n_axis).^0.5);
else %MA
disp('error - To be completed later');
%%%%TOOOO LAZZZZZYYYY
end
elseif ((size(A,2)==2) && A(2)==(-1)) %RW
allan_th(in,:)=((2*(n_axis.^2)+1)./(6*n_axis)).^0.5*B(1);
else %strictly proper ARMA
%convert model into correlation exponentials
exp_param=arma2cov_v000(B, A);
allan_th(in,:)=exp2allan(exp_param,n_axis);
end
end

% if (~isempty(allan))
%     figure;
%     loglog(t_axis,allan);grid;
%     hold on;
%     for in=1:nmod
%         loglog(t_axis,allan_th(in,:),'k');
%     end
%     loglog(t_axis,sum(allan_th.^2,1).^0.5,'r');
% end

if (~isempty(allan))
figure;
loglog(t_axis,allan,'linewidth',2);
hold on;
loglog(t_axis,sum(allan_th.^2,1).^0.5,'r','linewidth',2);   %%just for legend correction
for in=1:nmod
loglog(t_axis,allan_th(in,:),'k','linewidth',2);
end
loglog(t_axis,sum(allan_th.^2,1).^0.5,'r','linewidth',2);
end
xlabel('Cluster Length (sec)');
grid;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function allan_res=exp2allan(ep,n_axis)
nexp=size(ep,1);
allan_res=zeros(1,length(n_axis));
for in=1:length(n_axis)
n=n_axis(in);
var_cluster=sum(ep(:,2))/n;
cor_cluster=0;
for ii=1:nexp
c=ep(ii,2);
v=ep(ii,1);
sr_a=2/n*c*v*(n-1-n*v+v^n)/n/(1-v)^2;
var_cluster=var_cluster+sr_a;

sr_a=c*(v^n)/n;
sr_b=c/(n^2)*v*(1-n*(v^(n-1))+(n-1)*(v^n))/(1-v)^2;
sr_c=c/(n^2)*v^(2*n-1)*(1-n*(v^(1-n))+(n-1)*(v^(-n)))/(1-v^(-1))^2;

sr_d=sr_a+sr_b+sr_c;
if (isnan(sr_d))    %this method is very problematic in terms of numerical precison for high "n" values
cor_cluster=cor_cluster+0;
else
cor_cluster=cor_cluster+sr_d;
end
end

allan_res(in)=(var_cluster-cor_cluster)^0.5;
end

```