static mrs_real Std (const realvec& beatHistogram)
{
  return sqrt (beatHistogram.var ());
}

void 
SOM::myProcess(realvec& in, realvec& out)
{
	mrs_string mode = getctrl("mrs_string/mode")->to<mrs_string>();

	mrs_natural o,t;
	mrs_real geom_dist;
	mrs_real geom_dist_gauss;
  
	int px;
	int py;

	MarControlAccessor acc_grid(ctrl_gridmap_);
	realvec& grid_map = acc_grid.to<mrs_realvec>();  


	if (mode == "train")  
	{
		
		mrs_real dx;
		mrs_real dy;
		mrs_real adj;

		for (t=0; t < inSamples_; t++) 
		{

			
			px = (int) in(inObservations_-2, t);
			py = (int) in(inObservations_-1, t);


			
			if ((px == -1.0)&&(px == -1.0))
			{
				find_grid_location(in, t);
				px = (int) grid_pos_(0);
				py = (int) grid_pos_(1);				
			}
			out(0,t) = px;
			out(1,t) = py;
			out(2,t) = in(inObservations_-3,t);	  
			
			for (int x=0; x < grid_width_; x++) 
				for (int y=0; y < grid_height_; y++)
				{
					dx = px-x;
					dy = py-y;
					geom_dist = sqrt((double)(dx*dx + dy*dy));
					geom_dist_gauss = gaussian( geom_dist, 0.0, neigh_std_, false);
					
					// subtract map vector from training data vector 
					adj = alpha_ * geom_dist_gauss;
					for (o=0; o < inObservations_-3; o++) 
					{
						adjustments_(o) = in(o,t) - grid_map(x * grid_height_ + y, o);
						adjustments_(o) *= adj;
						grid_map(x * grid_height_ + y, o) += adjustments_(o);
					}
				}
		}
		
		alpha_ *= getctrl("mrs_real/alpha_decay_train")->to<mrs_real>();
		neigh_std_ *= getctrl("mrs_real/neighbourhood_decay_train")->to<mrs_real>();	  

		

	}
	if (mode == "init")
	{
		mrs_real dx;
		mrs_real dy;
		mrs_real adj;

		mrs_real std_ = getctrl("mrs_real/std_factor_init")->to<mrs_real>();
		neigh_std_ = ((0.5*(grid_width_+grid_height_)) * std_);

		for (t=0; t < inSamples_; t++) 
		{
			// no need to find grid locations, just read from the file
			px = (int) in( in.getRows() - 2, t);
			py = (int) in( in.getRows() - 1, t);
			for(int i =0; i < inObservations_ - 3; ++i)
			{
				grid_map(px * grid_height_ + py, i) = in(i);
			}	

			for (int x=0; x < grid_width_; x++) 
				for (int y=0; y < grid_height_; y++)
				{
					dx = px-x;
					dy = py-y;
					geom_dist = sqrt((double)(dx*dx + dy*dy));
					geom_dist_gauss = gaussian( geom_dist, 0.0, neigh_std_, false);

					// subtract map vector from training data vector 
					adj = alpha_ * geom_dist_gauss;
					for (o=0; o < inObservations_-3; o++) 
					{
						adjustments_(o) = in(o,t) - grid_map(x * grid_height_ + y, o);
						adjustments_(o) *= adj;
//.........这里部分代码省略.........

void 
GMMClassifier::myProcess(realvec& in, realvec& out)
{
	mrs_string mode = ctrl_mode_->to<mrs_string>();
	
	// reset 
	if ((prev_mode_ == "predict") && (mode == "train"))
	{
		//just drop all accumulated feature vectors and
		//copy take the new ones from the input
		trainMatrix_ = in;
	}
	
	if (mode == "train")  
	{
		MRSASSERT(trainMatrix_.getRows() == inObservations_);
		
		//stretch to acommodate input feat Vecs
		mrs_natural storedFeatVecs = trainMatrix_.getCols();
		trainMatrix_.stretch(inObservations_, storedFeatVecs + inSamples_);
		
		//append input data
		for(mrs_natural c=0; c < inSamples_; ++c)
			for(mrs_natural r = 0; r < inObservations_; ++r)
				trainMatrix_(r, c+storedFeatVecs) = in(r,c);
	}
	
	if (mode == "predict")
	{
		mrs_real maxProb = 0.0;
		mrs_natural maxClass = 0;
		mrs_real prob;
		mrs_real dist;
		realvec vec;
		realvec means;
		realvec covars;
		
		MRSASSERT(trainMatrix_.getRows() == inObservations_);
		
		for(mrs_natural t=0; t < inSamples_; ++t)
		{	
			
			in.getCol(t, vec);
			
			for (mrs_natural cl=0; cl < classSize_; cl++)
			{
				for (mrs_natural k=0; k < nMixtures_; k++)
				{
					means_[cl].getCol(k, means);
					covars_[cl].getCol(k, covars);
					dist = NumericLib::mahalanobisDistance(vec, means, covars);					
					likelihoods_(cl,k) = weights_[cl](k) / dist;
				}
				prob = 0.0;
				for (mrs_natural k=0; k < nMixtures_; k++)
				{
					prob += likelihoods_(cl,k);
				}
				if (prob > maxProb) 
				{
					maxProb = prob;
					maxClass = cl;
				}
			}
			out(0,t) = in(labelRow_, t);
			out(1,t) = (mrs_real)maxClass; //FIXME: what about he maxProb (i.e. Confidence)?
		}
	}
	
	prev_mode_ = mode;
}

void
ADRessSpectrum::myProcess(realvec& in, realvec& out)
{
	mrs_natural t;
	out.setval(0.0);

	//output spectrum of the "selected" source, given d and H
	mrs_natural H = (mrs_natural)(beta_* ctrl_H_->to<mrs_natural>());
	if(H < 0)
	{
		H = 0;
		ctrl_H_->setValue(0.0);
	}
	if(H > beta_)
	{
		H = beta_;
		ctrl_H_->setValue(1.0);
	}

	mrs_natural H2 = H/2;
	
	mrs_natural d = (mrs_natural)(beta_*ctrl_d_->to<mrs_real>());
	if(d < 0)
	{
		d = 0;
		ctrl_d_->setValue(0.0);
	}
	if(d > beta_)
	{
		d = beta_;
		ctrl_d_->setValue(1.0);
	}

	mrs_real mag = 0;
	mrs_real phase = 0;
	mrs_real azim = 0;
	
	
	for(mrs_natural k=0; k < N2_; ++k)
	{
		//get magnitude, phase and azimuth of bin k from input
		mag = 0.0;
		for(t=0; t <= beta_; ++t)
		{
			//search for non-zero values in azimuth plane
			azim = -1;
			if(in(k,t+1) > 0.0)
			{
				azim = t;
				mag = in(k,t+1);
				phase = in(k, 0);
				break;
			}
			if(in(k+N2_,t+1) > 0.0)
			{
				azim = beta_*2-t;
				mag = in(k+N2_,t+1);
				phase = in(k+N2_, 0);
				break;
			}
		}

		if(azim < 0)
		{
			//no sound at this bin,
			//so do not send anything to output
			continue;
		}
		
		//check if bin is inside specified range,
		//otherwise, send nothing to output
		if(abs(d-azim) <= H2)
		{
			//convert back to rectangular form
			re_ = mag*cos(phase);
			im_ = mag*sin(phase);

			//write bin to output
			if (k==0)
			{
				out(0,0) = re_; //DC
			}
			else if (k == N2_-1) 
			{
				out(1, 0) = re_; //Nyquist
			}
			else
			{
				out(2*k, 0) = re_;  //all other bins
				out(2*k+1, 0) = im_;
			}
		}
	}
}

void 
PeakConvert::myProcess(realvec& in, realvec& out)
{
	mrs_natural o;
	mrs_real a, c;
	mrs_real b, d;
	mrs_real phasediff;

	out.setval(0);
	peakView pkViewOut(out);

	for(mrs_natural f=0 ; f < inSamples_; ++f)
	{
		//we should avoid the first empty frames, 
		//that will contain silence and consequently create 
		//discontinuities in the signal, ruining the peak calculation!
		//only process if we have a full data vector (i.e. no zeros)
		if(frame_ >= skip_) 
		{
			// handle amplitudes from shifted spectrums at input
			for (o=0; o < size_; o++)
			{
				if (o==0) //DC bins
				{
					a = in(0,f); 
					b = 0.0; 
					c = in(N_, f);
					d = 0.0;
				}
				else if (o == size_-1) //Nyquist bins
				{
					a = in(1, f);
					b = 0.0;
					c = in(N_+1, f);
					d = 0.0;
				}
				else //all other bins
				{
					a = in(2*o, f);
					b = in(2*o+1, f);
					c = in(N_+2*o, f);
					d = in(N_+2*o+1, f);
				}

				// computer magnitude value 
				//mrs_real par = lobe_value_compute (0, 1, 2048); //[?]

				// compute phase
				phase_(o) = atan2(b,a);

				// compute precise frequency using the phase difference
				lastphase_(o)= atan2(d,c);
				if(phase_(o) >= lastphase_(o))
					phasediff = phase_(o)-lastphase_(o);
				else
					phasediff = phase_(o)-lastphase_(o)+TWOPI;
				if(prec_)
					frequency_(o) = phasediff * factor_ ;
				else
					frequency_(o) = o*fundamental_;

				// compute precise amplitude
				mag_(o) = sqrt((a*a + b*b))*2; //*4/0.884624;//*50/3); // [!]
				mrs_real mag = lobe_value_compute((o * fundamental_-frequency_(o))/factor_, 1, N_);
				magCorr_(o) = mag_(o)/mag;

				// computing precise frequency using the derivative method // use at your own risk	[?]
				/*	mrs_real lastmag = sqrt(c*c + d*d);
				mrs_real rap = (mag_(o)-lastmag)/(lastmag*2);
				f=asin(rap);
				f *= (getctrl("mrs_real/israte")->to<mrs_real>()*inObservations/2.0)/PI;
				*/
				// rough frequency and amplitude
				//frequency_(o) = o * fundamental_;
				//magCorr_(o) = mag_(o);

				if(lastfrequency_(o) != 0.0)
					deltafrequency_(o) = (frequency_(o)-lastfrequency_(o))/(frequency_(o)+lastfrequency_(o));

				deltamag_(o) = (mag_(o)-lastmag_(o))/(mag_(o)+lastmag_(o));

				// remove potential peak if frequency too irrelevant
				if(abs(frequency_(o)-o*fundamental_) > 0.5*fundamental_)
					frequency_(o)=0.0;

				lastfrequency_(o) = frequency_(o);
				lastmag_(o) = mag_(o);
			}

			// select bins with local maxima in magnitude (--> peaks)
			realvec peaks_ = mag_;
			realvec tmp_;
			peaker_->updControl("mrs_real/peakStrength", 0.2);// to be set as a control [!]
			peaker_->updControl("mrs_natural/peakStart", downFrequency_);   // 0
			peaker_->updControl("mrs_natural/peakEnd", upFrequency_);  // size_
			peaker_->updControl("mrs_natural/inSamples", mag_.getCols());
			peaker_->updControl("mrs_natural/inObservations", mag_.getRows());
			peaker_->updControl("mrs_natural/onSamples", peaks_.getCols());
			peaker_->updControl("mrs_natural/onObservations", peaks_.getRows());
			if(pick_)
//.........这里部分代码省略.........

void 
BeatHistogram::myProcess(realvec& in, realvec& out)
{
  
  if (reset_)
  {
    out.setval(0.0);
    reset_ = false;
    setctrl("mrs_bool/reset", false);
  }
  

  //out.setval(0.0);
  
  mrs_natural bin=0;
  mrs_real amp;
  mrs_real srate = getctrl("mrs_real/israte")->to<mrs_real>();
  mrs_natural count = 1;
  mrs_natural prev_bin =endBin_-1;
  mrs_natural pprev_bin =endBin_-1;
  mrs_real sumamp = 0.0;
  mrs_real tempo_weight = 0.0;


#ifdef MARSYAS_MATLAB
#ifdef MTLB_DBG_LOG
	MATLAB_PUT(in, "acr");
	MATLAB_EVAL("figure(1);plot(acr),grid on");
#endif
#endif

	
	for (mrs_natural o=0; o < inObservations_; o++)
	{
	  for (mrs_natural t = 1; t < inSamples_; t++)
	  {
		bin = (mrs_natural)(( (2*srate) * 60.0 * factor_ / (t+1)) + 0.5);
		
		amp = in(o,t);


		// optional tempo weight 
		
		if (getctrl("mrs_bool/tempoWeighting")->to<mrs_bool>())
		  {
		     tempo_weight = 5.0 * 
		       log10((t+1) * 400.0 / (srate * 60.0 * factor_)) * 
		       log10((t+1) * 400.0 / (srate * 60.0 * factor_));
		     tempo_weight = exp(0.5 * tempo_weight * tempo_weight);
		  }
		else 
		  {
		    tempo_weight = 1.0;
		  }
		
		amp = tempo_weight * amp;
		
		if (amp < 0.0) 
		  amp = 0.0;

		
		if ((bin > 40)&&(bin < endBin_))
		  {
		  if (prev_bin == bin) 
		  {
			sumamp += amp;
			count++;
		  }
		  else 
		  {
			sumamp += amp;
			out(o,prev_bin) += ((sumamp / count));
			count = 1;
			sumamp = 0.0;
		  }

		  // linear interpolation of the "not-set" bins...
		  if (pprev_bin-prev_bin > 1)
		  {
            // prev is x0, pprev is x1
            mrs_real x0 = prev_bin;
            mrs_real x1 = pprev_bin;
            mrs_real y0 = out(o, prev_bin);
            mrs_real y1 = out(o, pprev_bin);
			for (mrs_natural k = prev_bin+1; k < pprev_bin; k++)
			  out (o,k)	= y0 + (y1-y0)*(k-x0)/(x1-x0);
		  }

		  pprev_bin = prev_bin;
		  prev_bin = bin;
		}
	  }
	  

	  
	}
	

#ifdef MARSYAS_MATLAB
#ifdef MTLB_DBG_LOG
//.........这里部分代码省略.........

void
AutoCorrelation::myProcess(realvec& in, realvec& out)
{





  mrs_natural t,o;
  k_ = ctrl_magcompress_->to<mrs_real>();

  // Copy to output to perform inplace fft and zeropad to double size

  scratch_.create(fftSize_); //scratch_ needs to be reset every time
  for (o=0; o < inObservations_; o++)
  {
    for (t=lowSamples_; t < (lowSamples_+numSamples_); t++)
    {
      scratch_(t-(lowSamples_)) = in(o,t);
    }





    //zeropad
    for(t=(lowSamples_+numSamples_); t < fftSize_; t++)
      scratch_(t) = 0.0;


    //get pointer to data (THIS BREAKS ENCAPSULATION! FIXME [!])
    mrs_real *tmp = scratch_.getData();

    //compute forward FFT (of size fftSize_)
    myfft_->rfft(tmp, fftSize_/2, FFT_FORWARD); //rfft() takes as second argument half of the desired FFT size (see fft.cpp)

    // Special case for zero and Nyquist/2,
    // which only have real part
    if (k_ == 2.0)
    {
      re_ = tmp[0];
      tmp[0] = re_ * re_;
      re_ = tmp[1];
      tmp[1] = re_ * re_;
    }
    else
    {
      re_ = tmp[0];
      re_ = sqrt(re_ * re_);
      tmp[0] = pow(re_, k_);
      re_ = tmp[1];
      re_ = sqrt(re_ * re_);
      tmp[1] = pow(re_, k_);
    }

    // Compress the magnitude spectrum and zero
    // the imaginary part.
    for (t=1; t < fftSize_/2; t++)
    {
      re_ = tmp[2*t];
      im_ = tmp[2*t+1];
      if (k_ == 2.0)
        am_ = re_ * re_ + im_ * im_;
      else
      {
        am_ = sqrt(re_ * re_ + im_ * im_);
        am_ = pow(am_, k_);
      }
      tmp[2*t] = am_;
      tmp[2*t+1] = 0;
    }

    // Take the inverse Fourier Transform (of size fftSize_)
    myfft_->rfft(tmp, fftSize_/2, FFT_INVERSE);

    // Copy result to output
    if(normalize_)
    {
      cout << "NORM Normalization happening" << endl;
      for (t=0; t < onSamples_; t++)
      {
        out(o,t) = scratch_(t)*norm_(t);
      }
    }
    else
      for (t=0; t < onSamples_; t++)
      {
        // out(o,t) = 0.1 * scratch_(t) + 0.99 * out(o,t);
        out(o,t) = 1.0 * scratch_(t) + 0.0 * out(o,t);
        // out(o,t) = 0.5 * scratch_(t) + 0.5 * out(o,t);
        // out(o,t) +=  scratch_(t);

      }

  }


  if (ctrl_makePositive_->to<mrs_bool>())
  {
    out -= out.minval();
//.........这里部分代码省略.........

mrs_real
McAulayQuatieri::peakTrack(realvec& vec, mrs_natural frame, mrs_natural grpOne, mrs_natural grpTwo)
{
  mrs_real dist;
  mrs_natural candidate;
  mrs_natural lastMatched = -1;
  mrs_natural matchedTracks = 0;

  mrs_real delta = ctrl_delta_->to<mrs_real>();

  if(frame+1 >= vec.getCols())
  {
    MRSERR("McAulayQuatieri::peakTrack - frame index is bigger than the input vector!");
    return -1.0;
  }

  peakView tmpPeakView(vec);

  //get the trackID for any future track to be born (in STEP 3 - see below)
  mrs_natural nextTrack = tmpPeakView.getFrameNumPeaks(0, grpOne);

  //iterate over peaks in current frame
  for(mrs_natural n = 0; n < tmpPeakView.getFrameNumPeaks(frame, grpOne); ++n)
  {
    mrs_real lastdist = MAXREAL;
    candidate = -1;

    // STEP 1
    // find a candidate match on the next frame for each peak (i.e. track) in current frame
    for(mrs_natural m = lastMatched + 1; m < tmpPeakView.getFrameNumPeaks(frame+1, grpTwo); ++m)
    {
      //set track parameter of all peaks of next frame to -1 so we know later
      //which ones were not matched (=> BIRTH of new tracks)
      tmpPeakView(m, peakView::pkTrack, frame+1, grpTwo) = -1.0;

      dist = abs(tmpPeakView(n, peakView::pkFrequency, frame, grpOne) - tmpPeakView(m, peakView::pkFrequency, frame+1, grpTwo));
      if (dist < delta && dist < lastdist)
      {
        //found a candidate!
        lastdist  = dist;
        candidate = m;
      }
    }

    // STEP 2
    // must confirm candidate (if any)
    if(candidate >= 0) //check if a candidate was found
    {
      //confirm if this is not the last peak in current frame
      if(n < tmpPeakView.getFrameNumPeaks(frame, grpOne)-1)
      {
        //check the next remaining peak in current frame and see if it is a better match for the found candidate
        dist = abs(tmpPeakView(n+1, peakView::pkFrequency, frame, grpOne) - tmpPeakView(candidate, peakView::pkFrequency, frame+1, grpTwo));
        if(dist < lastdist)
        {
          // it is a better match! Check two additional conditions:
          // 1. an unmatched lower freq candidate should exist
          // 2. it is inside the frequency interval specified by delta
          if(candidate - 1 > lastMatched)
          {
            if(abs(tmpPeakView(n, peakView::pkFrequency, frame, grpOne) - tmpPeakView(candidate-1, peakView::pkFrequency, frame+1, grpTwo)) < delta)
            {
              //found a peak to continue the track -> confirm candidate!
              tmpPeakView(candidate-1, peakView::pkTrack, frame+1, grpTwo) = tmpPeakView(n, peakView::pkTrack, frame, grpOne);
              lastMatched = candidate-1;
              matchedTracks++;
            }
          }
        }
        else
        {
          //no better match than this one, so confirm candidate!
          tmpPeakView(candidate, peakView::pkTrack, frame+1, grpTwo) = tmpPeakView(n, peakView::pkTrack, frame, grpOne);
          lastMatched = candidate;
          matchedTracks++;
        }
      }
      else
      {
        //if this was the last peak in current frame, so inherently it was the best match.
        //Candidate is therefore automatically confirmed and can be propagated.
        tmpPeakView(candidate, peakView::pkTrack, frame+1, grpTwo) = tmpPeakView(n, peakView::pkTrack, frame, grpOne);
        lastMatched = candidate;
        matchedTracks++;
      }
    }
  } //end of loop on peaks of current frame

  // STEP 3
  // check for any unmatched peaks in the next frame and give BIRTH to new tracks!
  for(mrs_natural m = 0; m < tmpPeakView.getFrameNumPeaks(frame+1, grpTwo); ++m)
  {
    if(tmpPeakView(m, peakView::pkTrack, frame+1, grpTwo) == -1.0)
      tmpPeakView(m, peakView::pkTrack, frame+1, grpTwo) = nextTrack++; //BIRTH of new track
  }

  return matchedTracks;
}

void
Filter::myProcess(realvec& in, realvec& out)
{
  //checkFlow(in,out);

  mrs_natural i,j,c;
  mrs_natural size = in.getCols();
  mrs_natural stateSize = state_.getCols();
  mrs_natural channels = in.getRows();

  mrs_real gain = getctrl("mrs_real/fgain")->to<mrs_real>();

  // State array holds the various delays for the difference equation
  // of the filter. Similar implementation as described in the manual
  // for MATLAB Signal Processing Toolbox. State corresponds to
  // the z, num_coefs to the b and denom_coefs to the a vector respectively
  // in_window is the input x(n) and out_window is the output y(n)

  //dcoeffs_/=10;


  // state_.setval(0);
  if (norder_ == dorder_) {
    for (c = 0; c < channels; ++c) {
      for (i = 0; i < size; ++i) {
        out(c,i) = ncoeffs_(0) * in(c,i) + state_(c,0);
        for (j = 0; j < stateSize - 1; j++)
        {
          state_(c,j) = ncoeffs_(j+1) * in(c,i) + state_(c,j+1) - dcoeffs_(j+1) * out(c,i);
        }
        state_(c,stateSize - 1) = ncoeffs_(order_-1) * in(c,i) - dcoeffs_(order_-1) * out(c,i);
      }
    }
  }
  else if (norder_ < dorder_) {
    for (c = 0; c < channels; ++c) {
      for (i = 0; i < size; ++i) {
        out(c,i) = ncoeffs_(0) * in(c,i) + state_(c,0);
        for (j = 0; j < norder_ - 1; j++)
        {
          state_(c,j) = ncoeffs_(j+1) * in(c,i) + state_(c,j+1) - dcoeffs_(j+1) * out(c,i);
        }
        for (j = norder_ - 1; j < stateSize - 1; j++)
        {
          state_(c,j) = state_(c,j+1) - dcoeffs_(j+1) * out(c,i);
        }
        state_(c,stateSize - 1) = -dcoeffs_(order_ - 1) * out(c,i);
      }
    }
  }
  else {
    for (c = 0; c < channels; ++c) {
      for (i = 0; i < size; ++i) {
        out(c,i) = ncoeffs_(0) * in(c,i) + state_(c,0);
        for (j = 0; j < dorder_ - 1; j++)
        {
          state_(c,j) = ncoeffs_(j+1) * in(c,i) + state_(c,j+1) - dcoeffs_(j+1) * out(c,i);
        }
        for (j = dorder_ - 1; j < stateSize - 1; j++)
        {
          state_(c,j) = ncoeffs_(j+1) * in(c,i) + state_(c,j+1);
        }
        state_(c,stateSize - 1) = ncoeffs_(order_-1) * in(c,i);
      }
    }
  }
  out *= gain;


  //		MATLAB_PUT(in, "Filter_in");
  //	 	MATLAB_PUT(out, "Filter_out");
  //	 	MATLAB_PUT(ncoeffs_, "ncoeffs_");
  //	 	MATLAB_PUT(dcoeffs_, "dcoeffs_");
  //	 	MATLAB_EVAL("MAT_out = filter(ncoeffs_, dcoeffs_, Filter_in)");
  //
  //	 	MATLAB_EVAL("spec_in = abs(fft(Filter_in));");
  //	 	MATLAB_EVAL("spec_out = abs(fft(Filter_out));");
  //	 	MATLAB_EVAL("spec_mat = abs(fft(MAT_out));");
  //
  //	 	MATLAB_EVAL("subplot(2,1,1);plot(Filter_in);hold on; plot(Filter_out, 'r'); plot(MAT_out, 'g');hold off");
  //	 	MATLAB_EVAL("subplot(2,1,2);plot(spec_in(1:end/2));hold on; plot(spec_out(1:end/2),'r');plot(spec_mat(1:end/2),'g');hold off;");
  //	 	MATLAB_EVAL("h = abs(fft([1 -.97], length(Filter_in)));");
  //	 	MATLAB_EVAL("hold on; plot(h(1:end/2), 'k'); hold off");
  //		//MATLAB_GET("MAT_out", out)
  //
}

void
PeakerOnset::myProcess(realvec& in, realvec& out)
{
	mrs_natural o,t;
	(void) o;
	ctrl_onsetDetected_->setValue(false);
	ctrl_confidence_->setValue(0.0);
	out.setval(0.0);

	mrs_natural w = ctrl_lookAheadSamples_->to<mrs_natural>();

	if(w == 0)
		return;

	//point to check for an onset
	mrs_natural checkPoint = inSamples_-1-w;
	mrs_real checkPointValue = in(checkPoint);
	mrs_bool isOnset = true;

	//check first condition
	mrs_natural interval = mrs_natural(2.0/3.0*w);
	//for (t = inSamples_-1; t >= inSamples_-1-2*w ; t--)
	for(t=checkPoint-interval; t <= checkPoint+interval; t++)
	{
		if(checkPointValue < in(t))
		{
			isOnset = false;
			MRSDIAG("PeakerOnset::myProcess() - Failed 1st condition!");
			break;
		}
	}

// 	//new check proposed by Fabien Gouyon
// 	mrs_real ww = w/2;
// 	mrs_real maxVal = MINREAL;
// 	for(t = inSamples_-1-ww; t > checkPoint; t--)
// 	{
// 		if(in(t) > maxVal)
// 		{
// 			maxVal = in(t);
// 		}
// 		else
// 		{
// 			isOnset = false;
// 			//cout << "failed 1st condition!" << endl;
// 			break;
// 		}
// 	}
// 	maxVal = MINREAL;
// 	for(t = inSamples_-1-2*ww; t < checkPoint; t++)
// 	{
// 		if(in(t) > maxVal)
// 		{
// 			maxVal = in(t);
// 		}
// 		else
// 		{
// 			isOnset = false;
// 			//cout << "failed 1st condition!" << endl;
// 			break;
// 		}
// 	}

/* Last version (by lgmartins) -> corrected (below) for not being strict to the window size defined by the precede ShiftInput
	//check second condition
	mrs_real m = 0.0;
	for(t=0; t < inSamples_; t++)
		m += in(t);
	m /= inSamples_;
*/	
	
	mrs_natural mul = 3; //multiplier proposed in Dixon2006
	mrs_real m = 0.0;
	for(t=checkPoint-(mul*w); t < inSamples_; t++)
		m += in(t);
	m /= (w*4+1);
	
	//checkPoint value should be higher than the window mean and mean should
	//be a significant value (otherwise we most probably are in a silence segment,
	//and we do not want to detect onsets on segments!)
	if(checkPointValue <= (m * ctrl_threshold_->to<mrs_real>()) || m < 10e-20)
	{
		isOnset = false;
		MRSDIAG("PeakerOnset::myProcess() - Failed 2nd condition!");
	}

	//third condition from Dixon2006 (DAFx paper) is not implemented
	//since it was found on that paper that its impact is minimal...

	if(isOnset)
	{
		ctrl_onsetDetected_->setValue(true);
		//ctrl_confidence_->setValue(1.0); //[!] must still find a way to output a confidence...
		ctrl_confidence_->setValue(checkPointValue/100.0); // ad-hoc value which should still have more meaning than a pure 1.0 vs. 0.0.
		out.setval(1.0);
		MRSDIAG("PeakerOnset::myProcess() - Onset Detected!");
	}

	//used for toy_with_onsets.m (DO NOT DELETE! - COMMENT INSTEAD)

//.........这里部分代码省略.........

void
McAulayQuatieri::myProcess(realvec& in, realvec& out)
{
  mrs_natural t,o,c;
  t=0;
  o=0;
  c=0;

  realvec* outPtr;

  out(o,t) = in(o,t);          //    ??????

  //if we want to use memory and we already have data from
  //past inputs (i.e. memory is not empty)...
  if(ctrl_useMemory_->to<mrs_bool>() && memory_.getSize() != 0)
  {
    //concatenate memory column vector with current input
    //so we can continue peak tracking from previous input
    tmp_.stretch(onObservations_, onSamples_+1);
    for(o = 0; o < onObservations_; ++o)
      tmp_(o, 0) = memory_(o);
    for(o = 0; o < onObservations_; ++o)
      for(c = 0; c < onSamples_; ++c)
        tmp_(o,c+1) = in(o,c);
    outPtr = &tmp_;

    //attempt matching of groups between the frame in memory
    //and the first frame from current input
    if(ctrl_useGroups_->to<mrs_bool>())
    {
      peakView inPV(in);
      mrs_realvec inFirstFrame;
      in.getCol(0, inFirstFrame);
      peakView inFirstFramePV(inFirstFrame);
      peakView memPV(memory_);
      peakView tmpPV(tmp_);

      //mrs_natural numInGroups = inPV.getNumGroups();
      mrs_natural numInFirstFrameGroups = inFirstFramePV.getNumGroups();
      mrs_natural numMemGroups = memPV.getNumGroups();

      //we must update the group numbers of the groups
      //in tmp realvec (i.e. [mem|current input])
      //so they do not clash with the previous ones
      if(nextGroup_ > 0)
        for(mrs_natural f=1; f < tmpPV.getNumFrames(); ++f)
          for(mrs_natural p = 0; p < tmpPV.getFrameNumPeaks(f); ++p)
            tmpPV(p, peakView::pkGroup, f) = tmpPV(p, peakView::pkGroup, f) + nextGroup_;

      // Try matching previous groups (in memory from last input)
      // with groups in current input

      // create a tmp copy of the frame in memory and the first frame
      // of current input, so we can do the group matching without
      // destroying the input values
      realvec frames2Match(inObservations_, 2);

      // calculate the matching score for all pairs of groups under matching
      realvec matchScores(numInFirstFrameGroups, numMemGroups);
      for(mrs_natural mg=0; mg < numMemGroups; ++mg)
      {
        for(mrs_natural ig = nextGroup_; ig < nextGroup_ + numInFirstFrameGroups; ++ig)
        {
          //since peakTrack(...) is destructible, we must reset frames2Match everytime... [!]
          for(o=0; o<inObservations_; ++o)
            for(c=0; c < 2; ++c)
              frames2Match(o, c) = tmp_(o, c);

          //use McAulay-Quatieri num of successful peak continuations as a score
          //for the group matching (may be replaced by some other metric in future)
          matchScores(ig-nextGroup_, mg) = peakTrack(frames2Match, 0, ig, mg);
        }
      }

      //Given the matchScores, try to find the optimal assignment
      //of the groups coming from previous input (stored in memory)
      //and the new groups in the current input
      //(using, for e.g. the hungarian method)
      realvec assignedGrp(numInFirstFrameGroups);

      //convert matchScores to costs
      mrs_real maxScore = matchScores.maxval();
      for(o=0; o < matchScores.getRows(); ++o)
        for(c=0; c < matchScores.getCols(); ++ c)
          matchScores(o,c) = maxScore - matchScores(o,c);

      NumericLib::hungarianAssignment(matchScores, assignedGrp); //!!!!!!!!!!!!!!! [TODO][!]

      // given the assignments, try to propagate the group IDs
      // to the groups in the current input
      mrs_natural ig;
      for(mrs_natural f=1; f < tmpPV.getNumFrames(); ++f)
      {
        for(mrs_natural p = 0; p < tmpPV.getFrameNumPeaks(f); ++p)
        {
          //get input group ID (converted to the range [0:...])
          ig = (mrs_natural)(tmpPV(p, peakView::pkGroup, f)) - nextGroup_;

          if(assignedGrp(ig) > -1) //a match was found for this group (ig)
          {
//.........这里部分代码省略.........

void
SelfSimilarityMatrix::myProcess(realvec& in, realvec& out)
{
    if(this->getctrl("mrs_natural/mode")->to<mrs_natural>() ==  SelfSimilarityMatrix::outputDistanceMatrix)
    {
        //check if there are any elements to process at the input
        //(in some cases, they may not exist!) - otherwise, do nothing
        //(i.e. output will be zeroed out)
        if(inSamples_ > 0)
        {
            unsigned int child_count = marsystems_.size();
            if(child_count == 1)
            {
                mrs_natural nfeats = in.getRows();

                //normalize input features if necessary
                if(ctrl_normalize_->to<mrs_string>() == "MinMax")
                    in.normObsMinMax(); // (x - min)/(max - min)
                else if(ctrl_normalize_->to<mrs_string>() == "MeanStd")
                    in.normObs(); // (x - mean)/std

                //calculate the Covariance Matrix from the input, if defined
                if(ctrl_calcCovMatrix_->to<mrs_natural>() & SelfSimilarityMatrix::fixedStdDev)
                {
                    //fill covMatrix diagonal with fixed value (remaining values are zero)
                    MarControlAccessor acc(ctrl_covMatrix_);
                    realvec& covMatrix = acc.to<mrs_realvec>();
                    covMatrix.create(inObservations_, inObservations_);
                    mrs_real var = ctrl_stdDev_->to<mrs_real>();
                    var *= var;
                    for(mrs_natural i=0; i< inObservations_; ++i)
                    {
                        covMatrix(i,i) = var;
                    }
                }
                else if(ctrl_calcCovMatrix_->to<mrs_natural>() & SelfSimilarityMatrix::diagCovMatrix)
                {
                    in.varObs(vars_); //FASTER -> only get the vars for each feature
                    mrs_natural dim = vars_.getSize();
                    //fill covMatrix diagonal with var values (remaining values are zero)
                    MarControlAccessor acc(ctrl_covMatrix_);
                    realvec& covMatrix = acc.to<mrs_realvec>();
                    covMatrix.create(dim, dim);
                    for(mrs_natural i=0; i< dim; ++i)
                    {
                        covMatrix(i,i) = vars_(i);
                    }
                }
                else if(ctrl_calcCovMatrix_->to<mrs_natural>() & SelfSimilarityMatrix::fullCovMatrix)
                {
                    MarControlAccessor acc(ctrl_covMatrix_);
                    realvec& covMatrix = acc.to<mrs_realvec>();
                    in.covariance(covMatrix); //SLOWER -> estimate the full cov matrix
                }
                else if(ctrl_calcCovMatrix_->to<mrs_natural>() == SelfSimilarityMatrix::noCovMatrix)
                {
                    ctrl_covMatrix_->setValue(realvec());
                }

                for(mrs_natural i=0; i < in.getCols(); ++i)
                {
                    in.getCol(i, i_featVec_);
                    for(mrs_natural j=0; j <= i; ++j)
                    {
                        in.getCol(j, j_featVec_);

                        //stack i and j feat vecs
                        for(mrs_natural r=0; r < nfeats; ++r)
                        {
                            stackedFeatVecs_(r, 0) = i_featVec_(r);
                            stackedFeatVecs_(r+nfeats, 0) = j_featVec_(r);
                        }
                        //do the metric calculation for these two feat vectors
                        //and store it in the similarity matrix (which is symmetric)
                        marsystems_[0]->process(stackedFeatVecs_, metricResult_);
                        out(i,j) = metricResult_(0,0);
                        //metric should be symmetric!
                        out(j, i) = out(i, j);
                    }
                }
            }
            else
            {
                out.setval(0.0);
                if(child_count == 0)
                {
                    MRSWARN("SelfSimilarityMatrix::myProcess - no Child Metric MarSystem added - outputting zero similarity matrix!");
                }
                else
                {
                    MRSWARN("SelfSimilarityMatrix::myProcess - more than one Child MarSystem exists (i.e. invalid metric) - outputting zero similarity matrix!");
                }
            }
        }

        //MATLAB_PUT(out, "simMatrix");
        //MATLAB_EVAL("figure(1);imagesc(simMatrix);");

        //MATLAB_PUT(out, "simMat");
        //MATLAB_EVAL(name_+"=["+name_+",simMat(:)'];");
//.........这里部分代码省略.........

void 
OneRClassifier::myProcess(realvec& in, realvec& out)
{
  cout << "OneRClassifier::myProcess" << endl;
  cout << "in.getCols() = " << in.getCols() << endl;
  cout << "in.getRows() = " << in.getRows() << endl;
  //get the current mode, either train of predict mode
  bool trainMode = (getctrl("mrs_string/mode")->to<mrs_string>() == "train");
  row_.stretch(in.getRows());
  if (trainMode)
    {
      if(lastModePredict_ || instances_.getCols()<=0)
	{
	  mrs_natural nAttributes = getctrl("mrs_natural/inObservations")->to<mrs_natural>();
	  cout << "nAttributes = " << nAttributes << endl;
	  instances_.Create(nAttributes);
	}
      
      lastModePredict_ = false;
      
      //get the incoming data and append it to the data table
      for (mrs_natural ii=0; ii< inSamples_; ++ii)
	{
	  mrs_real label = in(inObservations_-1, ii);
	  instances_.Append(in);
	  out(0,ii) = label;
	  out(1,ii) = label;
	}//for t
    }//if
  else
    {//predict mode

	  cout << "OneRClassifier::predict" << endl;
	  if(!lastModePredict_)
	    {
	      //get the number of class labels and build the classifier
	      mrs_natural nAttributes = getctrl("mrs_natural/inObservations")->to<mrs_natural>();
	      cout << "BUILD nAttributes = " << nAttributes << endl;
	      Build(nAttributes);
	    }//if
	  lastModePredict_ = true;
	  cout << "After lastModePredict" << endl;


      //foreach row of predict data, extract the actual class, then call the
      //classifier predict method. Output the actual and predicted classes.
      for (mrs_natural ii=0; ii<inSamples_; ++ii)
	{
	  //extract the actual class
	  mrs_natural label = (mrs_natural)in(inObservations_-1, ii);
	      
	  //invoke the classifier predict method to predict the class
	  in.getCol(ii,row_);
	  mrs_natural prediction = Predict(row_);
	  cout << "PREDICTION = " << prediction << endl;
	  cout << "row_ " << row_ << endl;

	  //and output actual/predicted classes
	  out(0,ii) = (mrs_real)prediction;
	  out(1,ii) = (mrs_real)label;
	}//for t
    }//if
	
}//myProcess

void PeakInObservation::myProcess(realvec& inVec, realvec& outVec)
{
	// (!!) Should be simplified
	outVec.setval(0.f);

	//int nmin = 0;
	mrs_real vmin = inVec(0);
	int nmax = 0;
	mrs_real vmax = inVec(0);

	int nthresh = 0;
	bool theValid = true;
	bool theMaxFlag = true;
	
	for (mrs_natural n = 1; n < inVec.getSize(); n++){
		if (theMaxFlag)
			if (inVec(n) > vmax){
				// Zone 1: [min hysteresis, max]
				vmax = inVec(n);
				nmax = n;
				nthresh = n;
				theValid = true;

				vmin = vmax;
				//nmin = nmax;
			}else{
				if (inVec(n)<vmax/HystFactor_ && nmax!=0){
					// Zone 3: [max hysteresis, min]
					
					if ((mrs_natural)n > nthresh + HystLength_){
						// Maximum was WIDE ENOUGH
						if (theValid){
							outVec(nmax) = vmax;
							theMaxFlag = false;
						}else{
							//Search for new maximum
							vmax = inVec(n);
							nmax = n;
							nthresh = n;
							theValid = true;

							vmin = vmax;
							//nmin = nmax;
						}

					}else{
						// Maximum was TOO SMALL
						if (inVec(n) < vmin){ 
							vmin = inVec(n);
							//nmin = n;
						}
					}
				}else{
					// Zone 2: [max, max hysteresis]
					if (nthresh != (mrs_natural)n-1){
						theValid = false;
						if ((mrs_natural)n > nthresh + HystLength_){
							// Search for new maximum
							vmax = inVec(n);
							nmax = n;
							nthresh = n;
							theValid = true;

							vmin = vmax;
							//nmin = nmax;
						}
					}else
						nthresh = n;
				}
			}
		else
			if (inVec(n) < vmin){ 
				vmin = inVec(n);
				//nmin = n; 
			}else
				if (inVec(n) > vmin*HystFactor_){
					vmax = inVec(n); 
					nmax = n;
					nthresh = 0;

					vmin = vmax;
					//nmin = nmax;
					theValid = true;
					theMaxFlag = true; 
				}
	}
}

void
BeatHistoFeatures::beatHistoFeatures(realvec& in, realvec& out)
{

  mrs_real sum = 0;

  for (mrs_natural o=0; o < inObservations_; o++)
    for (mrs_natural t = 0; t < inSamples_; t++)
    {
      sum += in(o,t);
    }


  mrs_real result[2];
  mrs_natural i,startIdx = 200;
  // zero-out below 50BPM
  for (i=0; i < startIdx; i++)
    in(i) = 0;

  for (i = startIdx; i < in.getCols (); i++)
    if (in(i) < 0)
      in(i) = 0;





  pkr1_->process(in, pkres1_);
  mxr_->process(pkres1_,mxres_);



  vector<mrs_real> bpms;
  bpms.push_back(mxres_(0,1));
  bpms.push_back(mxres_(0,3));
  bpms.push_back(mxres_(0,5));

  sort(bpms.begin(), bpms.end());

  out(0,0) = sum;
  for (unsigned int i=0; i<bpms.size(); i++)
    for (unsigned int j =0; j < bpms.size(); j++)
    {
      if (bpms[i] == mxres_(0,2*j+1))
        out(i+1,0) = mxres_(0,2*j);
    }



  out(4,0) = bpms[0] /4.0;
  out(5,0) = bpms[1] /4.0;
  out(6,0) = bpms[2] /4.0;
  out(7,0) = out(4,0) / out(5,0);



  NormInPlace (in);



#ifdef MARSYAS_MATLAB
#ifdef MTLB_DBG_LOG
  MATLAB_PUT(in, "beathist");
  MATLAB_EVAL("figure(1);plot((201:800)/4, beathist(201:800)),grid on");
#endif
#endif

  MaxAcf (result[0], result[1],in, flag_, startIdx, 600);
  out(8,0)	= result[0];
  out(9,0)	= result[1];
  out(10,0)	= MaxHps (in, startIdx);
  out(11,0)    = SpectralFlatness (in, startIdx);
  out(12,0)	= Std(in);
  out(13,0)	= PeriodicCentroid(in, false, startIdx);
  out(14,0)	= PeriodicCentroid(in, true, startIdx);
  out(15,0)	= PeriodicSpread(in, out(13,0), false, startIdx);
  out(16,0)	= PeriodicSpread(in, out(14,0), true, startIdx);
  out(17,0)	= NumMax(in);

}

void
TimeFreqPeakConnectivity::myProcess(realvec& in, realvec& out)
{
	// get pitch resolution
	const mrs_real	reso		= ctrl_reso_->to<mrs_real>();
	const mrs_bool	isInBark	= getctrl("mrs_bool/inBark")->to<mrs_bool>();

	// a matrix indicating where peaks are in the time frequency plane
	MRSASSERT(textWinSize_ >= in(1,inSamples_-1)-in(1,0));
	peakMatrix_.stretch(numBands_, textWinSize_);

	// init
	peakMatrix_.setval (1);
	for (t = 0; t < textWinSize_; t++)
		for (o=0; o < numBands_; o++)
			peakIndices_[o][t]	= -1;

	// initialized pseudo spectrogram representation
	for (t = 0; t < inSamples_; t++)
	{
		mrs_natural	row = (isInBark)? Freq2RowIdx (in(0,t),reso) : Freq2RowIdx (bark2hertz(in(0,t)),reso),	// input is in bark frequency, so we might have to convert it back
					col = (mrs_natural)(in(1,t)-in(1,0)+.1);
		MRSASSERT(col >= 0 && col < textWinSize_);
		MRSASSERT(row >= 0 && row < numBands_);
		// check whether more than one peak are at that matrix pos, i.e. we already have set the entry
		if (peakIndices_[row][col] != -1)
		{
			multipleIndices->AddIndex (row,col,peakIndices_[row][col]);
			multipleIndices->AddIndex (row,col,t);
			peakIndices_[row][col]	= -2;
		}
		else
			peakIndices_[row][col]	= t; // a matrix indicating which matrix bin corresponds to which input index
		peakMatrix_(row, col)	= 0;
	}

	// initialize output
	out.setval (costInit);

#ifdef MATLAB_DBG_OUT
#ifdef MARSYAS_MATLAB
	MATLAB_PUT(peakMatrix_, "peakMatrix");
	MATLAB_EVAL ("figure(1),imagesc(peakMatrix),colorbar");
#endif
#endif

	// iteration over all pairs
	for (t = 0; t < inSamples_; t++)
	{
		for (o=inSamples_-1; o >= t;o--)
		{
			// don't compute distance if we already have it
			if (out(t,o) != costInit)
				continue;

			// get peak matrix indices
			mrs_natural rowt = (isInBark)? Freq2RowIdx (in(0,t),reso) : Freq2RowIdx (bark2hertz(in(0,t)),reso),	// input is in bark frequency, so we might have to convert it back
						rowo = (isInBark)? Freq2RowIdx (in(0,o),reso) : Freq2RowIdx (bark2hertz(in(0,o)),reso),	// input is in bark frequency, so we might have to convert it back
						colt = (mrs_natural)(in(1,t)-in(1,0)+.1),
						colo = (mrs_natural)(in(1,o)-in(1,0)+.1),
						pathLength;

			MRSASSERT(colt >= 0 && colt < textWinSize_);
			MRSASSERT(colo >= 0 && colo < textWinSize_);
			MRSASSERT(rowt >= 0 && rowt < numBands_);
			MRSASSERT(rowo >= 0 && rowo < numBands_);


			// self similarity and similarity with overlapping peaks
			if ((t == o) || (rowt == rowo && colt == colo))
			{
				SetOutput(out, 0, rowt, colt, rowo, colo);
				continue;
			}

			// check if path calculation makes sense with the current dp step size
			if (abs(rowt - rowo) > abs(colt-colo))
			{
				SetOutput(out, 1, rowt, colt, rowo, colo);
				continue;
			}

			// let's calculate only from left to right
			if (colo < colt)
				continue;

			// dynamic programming
			CalcDp (peakMatrix_, rowt, colt, rowo, colo);
			pathLength	= colo-colt+1;

#ifdef MATLAB_DBG_OUT
#ifdef MARSYAS_MATLAB
			MATLAB_PUT(costMatrix_, "cost");
			MATLAB_EVAL ("figure(2),imagesc(cost,[0 10]),colorbar");
#endif
#endif

			// set cost for this path and all subpaths
			for (mrs_natural i = 0; i < pathLength; i++)
			{
//.........这里部分代码省略.........