%内含2种最优小波基计算
clear;
close all;
clc;
load('M.mat');
[max_M,index]=max(M);%找到每列最大值和每列最大值所在的行数
%[phi,g1,xval]=wavefun('haar',10);%获得参数，以得到小波时域图
for i=1:7
M(:,i)=M(:,i)/max_M(i);
end%归一化
[phi,g1,xval]=wavefun('db2',15);%建议使用10次以上迭代计算，比较精确
[m,p]=max(g1);%找到最大值和最大值所在的行数
d=fix(length(g1)/1000);%带余除法取商，确定采样间隔
a=mod(p,d);%确定能取到最大值的首个位置
if a==0
a=d;
end
l=(g1(a:d:length(g1))/m)';%采样后的小波,归一化
wavelet=l([1:1000]);%取前1000个点
[m,p]=max(wavelet);
plot(wavelet);
%相似性计算
%b=circshift(a',[2,0]);%列向量向下循环位移2个
% for i=1:7
% b=circshift(wavelet,[index(i)-p,0]);
% k=corrcoef(b,M(:,i));
% gamma(i)=k(1,2)%相关计算
% end
s=0;x1=0;x2=0;m1=mean(M);m2=mean(wavelet);
for i=1:7
b=circshift(wavelet,[index(i)-p,0]);
for k=1:1000
s=s+abs(M(k,i)-m1(i))*abs(wavelet(k)-m2);
x1=x1+(M(k,i)-m1(i))^2;
x2=x2+(wavelet(k)-m2)^2;
end
%方法1
gamma(i)=abs(s)/sqrt(x1*x2);
%方法2
%k=corrcoef(b,M(:,i));
%gamma(i)=k(1,2)%相关计算
end
% subplot(2,1,1);
% plot(xval,g1,'LineWidth',2);xlabel('t');title('db4?时域');
% g2=fft(g1);g3=abs(g2);
% subplot(2,1,2);
% plot(g3,'LineWidth',2);xlabel('f');title('db4?频域');
