localcorrelation.m

本人编写的利用小波算法进行滤波和局部相关运算

function y=positionalg()
load data2.mat;
lengths=50e3;
pos1=[];
sigma1=[];
corlength=1000;
corstep=corlength/2;
cornumber=261;

for i=14
s1=data2((i-1)*lengths+1:i*lengths,1);
s2=data2((i-1)*lengths+1:i*lengths,2);
s1=s1(1e4:corlength+1e4);
s2=s2(1e4:corlength+1e4);
% %-------------normal correlation ------------------- 
% s1=s1/sqrt(sum(s1.^2)/lengths);
% s2=s2/sqrt(sum(s2.^2)/lengths);
% [s1,cx,lx]=wden(s1,'minimaxi','s','one',5,'haar');
% s2=wden(s2,'minimaxi','s','one',5,'haar');
% 
% for k=corlength:corstep:lengths
%     local1=s1(k-corlength+1:k);
%     local2=s2(k-corlength+1:k);
%     y=xcorr(local1,local2);
%     [out,n]=max(abs(y));
% subplot(3,1,1), plot([local1 -local2]);
% grid;    
% subplot(3,1,2), plot(y);
% grid;    
% 
%  pos1=[pos1 (corlength-n)]
%  hold off;
% subplot(3,1,3), plot(pos1,'o');
% pos1;
% end
%--normal end ---------------------------------------

%---- wavelet correlation ------------------------------
scalea=100*2.^(0:0.5:10)';
s1wavelet=cwt(s1,scalea,'haar');
s2wavelet=cwt(s2,scalea,'haar');
%plot([s1wavelet(11,:) s2wavelet(11,:)])
for j=1:corlength
    s1factor(j)=trapz(scalea,(s1wavelet(:,j).*conj(s1wavelet(:,j)))./scalea); %integrate the normalze factor
    s2factor(j)=trapz(scalea,(s2wavelet(:,j).*conj(s2wavelet(:,j)))./scalea); %integrate the normalze factor
end
for j=1:length(scalea)
    s1nomalized(j,:)=s1wavelet(j,:)./s1factor(j);
    s2nomalized(j,:)=s2wavelet(j,:)./s2factor(j);
end
for j=1:corlength
    for k=max(1,j-cornumber):min(corlength,j+cornumber)
        scorrelated(j,k)=trapz(scalea,(s1nomalized(:,j).*s2nomalized(:,k))./scalea);
    end
end
for j=1:corlength
    [a,n]=max(scorrelated(:,j));
    pos(j)=j-n;
end
subplot(2,1,1),plot([s1 s2]);
subplot(2,1,2),plot(pos);
% for k=corlength:corstep:lengths
%     local1=s1(k-corlength+1:k);
%     local2=s2(k-corlength+1:k);
%     y=xcorr(local1,local2);
%     [out,n]=max(abs(y));
% subplot(3,1,1), plot([local1 -local2]);
% grid;    
% subplot(3,1,2), plot(y);
% grid;    
% 
%  pos1=[pos1 (corlength-n)]
%  hold off;
% subplot(3,1,3), plot(pos1,'o');
% pos1;
% end

%-----------Wavelet end ------------------------

pos=max(abs(pos1))
sigma1
subplot(2,1,1), plot(pos1);
%axis([0 140 -280 280]);
grid;
subplot(2,1,2), plot([s1 -s2]);
grid;
end