ormsby.m

这是一个关于hht变换很有用的工具箱

function [XX]=ormsby(X,RF,TF,DT,N)
 
% The function ORMSBY returns the time domain filtered data.
% Non MATLAB Library routine used: MBSSEL.
%
% Example-
%	[acc]=ormsby(xa,0.12,0.14,dt,m).
%
% Calling sequence-
% [XX]=ormsby(X,RF,TF,DT,N)
%
% Input-
%	x	    - column vector data x(n)
%	rf,tf	- numbers describing the window
%	dt	    - the true delta t
%	n       - the dimension of data x
% Output-
%	xx	    - column vector of filtered data xx(n) 
%
% Used by-
%	ACC2DIS

%     TIME DOMAIN ORMSBY FILTER
%
      NN=100 ;
      LL=100 ;
      PI=3.141592653 ;
      WP=2.*PI*RF ;
      WS=2.*PI*TF ;
      DW=WS-WP    ;
      for I=1:NN  ;
          WPN=WP*I*DT ;
          WSN=WS*I*DT ;
          AA=cos(WPN)-cos(WSN) ;
          HD(I)=AA/(DW*PI*I*I*DT*DT) ;
      end
   
%     READ TIME WINDOW BAND
      for I=1:LL ;
          W(I+1)=HD(I) ;
      end
      
      W(1)=2*(TF+RF)*0.5 ;
      LL1=LL+1           ;
      X1=5               ;
      W1(1)=1				 ;
      [F0]=MBSSEL(X1)    ;
      for I=2:LL1        ;
          X11=X1*sqrt(1-((I-1)/(LL1-1))^2)  ;
          [F1]=MBSSEL(X11)   ;
          W1(I)=F1/F0        ;
      end
%
%     SET WINDOW AREA EQUAL TO ONE
%
      for I=1:LL1        ;
          W1(I)=W(I)*W1(I)   ;
      end
      
      SM=W1(1)           ;
      for I=1:LL         ;
          SM=SM+W1(I+1)*2    ;
      end
      
      for I=2:LL1        ;
         W(LL-I+2)=W1(I)/SM ;
         W(LL+I)=W(LL-I+2)  ;
      end
   
      W(LL+1)=W1(1)/SM   ;
%
%     REFLECT DATA AT END AND BEGINNEING OF THE RECORD
%
      for I=1:N          ;
          Y(I+LL)=X(I)   ;
      end
      
      for I=1:LL         ;
      Y(I)=X(LL1-I)      ;
      Y(I+N+LL)=X(N-I+1) ;
      end
   
      for I=1:N                  ;
          XM=0.                  ;
          ILL=I+2*LL             ;
          for J=I:ILL	         ;
              M=J-I+1            ;
              XM=XM+Y(J)*W(M)    ;
          end
          X(I)=XM                ;
      end     
      
      XX=X ;