SVD-TLS-ARMA-Code

11 views (last 30 days)
Coo Boo
Coo Boo on 5 Sep 2012
Hi to all; I'm confused with some parts of the code which is about the use of singular value decomposition- total least squares estimation of AR parameters of ARMA model. The code and related discriptions are as follows:
%%%Initializing
clear
loops = 20;
M = 100;
pe =30;
p=0;
for loop=1:loops
for k = 1:50
n = [1:128];
w = randn(1,128);
x = sqrt(20)*sin(2*pi*0.2*n) + sqrt(2)*sin(2*pi*0.213*n)+w;
randn(1,128) x = sqrt(20)*sin(2*pi*0.2*n) + sqrt(2)*sin(2*pi*0.213*n)+randn(1,128)
Rxx = xcorr(x,'unbiased');
%%%%Constructing Hankel Matrix
for i = 1:M
for j = 1:pe+1
Re(i,j) = Rxx(pe+i+1-j);
end
end
%%%%%%applying SVD
[U,S,V]=svd(Re);
Ak = 0;
for i=1:pe+1
Ak = Ak + S(i,i)^2;
end
%%%%%Calculating the true order p
Akf=0;
v = Akf/Ak;
p=0;
while v<0.9979
p = p + 1;
Akf = Akf + S(p,p)^2;
v = Akf/Ak;
end
P(k) = p;
end
P;
p = fix(mean(P));
%%%%%%!!!!!! I don't know what does the following part do:!!!!!
%%%%Is there any refernce about this part? What is the basis?
Sp = 0;
for j=1:p
for i=1:pe+1-p
Vj=V(:,j);
Sp=Sp+S(j,j)^2*(Vj([i:i+p]))*(Vj([i:i+p]))';
end
end
invSp=inv(Sp);
for i = 1:p
a(i)=invSp(i+1,1)/invSp(1,1);
end
%%%%%%%the following parts are about finding roots and ARMA parameters
ARMA_LS(:,loop) = a;
Az(1) = 1;
for i = 2:p+1
Az(i) = a(i-1);
end
z = roots(Az);
num = 1;
for i = 1:2:p
f_ls(num) = atan(imag(z(i))/real(z(i)));
f_ls(num) = abs(f_ls(num))/(2*pi);
num = num +1;
end
f_ls = sort(f_ls);
Fz(:,loop) = f_ls;
end
p
ARMA_LS
ARMA_LS_mean = mean(ARMA_LS,2)
ARMA_LS_std = std(ARMA_LS')'
Fz
Fz_mean = mean(Fz,2)
Fz_var = var(Fz')'

Answers (0)

Categories

Find more on Discrete Multiresolution Analysis in Help Center and File Exchange

Tags

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!