C++ realvec::getCols方法代码示例

static void MaxAcf (mrs_real& max, mrs_real& mean, const realvec& beatHistogram, realvec& res,mrs_natural startSearchAt, mrs_natural stopSearchAt)
{
  mrs_natural k,len = beatHistogram.getCols ();

  res.setval(0.);

  // compute ACF
  for (k = startSearchAt; k < stopSearchAt; k++) // this can be optimized
  {
    mrs_real val	= 0;

    for (mrs_natural i = k; i < len; i++)
    {
      val += beatHistogram(i) * beatHistogram(i-k);
    }


    res(k)	= val / (len-k);
  }


  //pkr_->process(in, pkres_);

  max = res.maxval ();

  mean = 1e6*res.mean ();
}

void
PeakClusterSelect::swap(realvec& rv, mrs_natural sample1, mrs_natural sample2, mrs_bool swapColumns)
{
  if( swapColumns ) // swap two columns
  {
    int rows = rv.getRows();
    mrs_real tmp;

    for( int i=0 ; i<rows ; ++i )
    {
      tmp = rv(i, sample1);
      rv(i,sample1) = rv(i,sample2);
      rv(i,sample2) = tmp;
    }
  }
  else // swap two rows
  {
    int cols = rv.getCols();
    mrs_real tmp;

    for( int i=0 ; i<cols ; ++i )
    {
      tmp = rv(sample1,i);
      rv(sample1,i) = rv(sample2,i);
      rv(sample2,i) = tmp;
    }
  }
}

void Map::myProcess(realvec & in, realvec & out)
{
  {
    MarControlAccessor input_access(m_input_ctl);
    realvec & input_data = input_access.to<realvec>();

    assert(input_data.getRows() == in.getRows() &&
           input_data.getCols() == in.getCols());

    input_data = in;
  }

  const realvec & output_data = m_output_ctl->to<realvec>();

  assert(output_data.getRows() == out.getRows() &&
         output_data.getCols() == out.getCols());

  out = output_data;
}

static mrs_real SpectralFlatness (const realvec& beatHistogram, mrs_natural startIdx = 200)
{
  mrs_real    res = 0;
  mrs_natural len = beatHistogram.getCols ();
  //mrs_natural sum = beatHistogram.sum ();

  //beatHistogram	/= sum;
  for (mrs_natural i = startIdx; i < len; i++)
    res		+= log(beatHistogram(i)+1e-6);

  return exp(res/(len-startIdx));
}

mrs_real
QGMMModel::deltaBIC(realvec C1, mrs_natural N1, realvec C2, mrs_natural N2, realvec C, mrs_real lambda)
{
	//matrices should be square and equal sized
	if(C1.getCols() != C2.getCols() && C1.getCols() != C.getCols() &&
		C1.getCols()!= C1.getRows())
	{
		MRSERR("QGMMModel:deltaBIC: matrices should all be squared and equal sized...");
		return MAXREAL; //just a way to sinalize the above error... [!]
	}

	mrs_real res;
	mrs_real N = (mrs_real)(N1 + N2);
	mrs_real d = (mrs_real)C1.getCols();


	res = N * log(C.det());
	res -= (mrs_real)N1 * log(C1.det());
	res -= (mrs_real)N2 * log(C2.det());
	res *= 0.5f;

	res -= 0.5f * lambda * (d + 0.5f*d*(d+1.0f))* log(N);

	return res;
}

static mrs_real PeriodicSpread (realvec& vector, mrs_real centroid, mrs_bool isLog = false, mrs_natural startIdx = 200)
{
  mrs_natural len = vector.getCols ();
  mrs_real res = 0,
           norm = 0;

  for (mrs_natural i = startIdx; i < len; i++)
  {
    mrs_real theta = (isLog)? log(i*1./startIdx)* TWOPI :(i * TWOPI) / startIdx;
    res		+= vector(i) * abs(.5*(cos(theta)+1)-centroid);
    norm	+= vector(i);
  }

  return res/norm;
}

void
MATLABengine::putVariable(realvec value, mrs_string MATLABname)
{
  //----------------------------------
  // send a realvec to a MATLAB matrix
  //----------------------------------
  mwSize dims[2]; //realvec is 2D
  dims[0] = value.getRows();
  dims[1] = value.getCols();

  //realvec are by default double precision matrices => mxDOUBLE_CLASS
  mxArray *mxMatrix = mxCreateNumericArray(2, dims, mxDOUBLE_CLASS, mxREAL);
  mrs_real *data = value.getData();
  memcpy((void *)mxGetPr(mxMatrix), (void *)(data), dims[0]*dims[1]*mxGetElementSize(mxMatrix));
  engPutVariable(engine_, MATLABname.c_str(), mxMatrix);

  mxDestroyArray(mxMatrix);
}

void PeakConvert2::ComputePeaker (mrs_realvec in, realvec& out)
{
#ifdef ORIGINAL_VERSION
  peaker_->updControl("mrs_real/peakStrength", 0.2);// to be set as a control [!]
#else
  peaker_->updControl("mrs_real/peakStrength",1e-1);
  peaker_->updControl("mrs_real/peakStrengthRelMax" ,1e-2);
  peaker_->updControl("mrs_real/peakStrengthAbs",1e-10 );
  peaker_->updControl("mrs_real/peakStrengthTreshLpParam" ,0.95);
  peaker_->updControl("mrs_real/peakStrengthRelThresh" , 1.);
#endif

  peaker_->updControl("mrs_real/peakSpacing", 2e-3);   // 0
  peaker_->updControl("mrs_natural/peakStart", downFrequency_);   // 0
  peaker_->updControl("mrs_natural/peakEnd", upFrequency_);  // size_
  peaker_->updControl("mrs_natural/inSamples", in.getCols());
  peaker_->updControl("mrs_natural/inObservations", in.getRows());
  peaker_->updControl("mrs_natural/onSamples", out.getCols());
  peaker_->updControl("mrs_natural/onObservations", out.getRows());

  peaker_->process(in, out);
}

static mrs_real MaxHps (const realvec& beatHistogram, mrs_natural startIdx = 200)
{
  const mrs_natural order = 4;
  mrs_natural k,len = beatHistogram.getCols ();
  mrs_realvec	res = beatHistogram;	// make this a member

  //res.setval(-1e38);

  for (k =2; k < order; k++)
  {
    for (mrs_natural i = startIdx; i < len; i++)
    {
      if (k*i >= len)
        break;
      res(i)		+= log(beatHistogram(k*i)+1e-6);
    }
  }

  for (k = 0; k < startIdx; k++)
    res(k)	= -1e38;


  return exp (res.maxval ());
}

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

summaryStatistics ClassificationReport::computeSummaryStatistics(const realvec& mat)
{
  MRSASSERT(mat.getCols()==mat.getRows());

  summaryStatistics stats;

  mrs_natural size = mat.getCols();

  vector<mrs_natural>rowSums(size);
  for(int ii=0; ii<size; ++ii) rowSums[ii] = 0;
  vector<mrs_natural>colSums(size);
  for(int ii=0; ii<size; ++ii) colSums[ii] = 0;
  mrs_natural diagonalSum = 0;

  mrs_natural instanceCount = 0;
  for(mrs_natural row=0; row<size; row++)
  {
    for(mrs_natural col=0; col<size; col++)
    {
      mrs_natural num = (mrs_natural)mat(row,col);
      instanceCount += num;

      rowSums[row] += num;
      colSums[col] += num;

      if(row==col)
        diagonalSum += num;
    }
  }
  //printf("row1 sum:%d\n",rowSums[0]);
  //printf("row2 sum:%d\n",rowSums[1]);
  //printf("col1 sum:%d\n",colSums[0]);
  //printf("col2 sum:%d\n",colSums[1]);
  //printf("diagonal sum:%d\n",diagonalSum);
  //printf("instanceCount:%d\n",instanceCount);

  mrs_natural N = instanceCount;
  mrs_natural N2 = (N*N);
  stats.instances = instanceCount;
  stats.correctInstances = diagonalSum;

  mrs_natural sum = 0;
  for(mrs_natural ii=0; ii<size; ++ii)
  {
    sum += (rowSums[ii] * colSums[ii]);
  }
  mrs_real PE = (mrs_real)sum / (mrs_real)N2;
  mrs_real PA = (mrs_real)diagonalSum / (mrs_real)N;
  stats.kappa = (PA - PE) / (1.0 - PE);

  mrs_natural not_diagonal_sum = instanceCount - diagonalSum;
  mrs_real MeanAbsoluteError = (mrs_real)not_diagonal_sum / (mrs_real)instanceCount;
  //printf("MeanAbsoluteError:%f\n",MeanAbsoluteError);
  stats.meanAbsoluteError = MeanAbsoluteError;

  mrs_real RootMeanSquaredError = sqrt(MeanAbsoluteError);
  //printf("RootMeanSquaredError:%f\n",RootMeanSquaredError);
  stats.rootMeanSquaredError = RootMeanSquaredError;

  mrs_real RelativeAbsoluteError = (MeanAbsoluteError / 0.5) * 100.0;
  //printf("RelativeAbsoluteError:%f%%\n",RelativeAbsoluteError);
  stats.relativeAbsoluteError = RelativeAbsoluteError;

  mrs_real RootRelativeSquaredError = (RootMeanSquaredError / (0.5)) * 100.0;
  //printf("RootRelativeSquaredError:%f%%\n",RootRelativeSquaredError);
  stats.rootRelativeSquaredError = RootRelativeSquaredError;

  return stats;
}//computeSummaryStatistics

void 
PeakViewMerge::myProcess(realvec& in, realvec& out)
{
	peakView	*In[kNumMatrices],
				Out (out);
	mrs_natural i, rowIdx = 0,
				numPeaks[kNumMatrices],
				outputIdx	= 0;
	const mrs_bool discNegGroups	= ctrl_noNegativeGroups_->to<mrs_bool>();

	out.setval(0.);
	
	for (i = 0; i < kNumMatrices; i++)
	{
		mrs_natural	numRows		= (i==kMat1)? ctrl_frameMaxNumPeaks1_->to<mrs_natural>() :  ctrl_frameMaxNumPeaks2_->to<mrs_natural>();
		numRows					*= peakView::nbPkParameters;
		if (numRows == 0) // if the controls have not been set assume both matrixes to be of equal size	
			numRows	= in.getRows ()/kNumMatrices;
		peakViewIn_[i].stretch (numRows, in.getCols ());
		in.getSubMatrix (rowIdx, 0, peakViewIn_[i]);
		rowIdx		+= numRows;
		In[i]		= new peakView(peakViewIn_[i]);
		numPeaks[i]	= In[i]->getTotalNumPeaks ();
	}

	if (ctrl_mode_->to<mrs_string>() == "OR")
	{
		// write all entries of the second peakView to output
		for (i = 0; i < numPeaks[1]; i++)
		{
			if (discNegGroups && (*In[1])(i,peakView::pkGroup) < 0)
				continue;
			WriteOutput (Out, In[1], i, outputIdx);
			outputIdx++;
		}

		// write all entries of the first peakView to output except duplicates
		for (i = 0; i < numPeaks[0]; i++)
		{
			mrs_natural Idx;
			if (discNegGroups && (*In[0])(i,peakView::pkGroup) < 0)
				continue;
			for (mrs_natural k = 1; k < kNumMatrices; k++)
				Idx	= FindDuplicate (In[k], (*In[0])(i, peakView::pkFrequency), numPeaks[k]);

			if (Idx < 0)
			{
				WriteOutput (Out, In[0], i, outputIdx);
				outputIdx++;
			}
		}
	}
	else if (ctrl_mode_->to<mrs_string>() == "AND")
	{
		// find duplicates and write only them to output
		for (i = 0; i < numPeaks[0]; i++)
		{
			mrs_natural Idx;
			if (discNegGroups && (*In[0])(i,peakView::pkGroup) < 0)
				continue;
			for (mrs_natural k = 1; k < kNumMatrices; k++)
				Idx	= FindDuplicate (In[k], (*In[0])(i, peakView::pkFrequency), numPeaks[k]);

			if (Idx >= 0)
			{
				if (discNegGroups && (*In[1])(Idx,peakView::pkGroup) < 0)
					continue;
				WriteOutput (Out, In[0], i, outputIdx);
				outputIdx++;
			}
		}
	}
	else if (ctrl_mode_->to<mrs_string>() == "ANDOR")
	{
		// keep the input[0] peaks that are not in input[1]
		for (i = 0; i < numPeaks[0]; i++)
		{
			mrs_natural Idx;
			if (discNegGroups && (*In[0])(i,peakView::pkGroup) < 0)
				continue;
			for (mrs_natural k = 1; k < kNumMatrices; k++)
				Idx	= FindDuplicate (In[k], (*In[0])(i, peakView::pkFrequency), numPeaks[k]);

			if (Idx < 0)
			{
				WriteOutput (Out, In[0], i, outputIdx);
				outputIdx++;
			}
		}
	}
	else if (ctrl_mode_->to<mrs_string>() == "XOR")
	{
		// find duplicates and write only residual to output
		for (i = 0; i < numPeaks[0]; i++)
		{
			if (discNegGroups && (*In[0])(i,peakView::pkGroup) < 0)
				continue;
			mrs_natural Idx	= FindDuplicate (In[1], (*In[0])(i, peakView::pkFrequency), numPeaks[1]);

			if (Idx < 0)
//.........这里部分代码省略.........

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
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
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)
  //
}