Документ взят из кэша поисковой машины. Адрес оригинального документа : http://lnfm1.sai.msu.ru/grav/russian/lecture/filtr/labs/Lab6.m Дата изменения: Thu Mar 24 18:01:25 2011 Дата индексирования: Mon Oct 1 23:41:10 2012 Кодировка: Поисковые слова: http astrokuban.info astrokuban

%singular spectral analysis
% program is written 17.02.2009 by L.V. Zotov
% do not forget about hank.m file


clear;
addpath('/home/leonid/lectures/labs/');% path to hank.m file
N_signal=1024;
% generating two-sin signal
for (k=1:1:N_signal)
garm(k)=k/50+5*sin(2*pi/20*(k-1))+k/100*sin(2*pi/200*(k-1));
end;

plot(garm);


% ARMA process generating
ar(1)=0.5*randn(1);
ar(2)=-0.2*ar(1)+0.5*randn(1);

for (i=3:1:N_signal)
ar(i)=0.9*ar(i-1)-0.7*ar(i-2)+0.5*randn(1);
end;

plot(ar);

%making a signal
signal=garm+5*ar;

plot(signal);

%--------------------------------------------------------------------------
cd '/home/leonid/lectures/labs/lab6';

fout=fopen('signal.dat','wt');
for(i=1:1:N_signal)
fprintf(fout,'%12.6f%12.6f\n',i,signal(i));


end;
fclose(fout);

% parameter - caterpillar length
L=220;

%matrix
K=N_signal-L;
A=zeros(L,K);
for ii=1:1:L
for j=1:1:K
A(ii,j)=signal(ii+j-1);
end;
end;

%decomposition
[U,S,V] = svd(A);
K_max=max(L,K)
L_min=min(L,K)


%10 components evaluation
for(i=1:1:10)
clear SX;
clear G;
SX=zeros(L,N_signal-L);
SX(i,i)=(S(i,i));
G=U*SX*V';
RX(i,:)=hank(G,L,N_signal,L_min,K_max);
end;

%weights calculation
for(i=1:1:N_signal)
if(i<=L_min)
w(i)=i;
elseif(i<=K_max)
w(i)=L_min;
else
w(i)=N_signal-i;
end;
end;

WK=zeros(10,10);
for(i=1:1:10)
for(j=1:1:10)
for(k=1:1:N_signal)
WK(i,j)=WK(i,j)+RX(i,k)*w(k)*RX(j,k);
end;
end;
end;

%WK(1,1)=0;
[I,J] = meshgrid([1:1:10]);
pcolor(I,J,WK)

%plot components
subplot(10,1,1); plot(RX(1,:)); title('1')
subplot(10,1,2); plot(RX(2,:)); title('2')
subplot(10,1,3); plot(RX(3,:)); title('3')
subplot(10,1,4); plot(RX(4,:)); title('4')
subplot(10,1,5); plot(RX(5,:)); title('5')
subplot(10,1,6); plot(RX(6,:)); title('6')
subplot(10,1,7); plot(RX(7,:)); title('7')
subplot(10,1,8); plot(RX(8,:)); title('8')
subplot(10,1,9); plot(RX(9,:)); title('9')
subplot(10,1,10); plot(RX(10,:)); title('10')

%grouping
plot(signal,'black');
hold on;
plot(RX(1,:),'red')
plot(RX(2,:)+RX(3,:))
plot(RX(4,:)+RX(5,:),'green')
hold off;

all_ssa=RX(1,:);
for(i=2:1:10)
all_ssa=all_ssa+RX(i,:);
end;

ost=signal-all_ssa;


hold on
plot(ost,'black')
plot(RX(6,:)+RX(7,:)+RX(8,:)+RX(9,:)+RX(10,:),'cyan')
plot(ar,'magenta')