taoboy 发表于 2008-11-12 21:48

麻烦大家告诉我nspab.m文件在哪才能找到啊??

我是个初学者,时频分析工具箱,EMD工具箱里面就是找不到这个文件,论坛里也转了很长时间也没找到,
哪位大侠帮我找找啊,万分感谢啊!~~

xubin 发表于 2008-11-13 00:48

nspab.m文件

function h1= nspab(data,nyy,minw,maxw,dt)
%----- Get dimensions (number of time points and components)
= size(data);

%----- Get time interval


%----- Apply Hilbert Transform
data=hilbert(data);
a=abs(data);
omg=abs(diff(unwrap(angle(data))))/(2*pi*dt);

%----- Smooth amplitude and frequency   
filtr=fir1(8,.1);
for i=1:knb
    a(:,i)=filtfilt(filtr,1,a(:,i));
   omg(:,i)=filtfilt(filtr,1,omg(:,i));
end
for i=1:knb
   a(:,i)=filtfilt(filtr,1,a(:,i));
omg(:,i)=filtfilt(filtr,1,omg(:,i));
end

%----- Limit frequency and amplitude
for i=1:knb
   for i1=1:npt-1
      if omg(i1,i) >=maxw,
         omg(i1,i)=maxw;
         a(i1,i)=0;
      elseif omg(i1,i)<=minw,
         omg(i1,i)=minw;
         a(i1,i)=0;
      else
      end
   end
end
clear filtr data

%va=var(omg(200:1200))

%----- Get local frequency
dw=maxw - minw;
wmx=maxw;
wmn=minw;

%----- Construct the ploting matrix
clear p;
h1=zeros(npt-1,nyy+1);
p=round(nyy*(omg-wmn)/dw)+1;
for j1=1:npt-1
   for i1=1:knb
      ii1=p(j1,i1);      
      h1(j1,ii1)=h1(j1,ii1)+a(j1,i1);
   end
end

%----- Do 3-point to 1-point averaging
=size(h1);
%n1=fix(nx/3);
%h=zeros(n1,ny);
%for i1=1:n1
   %h(i1,:)=(h1(3*i1,:)+h1(3*i1-1,:)+h1(3*i1-2,:));
%end
%clear h1;

%----- Do 3-points smoothing in x-direction
fltr=1./3*ones(3,1);
for j1=1:ny
h1(:,j1)=filtfilt(fltr,1,h1(:,j1));
end
clear fltr;


%----- Define the results
%w=linspace(wmn,wmx,ny-1)';
%xs=linspace(t0,t1,nx)';
h1=flipud(rot90(h1));
h1=h1(1:ny-1,:);

taoboy 发表于 2008-11-13 12:45

太谢谢你了啊,不知道说什么了。。送你鲜花吧,呵呵

taoboy 发表于 2008-11-13 14:43

提示这个错误是什么意思啊??

??? Error using ==> filtfilt
Data must have length more than 3 times filter order.

Error in ==> nspab at 16
    a(:,i)=filtfilt(filtr,1,a(:,i));

Error in ==> egg1 at 11
h=nspab(imf,250,0,250,dt);
页: [1]
查看完整版本: 麻烦大家告诉我nspab.m文件在哪才能找到啊??