emd.cpp

vc代用matcom编程

  #include "matlib.h"
  #pragma hdrstop
  
  #include "emd.h"
  #include "extr.h"
  #include "io.h"
  
  
  
  Mm emd(Mm x, i_o_t, Mm& imf__o, Mm& ort__o, Mm& nbits__o) {
    begin_scope
    x.setname("x"); 
    dMm(imf); dMm(ort); dMm(nbits); dMm(defstop); dMm(t); dMm(stop); dMm(tst); dMm(S); dMm(sd); dMm(sd2); dMm(tol) \
      ; dMm(MAXITERATIONS); dMm(NBSYM); dMm(lx); dMm(sdt); dMm(sd2t); dMm(ner); dMm(nzr); dMm(r); dMm(k); dMm(nbit) \
      ; dMm(NbIt); dMm(m); dMm(mp); dMm(sx); dMm(test); dMm(indmin); dMm(indmax); dMm(indzer); dMm(lm); dMm(lM); dMm( \
      nem); dMm(nzm); dMm(j_); dMm(lmax); dMm(lmin); dMm(lsym); dMm(rmax); dMm(rmin); dMm(rsym); dMm(tlmin); dMm(tlmax) \
      ; dMm(trmin); dMm(trmax); dMm(mlmax); dMm(mlmin); dMm(mrmax); dMm(mrmin); dMm(envmax); dMm(envmin); dMm(envmoy) \
      ; dMm(s); 
    
    #line 1 "d:/matcom45/samples/emd/emd.m"
    call_stack_begin;
    #line 1 "d:/matcom45/samples/emd/emd.m"
    // nargin, nargout entry code
    double old_nargin=nargin_val; if (!nargin_set) nargin_val=1.0;
    nargin_set=0;
    double old_nargout=nargout_val; if (!nargout_set) nargout_val=3.0;
    nargout_set=0;
    
    // translated code
    
    
    #line 3 "d:/matcom45/samples/emd/emd.m"
    // default for stopping
    #line 4 "d:/matcom45/samples/emd/emd.m"
_   defstop = (BR(0.05),0.5,0.05);
    #line 5 "d:/matcom45/samples/emd/emd.m"
_   t = colon(1.0,1.0,length(x));
    #line 6 "d:/matcom45/samples/emd/emd.m"
_   stop = defstop;
    #line 7 "d:/matcom45/samples/emd/emd.m"
_   tst = 0.0;
    
    #line 9 "d:/matcom45/samples/emd/emd.m"
_   S = size(x);
    #line 10 "d:/matcom45/samples/emd/emd.m"
_   if (istrue(((S(1.0)>1.0)&&(S(2.0)>1.0)))||istrue((length(S)>2.0))) {
      #line 11 "d:/matcom45/samples/emd/emd.m"
_     display( error(TM("x must have only one row or one column")) );
      #line 12 "d:/matcom45/samples/emd/emd.m"
    }
    
    #line 14 "d:/matcom45/samples/emd/emd.m"
_   if (istrue(S(1.0)>1.0)) {
      #line 15 "d:/matcom45/samples/emd/emd.m"
_     x = ctranspose(x);
      #line 16 "d:/matcom45/samples/emd/emd.m"
    }
    
    #line 18 "d:/matcom45/samples/emd/emd.m"
_   S = size(t);
    #line 19 "d:/matcom45/samples/emd/emd.m"
_   if (istrue(((S(1.0)>1.0)&&(S(2.0)>1.0)))||istrue((length(S)>2.0))) {
      #line 20 "d:/matcom45/samples/emd/emd.m"
_     display( error(TM("t must have only one row or one column")) );
      #line 21 "d:/matcom45/samples/emd/emd.m"
    }
    
    #line 23 "d:/matcom45/samples/emd/emd.m"
_   if (istrue(S(1.0)>1.0)) {
      #line 24 "d:/matcom45/samples/emd/emd.m"
_     t = ctranspose(t);
      #line 25 "d:/matcom45/samples/emd/emd.m"
    }
    
    #line 27 "d:/matcom45/samples/emd/emd.m"
_   if (istrue(length(t)!=length(x))) {
      #line 28 "d:/matcom45/samples/emd/emd.m"
_     display( error(TM("x and t must have the same length")) );
      #line 29 "d:/matcom45/samples/emd/emd.m"
    }
    
    #line 31 "d:/matcom45/samples/emd/emd.m"
_   S = size(stop);
    #line 32 "d:/matcom45/samples/emd/emd.m"
_   if (istrue(((S(1.0)>1.0)&&(S(2.0)>1.0)))||istrue((S(1.0)>3.0))||istrue((S(2.0)>3.0))||istrue((length(S)>2.0) \
      )) {
      #line 33 "d:/matcom45/samples/emd/emd.m"
_     display( error(TM("stop must have only one row or one column of max three elements")) );
      #line 34 "d:/matcom45/samples/emd/emd.m"
    }
    
    #line 36 "d:/matcom45/samples/emd/emd.m"
_   if (istrue(S(1.0)>1.0)) {
      #line 37 "d:/matcom45/samples/emd/emd.m"
_     stop = ctranspose(stop);
      #line 38 "d:/matcom45/samples/emd/emd.m"
_     S = size(stop);
      #line 39 "d:/matcom45/samples/emd/emd.m"
    }
    
    #line 41 "d:/matcom45/samples/emd/emd.m"
_   if (istrue(S(2.0)<3.0)) {
      #line 42 "d:/matcom45/samples/emd/emd.m"
_     stop(3.0) = defstop(3.0);
      #line 43 "d:/matcom45/samples/emd/emd.m"
    }
    
    #line 45 "d:/matcom45/samples/emd/emd.m"
_   if (istrue(S(2.0)<2.0)) {
      #line 46 "d:/matcom45/samples/emd/emd.m"
_     stop(2.0) = defstop(2.0);
      #line 47 "d:/matcom45/samples/emd/emd.m"
    }
    
    #line 49 "d:/matcom45/samples/emd/emd.m"
_   sd = stop(1.0);
    #line 50 "d:/matcom45/samples/emd/emd.m"
_   sd2 = stop(2.0);
    #line 51 "d:/matcom45/samples/emd/emd.m"
_   tol = stop(3.0);
    
    #line 53 "d:/matcom45/samples/emd/emd.m"
_   if (istrue(tst)) {
      #line 54 "d:/matcom45/samples/emd/emd.m"
_     display(figure());
      #line 55 "d:/matcom45/samples/emd/emd.m"
    }
    
    #line 57 "d:/matcom45/samples/emd/emd.m"
    // maximum number of iterations
    #line 58 "d:/matcom45/samples/emd/emd.m"
_   MAXITERATIONS = 2000.0;
    
    #line 60 "d:/matcom45/samples/emd/emd.m"
    // maximum number of symmetrized points for interpolations
    #line 61 "d:/matcom45/samples/emd/emd.m"
_   NBSYM = 2.0;
    
    #line 63 "d:/matcom45/samples/emd/emd.m"
_   lx = length(x);
    
    #line 65 "d:/matcom45/samples/emd/emd.m"
_   sdt(lx) = 0.0;
    #line 66 "d:/matcom45/samples/emd/emd.m"
_   sdt = sdt+sd;
    #line 67 "d:/matcom45/samples/emd/emd.m"
_   sd2t(lx) = 0.0;
    #line 68 "d:/matcom45/samples/emd/emd.m"
_   sd2t = sd2t+sd2;
    
    #line 70 "d:/matcom45/samples/emd/emd.m"
    // maximum number of extrema and zero-crossings in residual
    #line 71 "d:/matcom45/samples/emd/emd.m"
_   ner = lx;
    #line 72 "d:/matcom45/samples/emd/emd.m"
_   nzr = lx;
    
    #line 74 "d:/matcom45/samples/emd/emd.m"
_   r = x;
    #line 75 "d:/matcom45/samples/emd/emd.m"
_   imf = nop_M;
    #line 76 "d:/matcom45/samples/emd/emd.m"
_   k = 1.0;
    
    #line 78 "d:/matcom45/samples/emd/emd.m"
    // iterations counter for extraction of 1 mode
    #line 79 "d:/matcom45/samples/emd/emd.m"
_   nbit = 0.0;
    
    #line 81 "d:/matcom45/samples/emd/emd.m"
    // total iterations counter
    #line 82 "d:/matcom45/samples/emd/emd.m"
_   NbIt = 0.0;
    
    #line 84 "d:/matcom45/samples/emd/emd.m"
_   while (istrue(ner>2.0)) {
      
      #line 86 "d:/matcom45/samples/emd/emd.m"
      // current mode
      #line 87 "d:/matcom45/samples/emd/emd.m"
_     m = r;
      
      #line 89 "d:/matcom45/samples/emd/emd.m"
      // mode at previous iteration
      #line 90 "d:/matcom45/samples/emd/emd.m"
_     mp = m;
      
      #line 92 "d:/matcom45/samples/emd/emd.m"
_     sx = sd+1.0;
      
      #line 94 "d:/matcom45/samples/emd/emd.m"
      // tests if enough extrema to proceed
      #line 95 "d:/matcom45/samples/emd/emd.m"
_     test = 0.0;
      
      #line 97 "d:/matcom45/samples/emd/emd.m"
_     /*[indmin,indmax,indzer] = */extr(m,i_o,indmin,indmax,indzer);
      #line 98 "d:/matcom45/samples/emd/emd.m"
_     lm = length(indmin);
      #line 99 "d:/matcom45/samples/emd/emd.m"
_     lM = length(indmax);
      #line 100 "d:/matcom45/samples/emd/emd.m"
_     nem = lm+lM;
      #line 101 "d:/matcom45/samples/emd/emd.m"
_     nzm = length(indzer);
      
      #line 103 "d:/matcom45/samples/emd/emd.m"
_     j_ = 1.0;
      
      #line 105 "d:/matcom45/samples/emd/emd.m"
      // sifting loop
      #line 106 "d:/matcom45/samples/emd/emd.m"
_     while (istrue((mean(sx>sd)>tol||any(sx>sd2)||(abs(nzm-nem)>1.0)))&&istrue((test==0.0))&&istrue(nbit<MAXITERATIONS) \
        ) {
        
        #line 108 "d:/matcom45/samples/emd/emd.m"
_       if (istrue(nbit>MAXITERATIONS/5.0)&&istrue(mod(nbit,floor(MAXITERATIONS/10.0))==0.0)) {
          #line 109 "d:/matcom45/samples/emd/emd.m"
_         display( disp((BR(TM("mode ")),int2str(k),TM(" nombre d iterations : "),int2str(nbit))) );
          #line 110 "d:/matcom45/samples/emd/emd.m"
_         display( disp((BR(TM("stop parameter mean value : ")),num2str(S))) );
          #line 111 "d:/matcom45/samples/emd/emd.m"
        }
        
        #line 113 "d:/matcom45/samples/emd/emd.m"
        // boundary conditions for interpolations :
        
        #line 115 "d:/matcom45/samples/emd/emd.m"
_       if (istrue(indmax(1.0)<indmin(1.0))) {
          #line 116 "d:/matcom45/samples/emd/emd.m"
_         if (istrue(m(1.0)>m(indmin(1.0)))) {
            #line 117 "d:/matcom45/samples/emd/emd.m"
_           lmax = fliplr(indmax(colon(2.0,1.0,min(length(indmax),NBSYM+1.0))));
            #line 118 "d:/matcom45/samples/emd/emd.m"
_           lmin = fliplr(indmin(colon(1.0,1.0,min(length(indmin),NBSYM))));
            #line 119 "d:/matcom45/samples/emd/emd.m"
_           lsym = indmax(1.0);
            #line 120 "d:/matcom45/samples/emd/emd.m"
          } else {
            
            #line 121 "d:/matcom45/samples/emd/emd.m"
_           lmax = fliplr(indmax(colon(1.0,1.0,min(length(indmax),NBSYM))));
            #line 122 "d:/matcom45/samples/emd/emd.m"
_           lmin = (BR(fliplr(indmin(colon(1.0,1.0,min(length(indmin),NBSYM-1.0))))),1.0);
            #line 123 "d:/matcom45/samples/emd/emd.m"
_           lsym = 1.0;
            #line 124 "d:/matcom45/samples/emd/emd.m"
          }
          #line 125 "d:/matcom45/samples/emd/emd.m"
        } else {
          
          
          #line 127 "d:/matcom45/samples/emd/emd.m"
_         if (istrue(m(1.0)<m(indmax(1.0)))) {
            #line 128 "d:/matcom45/samples/emd/emd.m"
_           lmax = fliplr(indmax(colon(1.0,1.0,min(length(indmax),NBSYM))));
            #line 129 "d:/matcom45/samples/emd/emd.m"
_           lmin = fliplr(indmin(colon(2.0,1.0,min(length(indmin),NBSYM+1.0))));
            #line 130 "d:/matcom45/samples/emd/emd.m"
_           lsym = indmin(1.0);
            #line 131 "d:/matcom45/samples/emd/emd.m"
          } else {
            
            #line 132 "d:/matcom45/samples/emd/emd.m"
_           lmax = (BR(fliplr(indmax(colon(1.0,1.0,min(length(indmax),NBSYM-1.0))))),1.0);
            #line 133 "d:/matcom45/samples/emd/emd.m"
_           lmin = fliplr(indmin(colon(1.0,1.0,min(length(indmin),NBSYM))));
            #line 134 "d:/matcom45/samples/emd/emd.m"
_           lsym = 1.0;
            #line 135 "d:/matcom45/samples/emd/emd.m"
          }
          #line 136 "d:/matcom45/samples/emd/emd.m"
        }
        
        #line 138 "d:/matcom45/samples/emd/emd.m"
_       if (istrue(indmax(indmax.size())<indmin(indmin.size()))) {
          #line 139 "d:/matcom45/samples/emd/emd.m"
_         if (istrue(m(m.size())<m(indmax(indmax.size())))) {
            #line 140 "d:/matcom45/samples/emd/emd.m"
_           rmax = fliplr(indmax(colon(max(length(indmax)-NBSYM+1.0,1.0),1.0,indmax.size())));
            #line 141 "d:/matcom45/samples/emd/emd.m"
_           rmin = fliplr(indmin(colon(max(length(indmin)-NBSYM,1.0),1.0,indmin.size()-1.0)));
            #line 142 "d:/matcom45/samples/emd/emd.m"
_           rsym = indmin(indmin.size());
            #line 143 "d:/matcom45/samples/emd/emd.m"
          } else {
            
            #line 144 "d:/matcom45/samples/emd/emd.m"
_           rmax = (BR(lx),fliplr(indmax(colon(max(length(indmax)-NBSYM+2.0,1.0),1.0,indmax.size()))));
            #line 145 "d:/matcom45/samples/emd/emd.m"
_           rmin = fliplr(indmin(colon(max(length(indmin)-NBSYM+1.0,1.0),1.0,indmin.size())));
            #line 146 "d:/matcom45/samples/emd/emd.m"
_           rsym = lx;
            #line 147 "d:/matcom45/samples/emd/emd.m"
          }
          #line 148 "d:/matcom45/samples/emd/emd.m"
        } else {
          
          #line 149 "d:/matcom45/samples/emd/emd.m"
_         if (istrue(m(m.size())>m(indmin(indmin.size())))) {
            #line 150 "d:/matcom45/samples/emd/emd.m"
_           rmax = fliplr(indmax(colon(max(length(indmax)-NBSYM,1.0),1.0,indmax.size()-1.0)));
            #line 151 "d:/matcom45/samples/emd/emd.m"
_           rmin = fliplr(indmin(colon(max(length(indmin)-NBSYM+1.0,1.0),1.0,indmin.size())));
            #line 152 "d:/matcom45/samples/emd/emd.m"
_           rsym = indmax(indmax.size());
            #line 153 "d:/matcom45/samples/emd/emd.m"
          } else {
            
            #line 154 "d:/matcom45/samples/emd/emd.m"
_           rmax = fliplr(indmax(colon(max(length(indmax)-NBSYM+1.0,1.0),1.0,indmax.size())));
            #line 155 "d:/matcom45/samples/emd/emd.m"
_           rmin = (BR(lx),fliplr(indmin(colon(max(length(indmin)-NBSYM+2.0,1.0),1.0,indmin.size()))));