zfile.asv

基于hht的大地电磁信号资料处理MATLAB方法实现

function [xsp,z]=zfile()
%获得一测点的292个频点的ZFILE,输出XFILE数据和ZFILE，放在不同的文件下，每个文件放一页(即一种采样频率的数据）
clear all;
clc;
%读取EH摊时间序列Y文件数据
%EHyfile=[pathname,filename];
%fid=fopen(EHyfile,'r');
fid=fopen('C:\MATLAB701\work\YMD206.018','r');
x=fread(fid,'int16');
fclose(fid);

%fid=fopen('C:\MATLAB701\work\YMD206.030','r');
%x=fread(fid,'int16');
%fclose(fid);
%将数据分段
L=length(x);
s5=x(5);s6=x(6);s7=x(7);
dic=find((x==s6));

j=1;
for ik=1:length(dic)
if((x(dic(ik)-1)==s5)&(x(dic(ik)+1)==s7))
  a(j)=dic(ik)-5;
   j=j+1;
   end
end
a(j)=L+1;
%分出每道12288个数据
s=zeros(j-1,7);
HY=zeros(s6,1);
EX=zeros(s6,1);
HX=zeros(s6,1);
EY=zeros(s6,1);
blocksj=zeros(s6,4,j-1);
%数据头信息
for i=1:j-1
s(i,:)=x(a(i):(a(i)+6));
HY=x((a(i)+7):4:(a(i+1)-4));
EX=x((a(i)+8):4:(a(i+1)-3));
HX=x((a(i)+9):4:(a(i+1)-2));
EY=x((a(i)+10):4:(a(i+1)-1));
%读出的数据块
blocksj(:,:,i)=[HY,EX,HX,EY];
end
clear x dic a HY EX HX EY L j i s5 s7 fid;
ns=size(s,1);
dec=diff(s(:,1));
ix=find(dec);
nix=length(ix);
    for i=1:nix
      dv(i)=s(ix(i),1);
      dv(i+1)=s(ix(i)+1,1);
    end
ndv=length(dv);
xsp=zeros(292,19, ndv);
cj=1;
for ii=1:ns
    if s(ii,1)~=dv(cj)
        cj=cj+1;
    end
    xp=zeros(292,19);
    for jj=1:3
        duan=blocksj((jj-1)*4096+1:jj*4096,1:4,ii);
        t1=1/(192000/2^dv(cj))*4096;
        t=linspace(0,t1,4096);
        xpc=getxfile(duan,t);
        xp(:,1:3)=xpc(:,1:3);
        xp(:,4:19)= xp(:,4:19)+xpc(:,4:19);
    end
    xsp(:,1:3,cj)=xp(:,1:3);
    xsp(:,4:19,cj)=xsp(:,4:19,cj)+xp(:,4:19);
end
%各频率求平均
xsp(:,4:19,1)=xsp(:,4:19,1)/24;
xsp(:,4:19,2)=xsp(:,4:19,2)/24;
xsp(:,4:19,3)=xsp(:,4:19,3)/18;
%XFILE放在不同的文件下，每个文件放一页
for j=1:cj
   s1='Xx84';
   s2=num2str(j);
   s3='.txt';
   s=[s1,s2,s3];
   fid=fopen(s,'w');
   for i=1:292
       fprintf(fid,'%15e',xsp(i,:,j));
       fprintf(fid,'\r\n');
   end
   fclose(fid);
end

%各频率的Z值
for j=1:cj
   pp=xsp(:,:,j);
   for i=1:292
       a=[pp(i,14) 0 pp(i,12) pp(i,6)
           0 pp(i,14) -pp(i,6) pp(i,12)
           pp(i,12) -pp(i,6) pp(i,4) 0
           pp(i,6) pp(i,12) 0 pp(i,4)];
       b=[pp(i,13) -pp(i,10) pp(i,8) pp(i,5)];
       [z1,flag]=gaujor(a,b);%z1保存ZXX,ZXY的实部和虚部
       
       a=[pp(i,14) 0 pp(i,12) pp(i,6)
           0 pp(i,14) -pp(i,6) pp(i,12)
           pp(i,12) -pp(i,6) pp(i,4) 0
           pp(i,6) pp(i,12) 0 pp(i,4)];
       b=[pp(i,18) pp(i,15) pp(i,16) pp(i,7)];
       [z2,flag]=gaujor(a,b); %z2保存ZXX,ZXY的实部和虚部
       z(i,:,j)=[xsp(i,1,j) z1' z2'];
   end
end
%ZFILE放在不同的文件下，每个文件放一页
for j=1:cj
    s1='Zz84';
    s2=num2str(j);
    s3='.txt';
    s=[s1,s2,s3];
    fid=fopen(s,'w');
    for i=1:292
        fprintf(fid,'%15e',z(i,:,j));
        fprintf(fid,'\r\n');
    end
    fclose(fid);
end
end