%模態(tài)聲發(fā)射定位:時(shí)空變換modeF(1,1)
clear
N=2048;    
%采樣率為kHz
sample=200;                                                              
dt=1/sample;
t=0:dt:(N-1)*dt;
df=sample/N;
f=0:df:(N-1)*df;
%A3-53;A3-17;A3-57;
fid1=fopen('E:\jiao\data\40 pipe\8-29\a3distance\A3leak45\da1.001','r');
fid2=fopen('E:\jiao\data\40 pipe\8-29\a3distance\A3leak45\da1.002','r');
x1=fscanf(fid1,'%g',[N,1]);
x2=fscanf(fid2,'%g',[N,1]);
fclose(fid1);
fclose(fid2);

%泄漏原始圖形
figure;
subplot(211);
plot(t,x1);
xlabel('時(shí)間（ms）');
ylabel('幅度V');
title('原始信號(hào)1波形');
subplot(212);
plot(t,x2);
xlabel('時(shí)間（ms）');
ylabel('幅度V');
title('原始信號(hào)2波形');

%頻譜分析
fx1=fft(x1);
fx2=fft(x2);
figure;
subplot(211);
plot(f(1:N/2),abs(fx1(1:N/2)));
xlabel('頻率（kHz）');
ylabel('幅度');
title(' 信號(hào)1的頻譜');
subplot(212);
plot(f(1:N/2),abs(fx2(1:N/2)));
xlabel('頻率（kHz）');
ylabel('幅度');
title(' 信號(hào)2的頻譜');

%zhx111=xcorr(x1,x1);
%figure;
%plot((zhx111));
%濾波器設(shè)計(jì)100,20
N1=120;
NL=10;
MM=N1-NL+1;
for ii=1:MM;
    pfx1(ii)=fx1(NL+ii-1);
    pfx2(ii)=fx2(NL+ii-1);
    fx(ii)=f(NL+ii-1);
end;
figure;
subplot(211);
plot(fx,abs(pfx1));
xlabel('頻率（kHz）');
ylabel('幅度');
title('信號(hào)1的頻譜p');
subplot(212);
plot(fx,abs(pfx2));
xlabel('頻率（kHz）');
ylabel('幅度');
title(' 信號(hào)2的頻譜p');

%彌散曲線的確定
NN=80;
fid1=fopen('E:\jiao\data\40 pipe\disper txt\ff11.txt','r');
fid2=fopen('E:\jiao\data\40 pipe\disper txt\kf11.txt','r');
fre1=fscanf(fid1,'%g',[NN,1]);
boshu=fscanf(fid2,'%g',[NN,1]);
fclose(fid1);
fclose(fid2);
figure;
subplot(211);
fre=fre1;
plot(fre,boshu);
title(' 理論計(jì)算波數(shù)');
p=polyfit(fre,boshu,3);
kf=polyval(p,fx);
subplot(212);
plot(fx,kf);
title(' 實(shí)際應(yīng)用波數(shù)');

%由接收傳感器波形,預(yù)算其他位置波形(導(dǎo)波傳播理論)
ddis=10;
dis=-200*ddis:ddis:200*ddis;
M=length(dis);
for ii=1:M;
    for jj=1:MM;
    ppfx1(ii,jj)=pfx1(jj)*exp(i*kf(jj)*dis(ii)*2*pi);
    ppfx2(ii,jj)=pfx2(jj)*exp(i*kf(jj)*dis(ii)*2*pi);
end;
end;


%模態(tài)提取算法
for ii=1:MM;
    pvfx1(ii)=0;
    pvfx2(ii)=0;
    for jj=1:M;
        pvfx1(ii)=pvfx1(ii)+ppfx1(jj,ii)*exp(-i*kf(ii)*dis(jj)*2*pi)*ddis;
        pvfx2(ii)=pvfx2(ii)+ppfx2(jj,ii)*exp(-i*kf(ii)*dis(jj)*2*pi)*ddis;
    end;
end;

%定位模型

ddis2=500;
dis2=-65*ddis2:ddis2:0*ddis2;
M2=length(dis2);
for ii=1:M2;
   ksh(ii)=0 ;
   kksh(ii)=0;
   kkksh(ii)=0;
    for j=1:MM;
       ksh(ii)=ksh(ii)+pvfx1(j)*conj(pvfx2(j))*exp(-i*kf(j)*dis2(ii))*df;
    end;
    kksh(ii)=real(ksh(ii))/imag(ksh(ii));
    rr=real(ksh(ii));
    mmk=imag(ksh(ii));
    kkksh(ii)=sqrt(rr*rr+mmk*mmk);
end;
dis1=dis2./1000;

figure;
plot(dis1,abs(ksh));
xlabel('距離（m）');
ylabel('相關(guān)程度1');
figure;
plot(dis1,abs(kksh));
xlabel('距離（m）');
ylabel('相關(guān)程度2');

figure;
plot(dis1,kkksh);
xlabel('距離（m）');
ylabel('相關(guān)程度3');

%直接相關(guān)結(jié)果
zhx=xcorr(x1,x2);.0
figure;
xtt=-(N-1)*dt/1000:dt/1000:(N-1)*dt/1000;
plot(xtt,zhx);
xlabel('時(shí)間（ms）');
ylabel('相關(guān)程度4');