x3file.m

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

function xsp=x3file()
%获得一测点的292个频点的XFILE,输出XFILE数据，放在不同的文件下，每个文件放一页(即一种采样频率的数据）
clear all;
clc;
%读取EH摊时间序列Y文件数据
[filename,pathname]=uigetfile('Y*.*','Please input a data file');
EHyfile=[pathname,filename];
fid=fopen(EHyfile,'r');
x=fread(fid,'int16');

%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               %每一块（12288点）计数一此XFILE
    if s(ii,1)~=dv(cj)
        cj=cj+1;
    end
    xp=zeros(292,19);
    t1=1/(192000/2^dv(cj))*4096;
    t=linspace(0,t1,4096);    %某一块某一段的时间
    for jj=1:3                %每一段（4096点）计数一此XFILE
        duan=blocksj((jj-1)*4096+1:jj*4096,1:4,ii);
        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;
%放在不同的文件下，每个文件放一页
for j=1:3
s1='Xx';
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


end