clear
n=1;
for t=0:0.1:2
    x=(exp(-0.8*t))*220e3*sin(2*pi*1.9*t)+(exp(0.3*t))*220e3*sin(2*pi*0.5*t);
    y(n)=x;
    n=n+1;
end
t=0:0.1:2;
x=y;
dt=0.1;
%去直流分量
L = length(x);
s=0;
for i=1:L
s=s+x(i);
end
xp=s/L;
for i=1:L
x(i)= x(i)-xp;
end
%去除奇異分量（階越，下跌，膨大）
y=x;
x=[];
for i=1:2
    x(i)=y(i);
end
for i=3:L
    x(i)=(y(i-2)-4*y(i-1)+3*y(i))/2/dt;
end
%%去噪聲
%
%y=[];
%y=x;
%x=[];
%x=hilbert(y);
%低通濾波
[n,fo,mo,w]=remezord([2.5 3],[1 0],[0.01 0.1],10);
b=remez(n,fo,mo,w);
y=fftfilt(b,x);
dt = 0.1;
s=x;
L = length(s);
pe = ceil(L/2);
R  = [];
k = 0;   
% 形成擴(kuò)展階矩陣
for i=1:pe
    for j=1:pe
        R(i,j)=0;
        for n=pe+1:L
            R(i,j)= R(i,j)+s(n-j)*conj(s(n-i));
        end
    end
end    
[U,S,V] = svd(R);  % 分解
[m,n] = size(R);
h = min(m,n);
cro = [];
for ks = 1:h;
    cro = [cro,S(ks,ks)];
end
ccro = cro/cro(1);
% 估計(jì)階數(shù)
for pp = 1:h
    if ccro(pp)>=1e-6
        pp = pp + 1;
    else 
        break;
    end
end
p = pp-1;
RE = R(2:end,2:p+1);
b  = R(:,1);
b = b(2:end);
a = pinv(RE)*(-b);
R1 = R(1,:);
R1 = R1(1:p+1);
a1 = [1 a'];
Ep = sum(R1.*a1);
P = [1 a'];
z = roots(P)
%% 估計(jì)序列X
kk = 1:p;
X(kk) = s(kk);
for pp = p+1:L
    ll = 1:p;
  X(pp) = sum(-a'.*s(pp-ll));
    pp = pp+1;
end
Zh = [];
for mm = 0:L-1
    Zh = [Zh,z.^mm];
end
Z = Zh';
Z = conj(Z);
Zhh  = Z';
b=2*inv(Zhh*Z)*Zhh*X';
A=abs(b)
f=angle(z)/2/pi/dt
a=log(abs(z))/dt
theta = angle(b)/pi*180
plot(X,':r');
hold on;
plot(s,'-.b');
hold on;