clear
clc
clf;
dx=0.01;
x=0:dx:40.96-dx;
load data.dat;
a=HCN(1501:5596,2);
b=HCN(1501:5596,5);
dal=dtrend(a);
ddal=dtrend(dal,1);
ddall=idfilt(ddal,50,[0.1/50 25/50]);
dbl=dtrend(b);
ddbl=dtrend(dbl,1);
ddbll=idfilt(ddbl,50,[0.1/50 25/50]);
data=iddata(ddbll,ddall,0.01);
na=32;nb=33;nk=delayest(data);
th=arx(data,'na',na,'nb',nb,'nk',nk);
data1=predict(th,data);
e1=pe(th,data);
[A,B]=polydata(th);
for i=1:length(B)
if B(1)==0
B=B(2:length(B));
end
end
[z,p,k]=tf2zp(B,A);
for i=1:length(p)
f(i)=(log(p(i)).*log(p(i)')).^0.5./(2*pi*dx);
h(i)=-1*(log(p(i))+log(p(i)'))./(4*pi*f(i)*dx);
end
f=smooth(f);
h=smooth(h);
subplot(3,1,1),plot(x,a,'-',x,b,':');
subplot(3,1,2),plot(f);
subplot(3,1,3),plot(h); |