eh4pl.m

来自「基于FFT的大地电磁信号资料处理MATLAB方法实现」· M 代码 · 共 74 行

function  qq=eh4pl(x)
%对EH4高频进行EMD分解并对各模态分量进行FFT变化
t=0.3/4096:0.3/4096:0.3;
f =emd(x);
imf=f';
tn=length(t);

t1=min(t(1),t(tn));
t2=max(t(1),t(tn));

%%%%%%%%%%%%%%%%%%% Plot the EMD %%%%%%%%%%%%%%%%%%%
figure;
clf;

imfStart = 1;
imfStop = size(imf,2);
imfLast = imfStop+1;

for j=imfStart:imfLast;
    subplot(imfLast-imfStart+1,1,j+1-imfStart);
    if j==imfStart
        fftpl(x,1);
            else
       fftpl(imf(:,j-1),1);
          end
      s1='imf';
   if j==1
       ylabel('signal')
   else
   s2=[s1 num2str(j-1)];
   ylabel(s2);
   end
   if( j == 1 )
       title( 'Spectrum of Signal and IMF Components' );
   end
   if( j < imfLast );
      set(gca,'xticklabel','');
   end;
end;
xlabel('f/Hz');

figure;
for j=imfStart:imfLast;
    subplot(imfLast-imfStart+1,1,j+1-imfStart);
    if j==imfStart
        z=sum(imf(:,imfStart:imfStop),2);
        plot(t,z);
        yMax = max(z);
        yMin = min(z);
    else
        plot(t,imf(:,j-1));
        yMax = max(imf(:,j-1));
        yMin = min(imf(:,j-1));
    end
    if( abs(yMax-yMin) > 1e-10 )
        axis([t1 t2 yMin yMax]);
    else
        axis([t1 t2 -1 1]);
    end
   s1='imf';
   if j==1
       ylabel('signal')
   else
   s2=[s1 num2str(j-1)];
   ylabel(s2);
   end
   if( j == 1 )
       title( 'Time Series and IMF Components' );
   end
   if( j < imfLast );
      set(gca,'xticklabel','');
   end;
end;
xlabel('Time');