findxchazhi.m

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

function [pa,pf]=findxchazhi(a,f)
%用插值法找与S匹配的频点，无效频点付0，先求出有效频点再一次插值，效率最高
  s=[6.250e+000 7.500e+000 8.750e+000 1.000e+001 1.125e+001 1.250e+001 1.375e+001 1.500e+001 1.625e+001 1.750e+001 1.875e+001 2.000e+001 2.125e+001 2.250e+001 2.375e+001 2.500e+001 2.625e+001 2.750e+001 2.875e+001 3.000e+001 3.125e+001 3.250e+001 3.375e+001 3.500e+001 3.625e+001 3.750e+001 3.875e+001 4.000e+001 4.125e+001 4.250e+001 4.375e+001 4.500e+001 4.625e+001  4.750e+001 4.875e+001 5.000e+001 5.125e+001 5.250e+001 5.375e+001 5.500e+001 5.625e+001 5.750e+001 5.875e+001 6.000e+001 6.125e+001 6.250e+001 6.375e+001 6.500e+001 6.625e+001 6.750e+001 6.875e+001 7.000e+001 7.125e+001 7.250e+001 7.375e+001 7.500e+001 7.625e+001 7.750e+001 7.875e+001 8.000e+001 8.125e+001 8.250e+001 8.375e+001 8.500e+001 8.625e+001 8.750e+001 8.875e+001 9.000e+001 9.125e+001 9.250e+001 9.375e+001 9.500e+001 9.625e+001 9.750e+001 9.875e+001 1.000e+002 1.125e+002 1.250e+002 1.375e+002 1.500e+002 1.625e+002 1.750e+002 1.875e+002 2.000e+002 2.125e+002 2.250e+002 2.375e+002 2.500e+002 2.625e+002 2.750e+002 2.875e+002 3.000e+002 3.125e+002 3.250e+002 3.375e+002 3.500e+002 3.625e+002 3.750e+002 3.875e+002 4.000e+002 4.125e+002 4.250e+002 4.375e+002 4.500e+002 4.625e+002 4.750e+002 4.875e+002 5.000e+002 5.125e+002 5.250e+002 5.375e+002 5.500e+002 5.625e+002 5.750e+002 5.875e+002 6.000e+002 6.125e+002 6.250e+002 6.375e+002 6.500e+002 6.625e+002 6.750e+002 6.875e+002 7.000e+002 7.125e+002 7.250e+002 7.375e+002 7.500e+002 7.625e+002 7.750e+002 7.875e+002 8.000e+002 8.125e+002 8.250e+002 8.375e+002 8.500e+002 8.625e+002 8.750e+002 8.875e+002 9.000e+002 9.125e+002 9.250e+002 9.375e+002 9.500e+002 9.625e+002 9.750e+002 9.875e+002 1.000e+003 1.125e+003 1.250e+003 1.375e+003 1.500e+003 1.625e+003 1.750e+003 1.875e+003 2.000e+003 2.125e+003 2.250e+003 2.375e+003 2.500e+003 2.625e+003 2.750e+003 2.875e+003 3.000e+003 3.125e+003 3.250e+003 3.375e+003 3.500e+003 3.625e+003 3.750e+003 3.875e+003 4.000e+003 4.125e+003 4.250e+003 4.375e+003 4.500e+003 4.625e+003 4.750e+003 4.875e+003 5.000e+003 5.125e+003 5.250e+003 5.375e+003 5.500e+003 5.625e+003 5.750e+003 5.875e+003 6.000e+003 6.125e+003 6.250e+003 6.375e+003 6.500e+003 6.625e+003 6.750e+003 6.875e+003 7.000e+003 7.125e+003 7.250e+003 7.375e+003 7.500e+003 7.625e+003 7.750e+003 7.875e+003 8.000e+003 8.125e+003 8.250e+003 8.375e+003 8.500e+003 8.625e+003 8.750e+003 8.875e+003 9.000e+003 9.125e+003 9.250e+003 9.375e+003 9.500e+003 9.625e+003 9.750e+003 9.875e+003 1.000e+004 1.125e+004 1.250e+004 1.375e+004 1.500e+004 1.625e+004 1.750e+004 1.875e+004 2.000e+004 2.125e+004 2.250e+004 2.375e+004 2.500e+004 2.625e+004 2.750e+004 2.875e+004 3.000e+004 3.125e+004 3.250e+004 3.375e+004 3.500e+004 3.625e+004 3.750e+004 3.875e+004 4.000e+004 4.125e+004 4.250e+004 4.375e+004 4.500e+004 4.625e+004 4.750e+004 4.875e+004 5.000e+004 5.125e+004 5.250e+004 5.375e+004 5.500e+004 5.625e+004 5.750e+004 5.875e+004 6.000e+004 6.125e+004 6.250e+004 6.375e+004 6.500e+004 6.625e+004 6.750e+004 6.875e+004 7.000e+004 7.125e+004 7.250e+004 7.375e+004 7.500e+004 7.625e+004 7.750e+004 7.875e+004 8.000e+004 8.125e+004 8.250e+004 8.375e+004 8.500e+004 8.625e+004 8.750e+004 8.875e+004 9.000e+004 9.125e+004 9.250e+004 9.375e+004 9.500e+004 9.625e+004 9.750e+004 9.875e+004 1.000e+005];
  pa=zeros(1,length(s));
   
  mi=min(f);
  m=find(s<mi);
  ma=max(f);
  n=find(s>ma);  %S中频点在F范围的则为有效频率 
  if length(m)~=0   %若M不为空，即F中的最小频率大于S中的最小频率
      ss=s(m(end)+1:n(1)-1); %先求出有效频点
      pa(1:m(end))=0;
      pa(m(end)+1:n(1)-1)=interp1(f',a,ss,'nearest'); %再一次性插值
  else
      ss=s(1:n(1)-1); %先求出有效频点
      pa(1:n(1)-1)=interp1(f',a,ss,'nearest'); %再一次性插值   
  end
  pa(n(1):end)=0;        %此频点的对应的数为0
  pf=s;     
 end