homef11x1copy.m

通过时空变换的方法对超声导波模态提取

%模态声发射定位:时空变换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('时间（ms）');
ylabel('幅度V');
title('原始信号1波形');
subplot(212);
plot(t,x2);
xlabel('时间（ms）');
ylabel('幅度V');
title('原始信号2波形');

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

%zhx111=xcorr(x1,x1);
%figure;
%plot((zhx111));
%滤波器设计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('信号1的频谱p');
subplot(212);
plot(fx,abs(pfx2));
xlabel('频率（kHz）');
ylabel('幅度');
title(' 信号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(' 理论计算波数');
p=polyfit(fre,boshu,3);
kf=polyval(p,fx);
subplot(212);
plot(fx,kf);
title(' 实际应用波数');

%由接收传感器波形,预算其他位置波形(导波传播理论)
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;


%模态提取算法
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('相关程度1');
figure;
plot(dis1,abs(kksh));
xlabel('距离（m）');
ylabel('相关程度2');

figure;
plot(dis1,kkksh);
xlabel('距离（m）');
ylabel('相关程度3');

%直接相关结果
zhx=xcorr(x1,x2);.0
figure;
xtt=-(N-1)*dt/1000:dt/1000:(N-1)*dt/1000;
plot(xtt,zhx);
xlabel('时间（ms）');
ylabel('相关程度4');