fftzzhijie.asv

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

function [xsp,z]=fftzfile()
%FFT直接最小二乘法
%获得一测点的292个频点的XFILE,输出XFILE数据，放在不同的文件下，每个文件放一页(即一种采样频率的数据）
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,9, ndv);
cj=1;
for ii=1:ns               %每一块（12288点）计数一此XFILE
    if s(ii,1)~=dv(cj)
        cj=cj+1;
    end
    xp=zeros(292,9);
   cy=192000/2^dv(cj);
                             %某一块某一段的时间
    for jj=1:3                %每一段（4096点）计数一此XFILE
        duan=blocksj((jj-1)*4096+1:jj*4096,1:4,ii);
        xpc=fftxzhijie(duan,cy);%得到各频点的复信号
        xp(:,1)=xpc(:,1);
        xp(:,2:9)= xp(:,2:9)+xpc(:,2:9);
    end
    xsp(:,1,cj)=xp(:,1);
    xsp(:,2:9,cj)=xsp(:,2:9,cj)+xp(:,2:9);
end
%各频率求平均
xsp(:,2:9,1)=xsp(:,2:9,1)/24;
xsp(:,2:9,2)=xsp(:,2:9,2)/24;
xsp(:,2:9,3)=xsp(:,2:9,3)/18;
%放在不同的文件下，每个文件放一页
for j=1:3
s1='Xfftzj';
s2=num2str(j);
s3='.txt';
s=[s1,s2,s3];
fid=fopen(s,'w');
for i=1:292
fprintf(fid,'%15e',real(xsp(i,:,j)));%复信号实部和虚部分行存储
fprintf(fid,'\r\n');
fprintf(fid,'%15e',imag(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,6) pp(i,7) 0 0
           pp(i,8) pp(i,9) 0 0
            0 0 pp(i,6) pp(i,7)
           0 0 pp(i,8) pp(i,9)];
       b=[pp(i,2) pp(i,3) pp(i,4) pp(i,5)];
       [z1,flag]=gaujor(a,b);
             
       z(i,:,j)=[xsp(i,1,j) z1']; %z中保存复数ZXX,ZXY ZYX ZYY
   end
end
%ZFILE放在不同的文件下，每个文件放一页
for j=1:cj
    s1='Zfftzj';
    s2=num2str(j);
    s3='.txt';
    s=[s1,s2,s3];
    fid=fopen(s,'w');
    for i=1:292
        fprintf(fid,'%15e',real(z(i,:,j)));%Z的实部和虚部分行存储
        fprintf(fid,'\r\n');
        fprintf(fid,'%15e',imag(z(i,:,j)));
        fprintf(fid,'\r\n');
    end
    fclose(fid);
end


end