emd.cpp

vc代用matcom编程

            #line 156 "d:/matcom45/samples/emd/emd.m"
_           rsym = lx;
            #line 157 "d:/matcom45/samples/emd/emd.m"
          }
          #line 158 "d:/matcom45/samples/emd/emd.m"
        }
        
        #line 160 "d:/matcom45/samples/emd/emd.m"
_       tlmin = 2.0*t(lsym)-t(lmin);
        #line 161 "d:/matcom45/samples/emd/emd.m"
_       tlmax = 2.0*t(lsym)-t(lmax);
        #line 162 "d:/matcom45/samples/emd/emd.m"
_       trmin = 2.0*t(rsym)-t(rmin);
        #line 163 "d:/matcom45/samples/emd/emd.m"
_       trmax = 2.0*t(rsym)-t(rmax);
        
        #line 165 "d:/matcom45/samples/emd/emd.m"
        // in case symmetrized parts do not extend enough
        #line 166 "d:/matcom45/samples/emd/emd.m"
_       if (istrue(tlmin(1.0)>t(1.0))||istrue(tlmax(1.0)>t(1.0))) {
          #line 167 "d:/matcom45/samples/emd/emd.m"
_         if (istrue(lsym==indmax(1.0))) {
            #line 168 "d:/matcom45/samples/emd/emd.m"
_           lmax = fliplr(indmax(colon(1.0,1.0,min(length(indmax),NBSYM))));
            #line 169 "d:/matcom45/samples/emd/emd.m"
          } else {
            
            #line 170 "d:/matcom45/samples/emd/emd.m"
_           lmin = fliplr(indmin(colon(1.0,1.0,min(length(indmin),NBSYM))));
            #line 171 "d:/matcom45/samples/emd/emd.m"
          }
          #line 172 "d:/matcom45/samples/emd/emd.m"
_         if (istrue(lsym==1.0)) {
            #line 173 "d:/matcom45/samples/emd/emd.m"
_           display( error(TM("bug")) );
            #line 174 "d:/matcom45/samples/emd/emd.m"
          }
          #line 175 "d:/matcom45/samples/emd/emd.m"
_         lsym = 1.0;
          #line 176 "d:/matcom45/samples/emd/emd.m"
_         tlmin = 2.0*t(lsym)-t(lmin);
          #line 177 "d:/matcom45/samples/emd/emd.m"
_         tlmax = 2.0*t(lsym)-t(lmax);
          #line 178 "d:/matcom45/samples/emd/emd.m"
        }
        
        #line 180 "d:/matcom45/samples/emd/emd.m"
_       if (istrue(trmin(trmin.size())<t(lx))||istrue(trmax(trmax.size())<t(lx))) {
          #line 181 "d:/matcom45/samples/emd/emd.m"
_         if (istrue(rsym==indmax(indmax.size()))) {
            #line 182 "d:/matcom45/samples/emd/emd.m"
_           rmax = fliplr(indmax(colon(max(length(indmax)-NBSYM+1.0,1.0),1.0,indmax.size())));
            #line 183 "d:/matcom45/samples/emd/emd.m"
          } else {
            
            #line 184 "d:/matcom45/samples/emd/emd.m"
_           rmin = fliplr(indmin(colon(max(length(indmin)-NBSYM+1.0,1.0),1.0,indmin.size())));
            #line 185 "d:/matcom45/samples/emd/emd.m"
          }
          #line 186 "d:/matcom45/samples/emd/emd.m"
_         if (istrue(rsym==lx)) {
            #line 187 "d:/matcom45/samples/emd/emd.m"
_           display( error(TM("bug")) );
            #line 188 "d:/matcom45/samples/emd/emd.m"
          }
          #line 189 "d:/matcom45/samples/emd/emd.m"
_         rsym = lx;
          #line 190 "d:/matcom45/samples/emd/emd.m"
_         trmin = 2.0*t(rsym)-t(rmin);
          #line 191 "d:/matcom45/samples/emd/emd.m"
_         trmax = 2.0*t(rsym)-t(rmax);
          #line 192 "d:/matcom45/samples/emd/emd.m"
        }
        
        #line 194 "d:/matcom45/samples/emd/emd.m"
_       mlmax = m(lmax);
        
        #line 195 "d:/matcom45/samples/emd/emd.m"
_       mlmin = m(lmin);
        #line 196 "d:/matcom45/samples/emd/emd.m"
_       mrmax = m(rmax);
        
        #line 197 "d:/matcom45/samples/emd/emd.m"
_       mrmin = m(rmin);
        
        #line 199 "d:/matcom45/samples/emd/emd.m"
        // definition of envelopes from interpolation
        
        #line 201 "d:/matcom45/samples/emd/emd.m"
_       envmax = interp1((BR(tlmax),t(indmax),trmax),(BR(mlmax),m(indmax),mrmax),t,TM("spline"));
        #line 202 "d:/matcom45/samples/emd/emd.m"
_       envmin = interp1((BR(tlmin),t(indmin),trmin),(BR(mlmin),m(indmin),mrmin),t,TM("spline"));
        
        #line 204 "d:/matcom45/samples/emd/emd.m"
_       envmoy = (envmax+envmin)/2.0;
        
        #line 206 "d:/matcom45/samples/emd/emd.m"
_       m = m-envmoy;
        
        #line 208 "d:/matcom45/samples/emd/emd.m"
_       /*[indmin,indmax,indzer] = */extr(m,i_o,indmin,indmax,indzer);
        #line 209 "d:/matcom45/samples/emd/emd.m"
_       lm = length(indmin);
        #line 210 "d:/matcom45/samples/emd/emd.m"
_       lM = length(indmax);
        #line 211 "d:/matcom45/samples/emd/emd.m"
_       nem = lm+lM;
        #line 212 "d:/matcom45/samples/emd/emd.m"
_       nzm = length(indzer);
        
        #line 214 "d:/matcom45/samples/emd/emd.m"
        // evaluation of mean zero
        #line 215 "d:/matcom45/samples/emd/emd.m"
_       sx = rdivide(2.0*(abs(envmoy)),(abs(envmax-envmin)));
        #line 216 "d:/matcom45/samples/emd/emd.m"
_       s = mean(sx);
        
        #line 218 "d:/matcom45/samples/emd/emd.m"
        // display
        
        #line 220 "d:/matcom45/samples/emd/emd.m"
_       if (istrue(tst)) {
          #line 221 "d:/matcom45/samples/emd/emd.m"
_         display( subplot(4.0,1.0,1.0) );
          #line 222 "d:/matcom45/samples/emd/emd.m"
_         plot((CL(t),mp));
          #line 222 "d:/matcom45/samples/emd/emd.m"
_         hold(TM("on"));
          #line 223 "d:/matcom45/samples/emd/emd.m"
_         plot((CL(t),envmax,TM("--k")));
          #line 223 "d:/matcom45/samples/emd/emd.m"
_         plot((CL(t),envmin,TM("--k")));
          #line 223 "d:/matcom45/samples/emd/emd.m"
_         plot((CL(t),envmoy,TM("r")));
          
          #line 225 "d:/matcom45/samples/emd/emd.m"
_         title((CL((BR(TM("IMF ")),int2str(k),TM(";   iteration "),int2str(nbit),TM(" before sifting")))));
          #line 226 "d:/matcom45/samples/emd/emd.m"
_         display( set(gca(),(CL(TM("XTick")),nop_M)) );
          #line 227 "d:/matcom45/samples/emd/emd.m"
_         hold(TM("off"));
          
          #line 229 "d:/matcom45/samples/emd/emd.m"
_         display( subplot(4.0,1.0,2.0) );
          #line 230 "d:/matcom45/samples/emd/emd.m"
_         display( plot((CL(t),sx)) );
          #line 231 "d:/matcom45/samples/emd/emd.m"
_         hold(TM("on"));
          #line 232 "d:/matcom45/samples/emd/emd.m"
_         display( plot((CL(t),sdt,TM("--r"))) );
          #line 233 "d:/matcom45/samples/emd/emd.m"
_         display( plot((CL(t),sd2t,TM(":k"))) );
          #line 234 "d:/matcom45/samples/emd/emd.m"
_         display( title((CL(TM("stop parameter")))) );
          #line 235 "d:/matcom45/samples/emd/emd.m"
_         display( set(gca(),(CL(TM("XTick")),nop_M)) );
          #line 236 "d:/matcom45/samples/emd/emd.m"
_         hold(TM("off"));
          
          #line 238 "d:/matcom45/samples/emd/emd.m"
_         display( subplot(4.0,1.0,3.0) );
          #line 239 "d:/matcom45/samples/emd/emd.m"
_         display( plot((CL(t),m)) );
          #line 240 "d:/matcom45/samples/emd/emd.m"
_         title((CL((BR(TM("IMF ")),int2str(k),TM(";   iteration "),int2str(nbit),TM(" after sifting")))));
          #line 241 "d:/matcom45/samples/emd/emd.m"
_         display( set(gca(),(CL(TM("XTick")),nop_M)) );
          
          #line 243 "d:/matcom45/samples/emd/emd.m"
_         subplot(4.0,1.0,4.0);
          #line 244 "d:/matcom45/samples/emd/emd.m"
_         display( plot((CL(t),r-m)) );
          #line 245 "d:/matcom45/samples/emd/emd.m"
_         title((CL(TM("residue"))));
          #line 246 "d:/matcom45/samples/emd/emd.m"
_         display( disp((BR(TM("stop parameter mean value : ")),num2str(s))) );
          #line 247 "d:/matcom45/samples/emd/emd.m"
_         if (istrue(tst==2.0)) {
            #line 248 "d:/matcom45/samples/emd/emd.m"
_           display( pauseM(0.01) );
            #line 249 "d:/matcom45/samples/emd/emd.m"
          } else {
            
            #line 250 "d:/matcom45/samples/emd/emd.m"
_           display(pauseM());
            #line 251 "d:/matcom45/samples/emd/emd.m"
          }
          
          #line 253 "d:/matcom45/samples/emd/emd.m"
        }
        
        #line 255 "d:/matcom45/samples/emd/emd.m"
        // end loop : stops if not enough extrema
        #line 256 "d:/matcom45/samples/emd/emd.m"
_       if (istrue(nem<3.0)) {
          #line 257 "d:/matcom45/samples/emd/emd.m"
_         test = 1.0;
          #line 258 "d:/matcom45/samples/emd/emd.m"
        }
        
        #line 260 "d:/matcom45/samples/emd/emd.m"
_       mp = m;
        #line 261 "d:/matcom45/samples/emd/emd.m"
_       nbit = nbit+1.0;
        #line 262 "d:/matcom45/samples/emd/emd.m"
_       NbIt = NbIt+1.0;
        
        #line 264 "d:/matcom45/samples/emd/emd.m"
_       if (istrue(nbit==(MAXITERATIONS-1.0))) {
          #line 265 "d:/matcom45/samples/emd/emd.m"
_         display( warning((BR(TM("forced stop of sifting : too many iterations... mode ")),int2str(k),TM(". stop parameter mean value : ") \
            ,num2str(s))) );
          #line 266 "d:/matcom45/samples/emd/emd.m"
        }
        
        #line 268 "d:/matcom45/samples/emd/emd.m"
      }
      #line 269 "d:/matcom45/samples/emd/emd.m"
_     imf(k,c_p) = m;
      #line 270 "d:/matcom45/samples/emd/emd.m"
_     if (istrue(tst)) {
        #line 271 "d:/matcom45/samples/emd/emd.m"
_       display( disp((BR(TM("mode ")),int2str(k),TM(" enregistre"))) );
        #line 272 "d:/matcom45/samples/emd/emd.m"
      }
      #line 273 "d:/matcom45/samples/emd/emd.m"
_     nbits(k) = nbit;
      #line 274 "d:/matcom45/samples/emd/emd.m"
_     k = k+1.0;
      #line 275 "d:/matcom45/samples/emd/emd.m"
_     r = r-m;
      #line 276 "d:/matcom45/samples/emd/emd.m"
_     /*[indmin,indmax,indzer] = */extr(r,i_o,indmin,indmax,indzer);
      #line 277 "d:/matcom45/samples/emd/emd.m"
_     ner = length(indmin)+length(indmax);
      #line 278 "d:/matcom45/samples/emd/emd.m"
_     nzr = length(indzer);
      #line 279 "d:/matcom45/samples/emd/emd.m"
_     nbit = 1.0;
      
      #line 281 "d:/matcom45/samples/emd/emd.m"
_     if (istrue((max(r)-min(r))<(1e-10)*(max(x)-min(x)))) {
        #line 282 "d:/matcom45/samples/emd/emd.m"
_       if (istrue(ner>2.0)) {
          #line 283 "d:/matcom45/samples/emd/emd.m"
_         display( warning(TM("forced stop of EMD : too small amplitude")) );
          #line 284 "d:/matcom45/samples/emd/emd.m"
        } else {
          
          
          #line 286 "d:/matcom45/samples/emd/emd.m"
_         display( disp(TM("forced stop of EMD : too small amplitude")) );
          #line 287 "d:/matcom45/samples/emd/emd.m"
        }
        #line 288 "d:/matcom45/samples/emd/emd.m"
_       break;
        #line 289 "d:/matcom45/samples/emd/emd.m"
      }
      
      #line 291 "d:/matcom45/samples/emd/emd.m"
    }
    
    #line 293 "d:/matcom45/samples/emd/emd.m"
_   imf(k,c_p) = r;
    
    #line 295 "d:/matcom45/samples/emd/emd.m"
_   ort = io(x,imf);
    
    #line 297 "d:/matcom45/samples/emd/emd.m"
_   if (istrue(tst)) {
      #line 298 "d:/matcom45/samples/emd/emd.m"
_     display(figure_close());
      #line 299 "d:/matcom45/samples/emd/emd.m"
    }
    
    call_stack_end;
    
    // nargin, nargout exit code
    nargin_val=old_nargin; nargout_val=old_nargout;
    
    // function exit code
    x.setname(NULL); 
    imf__o=imf; ort__o=ort; nbits__o=nbits; 
    return x_M;
    end_scope
  }
  
  
  Mm emd(Mm x) {
    begin_scope
    double old_nargin=nargin_val; nargin_val=1.0; nargin_set=1;
    double old_nargout=nargout_val; nargout_val=1.0; nargout_set=1;
    dMm(imf__o); dMm(ort__o); dMm(nbits__o); 
    emd(x, i_o, imf__o, ort__o, nbits__o);
    nargout_val=old_nargout;
    nargin_val=old_nargin;
    return(imf__o);
    end_scope
  }
  
  Mm emd(Mm x, i_o_t, Mm& imf__o, Mm& ort__o) {
    begin_scope
    double old_nargin=nargin_val; nargin_val=1.0; nargin_set=1;
    double old_nargout=nargout_val; nargout_val=2.0; nargout_set=1;
    dMm(nbits__o); 
    emd(x, i_o, imf__o, ort__o, nbits__o);
    nargout_val=old_nargout;
    nargin_val=old_nargin;
    return x_M;
    end_scope
  }