psbosubak.m

打开一原始图像

if fid==-1
            depth=depth+2;
            finm=strcat(str1,num2str(depth),str2);
            fid=fopen(finm);
            dd=(depth-1)/10;
            data(1:1024,dd)=fread(fid,'uint16')/2^15-dd;
        else
            depth=depth+9;
        end
        dd=(depth+1)/10;
        data(1:1024,dd)=fread(fid,'uint16')/2^15-dd;
%cc(1:5120)=0;
%cc(3584:(3584+1024-1))=a1;
ls=1:1024;
%ls=ls/51.200;
figure(1);
%plot(ls,a1,'k',ls,a2-1,'k',ls,a3-2,'k',ls,a4-3,'k',ls,a5-4,'k',ls,a6-5,'k');
title('1＃钻孔98.0m～99.0m深度地层波速测量实测曲线')
xlabel('采样时间（ms）');


%c1=cwt(a1,1:64,'db4');
wl='db4';
c1=cwt(a1,1:64,wl);
c2=cwt(a2,1:64,wl);
c3=cwt(a3,1:64,wl);
c4=cwt(a4,1:64,wl);
c5=cwt(a5,1:64,wl);
c6=cwt(a6,1:64,wl);


scale1=12;
figure(9);
%plot(ls,c1(scale1,:),'k',ls,c2(scale1,:)'-2,'k',ls,c3(scale1,:)'-4,'k',ls,c4(scale1,:)'-6,'k',ls,c5(scale1,:)'-8,'k',ls,c6(scale1,:)'-10,'k');
%plot(ls,abs(c1(scale1,:)),'k',ls,abs(c2(scale1,:))-2,'k',ls,c3(scale1,:)'-4,'k',ls,abs(c4(scale1,:))-6,'k',ls,abs(c5(scale1,:))-8,'k',ls,c6(scale1,:)'-10,'k');
title('原始数据经小波变换后的小波系数曲线(a=56)')
xlabel('采样时间（ms）');

scale2=12;
figure(10);
lsp=1024/2;
%lspt=(1:lsp)/51.2;
lspt=(1:lsp);
%plot(lspt,c1(scale2,1:lsp)',lspt,c2(scale2,1:lsp)'-2,lspt,c3(scale2,1:lsp)'-4,lspt,c4(scale2,1:lsp)'-6,lspt,c5(scale2,1:lsp)'-8,lspt,c6(scale2,1:lsp)'-10);
%plot(lspt,c1(scale2,1:lsp),'k',lspt,c2(scale2,1:lsp)-2,'k',lspt,c3(scale2,1:lsp)-4,'k',lspt,c4(scale2,1:lsp)-6,'k',lspt,c5(scale2,1:lsp)-8,'k',lspt,c6(scale2,1:lsp)-10,'k');
%plot(lspt,(c1(scale2,1:lsp)),'k',lspt,c2(scale2,1:lsp)-2,'k',lspt,c3(scale2,1:lsp)-4,'k',lspt,c4(scale2,1:lsp)-6,'k',lspt,c5(scale2,1:lsp)-8,'k',lspt,c6(scale2,1:lsp)-10,'k');
title('原始数据经小波变换后的小波系数曲线(a=8)')
xlabel('采样时间（ms）');

figure(11);
c1=cwt(a1,1:64,wl,'plot');
figure(12);
c1=cwt(a4,1:64,wl,'plot');


figure(2);
subplot(3,1,1);
%plot([a1,c1(scale1,:)']);
subplot(3,1,2);
%plot([a2,c2(scale2,:)']);
subplot(3,1,3);
%plot([a1,a2]);
%figure(3);
x=1:1024;
xx=x/51200;
%plot(x,c1(a,:)','r',x,c2(a,:)','b');
d1=[a1,c1(40,:)'];
d2=[a2,c2(40,:)'];
d3=[a1,a2];
d4=[c1(40,:)', c2(40,:)'];


sf=1024;
sf2=sf/2;
Y=fft(a1,sf);%参加DFT的采样点个数；
fa1=Y.*conj(Y)/sf;%功率谱密度；
f=51200*(0:(sf2-1))/sf;

/*figure(4);
%plot(f,fa1(1:sf2));
set(gca,'XTick',1000:1000:4000);
xqq=1000;xww=4000;yqq=0;yww=1;axis([xqq xww yqq yww]);%x轴和y轴的区间
x=zeros(1024,2);
x(:,1)=a1;
x(:,2)=a2;*/
%figure(5);
%plot(x);


save a1 a1;
save a2 a2;
save d1.dat d1 -ascii
save d2.dat d2 -ascii
save d3.dat d3 -ascii
save d4.dat d4 -ascii