findpointchazhi.m

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

function [pa,pf]=findpointchazhi(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];
  ma=max(f);
  for j=1:length(s)         %S中第J频点小于F中最大频率，则为有效频率 
    if s(j)<=ma
         ss(j)=s(j);
         pa(j)=interp1(f',a,ss(j),'nearest');
         pf(j)=s(j);
       
    else                %S中第J频点大于F中最大频率，则付0
               pa(j)=0;        %此频点的对应的数为0
               pf(j)=s(j);     %找出此频点
    end
 end
end