ytxh.m

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

function imf=ytxh(x,cy)
%取得X文件对应的频点，作频谱图(高频）
o=fft(x,length(x));
a=o.*conj(o)/length(x);
f=cy*(0:round(length(x)/2))/length(x);
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];
   for j=1:292
       for i=1:2048
           k(i)=abs(f(i)-s(j));
       end
           l=min(k);
           n(j)=find(k==l);
   end
   
   for i=1:292
pa(i)=a(n(i));
pf(i)=f(n(i));
   end
   plot(pf,pa);
end