diff --git a/VespucciLibrary/VespucciLibrary.pro b/VespucciLibrary/VespucciLibrary.pro
index 105cd0da..1406abf5 100644
--- a/VespucciLibrary/VespucciLibrary.pro
+++ b/VespucciLibrary/VespucciLibrary.pro
@@ -81,6 +81,7 @@ SOURCES +=\
     src/Math/Fitting/nonlinleastsq.cpp \
     src/Math/Accessory/Faddeeva.cpp \
     src/Math/Baseline/rollingball.cpp \
+    src/Math/Baseline/goldindec.cpp \
     src/Math/Smoothing/denoise.cpp \
     src/Math/Clustering/kmeanswrapper.cpp \
     src/vespucci.cpp
@@ -116,6 +117,8 @@ HEADERS  += \
     include/Math/Fitting/nonlinleastsq.h \
     include/Math/Accessory/Faddeeva.h \
     include/Math/Baseline/rollingball.h \
+    include/Math/Baseline/mconvert.h \
+    include/Math/Baseline/goldindec.h \
     include/Math/Smoothing/denoise.h \
     include/Math/Clustering/kmeanswrapper.h \
     include/vespucci.h \
diff --git a/VespucciLibrary/include/Math/Baseline/goldindec.h b/VespucciLibrary/include/Math/Baseline/goldindec.h
new file mode 100644
index 00000000..c29b2a69
--- /dev/null
+++ b/VespucciLibrary/include/Math/Baseline/goldindec.h
@@ -0,0 +1,15 @@
+#ifndef GOLDINDEC_H
+#define GOLDINDEC_H
+using namespace arma ;
+namespace Spec{
+    namespace Baseline{
+        vec GoldindecBaseline(vec X,
+                              vec Y,
+                              uword p,
+                              double peak_ratio,
+                              double eps);
+        vec LEGEND_C(vec y, vec n, uword ord, double s);
+        double T_rate(double x);
+    }
+}
+#endif // GOLDINDEC_H
diff --git a/VespucciLibrary/include/Math/Baseline/mconvert.h b/VespucciLibrary/include/Math/Baseline/mconvert.h
new file mode 100644
index 00000000..4af8de4e
--- /dev/null
+++ b/VespucciLibrary/include/Math/Baseline/mconvert.h
@@ -0,0 +1,445 @@
+#ifndef MCONVERT_H
+#define MCONVERT_H
+
+#include <armadillo>
+using namespace arma;
+
+namespace m2cpp {
+
+    template<typename eT>
+    inline arma::Col<eT> scol(eT a) {
+        return arma::Col<eT>(&a, 1, true);
+    }
+
+    template<typename eT>
+    inline arma::Row<eT> srow(eT a) {
+        return arma::Row<eT>(&a, 1, true);
+    }
+
+    template<typename eT>
+    inline arma::Mat<eT> smat(eT a) {
+        return arma::Mat<eT>(&a, 1, 1, true);
+    }
+
+
+    inline arma::uvec span(int a, int b) {
+        arma::uvec s;
+        int n = b - a;
+        if (n < 0) return s;
+        s.set_size(n + 1);
+        for (int ii = 0; ii <= n; ii++)
+            s(ii) = ii + a;
+        return s;
+    }
+
+    template <typename T>
+    inline T span(int a, int b) {
+        T s;
+        int n = b - a;
+        if (n < 0) return s;
+        s.set_size(n + 1);
+        for (int ii = 0; ii <= n; ii++)
+            s(ii) = ii + a;
+        return s;
+    }
+
+    inline arma::uvec span(int a, int step, int b)
+    {
+        arma::uvec s;
+        int n = (b - a + step) / step;
+        if (n < 0)
+        {
+            return s;
+        }
+        s.set_size(n);
+
+        for (int ii = 0; ii < n; ii++)
+        {
+            s(ii) = step * ii + a;
+        }
+
+        return s;
+    }
+
+    template <typename T>
+    inline T span(int a, int step, int b)
+    {
+        T s;
+        int n = (b - a + step) / step;
+        if (n < 0)
+        {
+            return s;
+        }
+        s.set_size(n);
+
+        for (int ii = 0; ii < n; ii++)
+        {
+            s(ii) = step * ii + a;
+        }
+
+        return s;
+    }
+
+
+    template <typename T>
+    inline arma::cx_mat fft(arma::Mat<typename T::elem_type> X, int dim)
+    {
+        if (dim == 1)
+            return arma::fft(X);
+        else
+            return arma::strans(arma::fft(arma::strans(X)));
+
+    }
+
+    template <typename T>
+    inline arma::cx_mat fft(arma::Mat<typename T::elem_type> X, int n, int dim)
+    {
+        if (dim == 1)
+            return arma::fft(X, n);
+        else
+            return arma::strans(arma::fft(arma::strans(X), n));
+
+    }
+
+
+    inline arma::cx_mat ifft(arma::cx_mat X, int dim)
+    {
+        if (dim == 1)
+            X = arma::ifft(X);
+        else
+            X = arma::strans(arma::ifft(arma::strans(X)));
+        return X;
+    }
+
+    inline arma::cx_mat ifft(arma::cx_mat X, int n, int dim)
+    {
+        if (dim == 1)
+            X = arma::ifft(X, n);
+        else
+            X = arma::strans(arma::ifft(arma::strans(X), n));
+        return X;
+    }
+    
+
+    //template<typename eT>
+    inline rowvec fspan(double a, double step, double b) {
+        //arma::Col<eT> s;
+        rowvec s;
+        int n = (int) ((b - a) / step);
+        if (n < 0) return s;
+
+        //s.set_size(n + 1);
+        s.set_size(n+1);
+        for (int ii = 0; ii <= n; ii++)
+            s(ii) = step * ii + a;
+        return s;
+    }
+
+    inline int nextpow2(int n) {
+        n = abs(n);
+        int p = 0;
+        int tmp = 1;
+        while (tmp < n)
+        {
+            tmp *= 2;
+            p++;
+        }
+        return p;
+    }
+
+    template<typename T1>
+    inline T1 square(T1 a) {
+        return a*a;
+    }
+
+    template<typename T1, typename T2>
+    inline  arma::Mat<typename T1::elem_type> hankel(const T1& c_, const T2& r_) {
+
+        typedef typename T1::elem_type eT;
+
+        int nc = r_.n_elem;
+        int nr = c_.n_elem;
+
+        const arma::Col<eT> c((eT*)c_.memptr(), nr, 0);
+        const arma::Row<eT> r((eT*)r_.memptr(), nc, 0);
+
+        if (r[0] != c[0]) {
+            //("hankel: differing diagonal element. Using the column one");
+        }
+
+        arma::Mat<typename T1::elem_type> retval(nr, nc);
+
+        for (int i = 1; i <= std::min(nr, nc); i++) {
+            retval.submat(1 - 1, i - 1, nr - i + 1 - 1, i - 1) = c.rows(i - 1, nr
+                                                                        - 1);
+        }
+        int tmp = 1;
+        if (nc <= nr) {
+            tmp = nr - nc + 2;
+        }
+
+        for (int i = nr; i >= tmp; i--) {
+            retval.submat(i - 1, 2 + nr - i - 1, i - 1, nc - 1) = r.cols(2
+                                                                         - 1, nc - nr + i - 1);
+        }
+        return retval;
+    }
+
+    template<typename T1>
+    inline arma::uword length(const T1& A) {
+        return A.n_elem;
+    }
+
+    /*
+    template<typename T>
+    inline arma::Mat <typename T::elem_type> convmtx(const T& v, int m) {
+
+        arma::Mat<typename T::elem_type> out = zeros(v.n_elem + m - 1, m);
+        arma::Col<typename T::elem_type> aux((typename T::elem_type*)v.memptr(), v.n_elem);
+
+        for (int ii = 0; ii < m; ii++) {
+            out.submat(ii, ii, ii + v.n_elem - 1, ii) = v;
+        }
+
+        if (v.n_rows == 1)
+            out = out.t();
+        return out;
+    }
+    */
+
+    template<typename T>
+    inline arma::Mat <typename T::elem_type> convmtx(const T& v, int m) {
+
+        arma::Mat<typename T::elem_type> out = zeros<arma::Mat<typename T::elem_type> >(v.n_elem + m - 1, m);
+        arma::Col<typename T::elem_type> aux((typename T::elem_type*)v.memptr(), v.n_elem);
+
+        if (v.n_rows == 1)
+        {
+            for (int ii = 0; ii < m; ii++) {
+                out.submat(ii, ii, ii + v.n_elem - 1, ii) = aux;
+            }
+        }
+        else
+        {
+            for (int ii = 0; ii < m; ii++) {
+                out.submat(ii, ii, ii + v.n_elem - 1, ii) = v;
+            }
+        }
+
+        if (v.n_rows == 1)
+            out = out.t();
+        return out;
+    }
+
+    template <typename eT>
+    inline typename arma::enable_if2< arma::is_real<typename eT::elem_type>::value, typename arma::Mat<typename eT::elem_type> >::result
+    conv2(const arma::Mat<typename eT::elem_type>& A, const arma::Mat<typename eT::elem_type>& B) {
+        uword n = A.n_rows + B.n_rows - 1;
+        uword m = A.n_rows + B.n_rows - 1;
+
+        arma::Mat<typename eT::elem_type> out = arma::real(arma::ifft2(fft2(A, n, m) % fft2(B, n, m)));
+        return out;
+    }
+
+    template <typename eT, typename fT>
+    inline typename arma::enable_if2 < arma::is_complex<typename eT::elem_type>::value || arma::is_complex<typename fT::elem_type>::value,
+                                       arma::Mat<typename std::complex< typename arma::get_pod_type<eT>::result > > >::result
+        conv2(const arma::Mat<typename eT::elem_type>& A, const arma::Mat<typename fT::elem_type>& B) {
+        uword n = A.n_rows + B.n_rows - 1;
+        uword m = A.n_rows + B.n_rows - 1;
+
+        arma::Mat<typename eT::elem_type> out = arma::ifft2(fft2(A, n, m) % fft2(B, n, m));
+        return out;
+    }
+
+    template <typename eT>
+    inline eT fix(const eT a) {
+       return a > eT(0) ? floor(a) : ceil(a);
+    }
+
+    template<typename T>
+    inline void intersect(arma::Col<typename T::elem_type>& C, arma::uvec& ia, arma::uvec& ib, const T& a, const T& b) {
+
+       typedef typename T::elem_type eT;
+
+       arma::uvec sa = arma::sort_index(a);
+       arma::uvec sb = arma::sort_index(b);
+
+       std::vector<eT> C_;
+       std::vector<arma::uword> ia_, ib_;
+
+       int na = int(a.n_elem);
+       int nb = int(b.n_elem);
+       int ja = 0, jb = 0;
+
+       for (;;) {
+
+          arma::uword sja = sa(ja);
+          arma::uword sjb = sb(jb);
+
+          eT ca = a(sja);
+          eT cb = b(sjb);
+
+          if (ca > cb) {
+             ja++;
+          }
+          else if (cb > ca) {
+             jb++;
+          }
+          else {
+             C_.push_back(ca);
+             ia_.push_back(sja);
+             ib_.push_back(sjb);
+             while (++ja < na && a(sa(ja)) == ca) {}
+             while (++jb < nb && b(sb(jb)) == cb) {}
+             if (ja == na || jb == nb)
+                break;
+          }
+
+       }
+
+       ia = arma::uvec(ia_) + 1;
+       ib = arma::uvec(ib_) + 1;
+       C = arma::Col<eT>(C_);
+    }
+
+    template<typename T>
+    inline void intersect(arma::Col<typename T::elem_type>& C, const T& a, const T& b) {
+       arma::uvec dum0, dum1;
+       intersect(C, dum0, dum1, a, b);
+    }
+
+    template<typename Tr, typename T>
+    inline void unique_rows_core(Tr& C, Tr& a_sorted, arma::uvec& ia, const T& a) {
+
+       typedef typename T::elem_type eT;
+
+       int ic_cur = 0;
+
+       auto cmp = [&a, &ic_cur](const int k, const int l) -> bool {
+          return a(k, ic_cur) < a(l, ic_cur);
+       };
+
+       int nr = int(a.n_rows);
+       int nc = int(a.n_cols);
+
+       arma::Col<eT> ac0(const_cast<eT*>(a.colptr(0)), nr, false);
+       arma::uvec ord = arma::sort_index(a.col(0));
+       std::vector<arma::uword> ia0, ia1;
+
+       std::vector<arma::uword>* ia_ = &ia0;
+       std::vector<arma::uword>* ia_prev_ = &ia1;
+
+       ia_->push_back(0);
+       for (int ii = 1; ii < nr; ii++) {
+          if (a(ord(ii), 0) != a(ord(ii - 1))) {
+             ia_->push_back(ii);
+          }
+       }
+       ia_->push_back(nr);
+
+       for (int ic = 1; ic < nc; ic++) {
+
+          ic_cur = ic;
+          int ir = 0, ir_prev = 0;
+          std::swap(ia_prev_, ia_);
+          int na = int(ia_prev_->size());
+          ia_->clear();
+
+          for (int ii = 0; ii < na - 1; ii++) {
+             ia_->push_back((*ia_prev_)[ii]);
+             int l = (*ia_prev_)[ii], u = (*ia_prev_)[ii + 1];
+             std::sort(&ord(l), &(ord(u - 1)) + 1, cmp);
+             for (int jj = l + 1; jj < u; jj++) {
+                if (a(ord(jj - 1), ic) != a(ord(jj), ic)) {
+                   ia_->push_back(jj);
+                }
+             }
+          }
+
+          ia_->push_back(nr);
+
+       }
+
+       ia = arma::uvec(*ia_);
+       int na = int(ia.n_elem);
+       C.set_size(na - 1, a.n_cols);
+
+       for (int ii = 0; ii < na - 1; ii++) {
+          C.row(ii) = a.row(ord(ia(ii)));
+       }
+
+       a_sorted.set_size(nr, nc);
+       for (int ir = 0; ir < nr; ir++) {
+          a_sorted.row(ir) = a.row(ord(ir));
+       }
+    }
+
+    template<typename T>
+    inline T sortrows(const T& a) {
+
+       typedef typename T::elem_type eT;
+       arma::uvec dum0;
+       arma::Mat<eT> dum1;
+       arma::Mat<eT> ret;
+       unique_rows_core(dum1, ret, dum0, a);
+       return ret;
+    }
+
+    template<typename T>
+    inline void unique_rows(T& C, const T& a) {
+
+       arma::uvec dum;
+       unique_rows(C, dum, a);
+
+    }
+
+    template<typename T>
+    inline T unique_rows(const T& a) {
+       T ret;
+       unique_rows(ret, a);
+       return ret;
+    }
+
+    template<typename T>
+    inline int isempty(const T& a) {
+        return a.n_elem == 0;
+    }
+
+    template<typename T>
+    inline void unique(T& a, const T& b, const char* m) {
+        T tmp;
+        const T& in = b;
+        if (&a == &b) {
+           tmp = b;
+           in = tmp;
+        }
+
+        if (strcmp(m, "rows") == 0) {
+           unique_rows(a, tmp);
+        }
+        else {
+           fprintf(stderr, "m2pp::unique(): Unrecognized option %s\n", m);
+        }
+    }
+
+    static arma::wall_clock timer_;
+
+    inline double tic() {
+       timer_.tic();
+       return timer_.toc();
+    }
+
+    inline double toc() {
+       return timer_.toc();
+    }
+
+    inline double toc(double start) {
+       return timer_.toc() - start;
+    }
+}
+
+
+#endif
+
+
diff --git a/VespucciLibrary/src/Math/Baseline/goldindec.cpp b/VespucciLibrary/src/Math/Baseline/goldindec.cpp
new file mode 100644
index 00000000..47b12397
--- /dev/null
+++ b/VespucciLibrary/src/Math/Baseline/goldindec.cpp
@@ -0,0 +1,107 @@
+#include "mconvert.h"
+#include "goldindec.h"
+///
+/// \brief Spec::Smoothing::Goldindec
+/// \param vec X Input arma::vec
+/// \param vec Y Input arma::vec
+/// \param p poly_order Polynomail order
+/// \param peak_ratio
+/// \param eps
+/// \return vec baseline
+/// 
+// References: Goldindec: A Novel Algorithm for Raman Spectrum Baseline Correction
+vec Spec::Baseline::GoldindecBaseline(vec X, vec Y, uword p, double peak_ratio, double eps)
+{
+  double OLD_s, OPT_s, s, udr, up_down_rate, a, b ;
+  uword t;
+  vec OPT_z, n, y, z ;
+  uvec i;
+  y = Y ;
+  n = X ;
+  a = 0.0 ;
+  b = 1.0 ;
+  udr = T_rate(peak_ratio) ;
+  s = a+0.382*(b-a) ;
+  z = LEGEND_C(y, n, p, s) ;
+  up_down_rate = m2cpp::length(uvec(find(y>=z)))*1.0/m2cpp::length(uvec(find(y<z))) ;
+  t = 0 ;
+  while (std::abs(up_down_rate-udr)>eps)
+  {
+    OLD_s = s ;
+    if (up_down_rate-udr>eps)
+    {
+      a = s ;
+      s = a+0.382*(b-a) ;
+      z = LEGEND_C(y, n, p, s) ;
+      up_down_rate = m2cpp::length(uvec(find(y>=z)))*1.0/m2cpp::length(uvec(find(y<z))) ;
+    }
+    else
+    {
+      b = s ;
+      s = a+0.382*(b-a) ;
+      z = LEGEND_C(y, n, p, s) ;
+      up_down_rate = m2cpp::length(uvec(find(y>=z)))*1.0/m2cpp::length(uvec(find(y<z))) ;
+    }
+    if (std::abs(OLD_s-s)<0.00001)
+      break ;
+    t = t+1 ;
+  }
+  OPT_s = s ;
+  OPT_z = LEGEND_C(y, n, p, OPT_s) ;
+  return OPT_z ;
+}
+
+vec Spec::Baseline::LEGEND_C(vec y, vec n, uword ord, double s)
+{
+  double alpha, dely, maxy ;
+  mat T, Tinv ;
+  uword N, num, it;
+  vec a, d, res, z, zp ;
+  rowvec p;
+  uvec i, j;
+  N = m2cpp::length(n) ;
+  n = sort(n) ;
+  i = sort_index(n) ;
+  y = y(i) ;
+  maxy = arma::max(y) ;
+  dely = (maxy-arma::min(y))/2.0 ;
+  n = 2*(n-n(N-1))/(n(N-1)-n(0))+1 ;
+  y = (y-maxy)/dely+1 ;
+  p = m2cpp::fspan(0, 1, ord);
+  mat n_temp = repmat(n, 1, ord+1);
+  for(num=0;num<ord+1;num++)
+      T.insert_cols(num, arma::pow(n_temp.col(num), p[num]));
+  //T = arma::pow(repmat(n, 1, ord+1), repmat(p, N, 1)) ;
+  Tinv = pinv(arma::trans(T)*T)*arma::trans(T) ;
+  a = Tinv*y ;
+  z = T*a ;
+  alpha = 0.99*1/2.0 ;
+  it = 0 ;
+  zp = arma::ones<vec>(N) ;
+  d = arma::zeros<vec>(N) ;
+  while (arma::as_scalar(arma::sum(arma::square((z-zp))))/arma::as_scalar(arma::sum(arma::square(zp)))>1e-9)
+  {
+    it = it+1 ;
+    zp = z ;
+    res = y-z ;
+    for (num=0; num<N; num++)
+    {
+      if (res[num]<s)
+        d[num] = res[num]*(2*alpha-1) ;
+      else
+        d[num] = -res[num]-alpha*(pow(s, 3))/(2*res[num]*res[num]) ;
+    }
+    a = Tinv*(y+d) ;
+    z = T*a ;
+  }
+  j = sort_index(i) ;
+  z = (z(j)-1)*dely+maxy ;
+  return z ;
+}
+
+double Spec::Baseline::T_rate(double x)
+{
+  double y ;
+  y = 0.7679+11.2358*x-39.7064*pow(x, 2)+92.3583*pow(x, 3) ;
+  return y ;
+}
diff --git a/VespucciLibrary/src/Math/PeakFinding/peakfinding.cpp b/VespucciLibrary/src/Math/PeakFinding/peakfinding.cpp
index b98a5149..20b61992 100644
--- a/VespucciLibrary/src/Math/PeakFinding/peakfinding.cpp
+++ b/VespucciLibrary/src/Math/PeakFinding/peakfinding.cpp
@@ -54,20 +54,21 @@ arma::vec Vespucci::Math::PeakFinding::FindPeaks(arma::vec X,
     dX.elem( find(dX == 0) ).fill(-arma::datum::eps); //to find first of repeated values
     double minimum_magnitude = X.min();
     double left_min = X(0);
-    arma::uword length = X.n_elem;
+    arma::uword length = dX.n_elem; // the length of dX is less one than X
     arma::uword temporary_location = 0;
 
     //Find where derivative changes sign:
     arma::uvec extrema_indices = find( (dX.subvec(0, length - 2) % dX.subvec(1, length - 1) )< 0);
+    extrema_indices = extrema_indices + 1; //from 1 start
     arma::vec X_extrema = X.elem(extrema_indices);
     length = X_extrema.n_elem;
     double temporary_magnitude = minimum_magnitude;
     bool peak_found = false;
-    arma::uword ii;
-    if (X(1) >= X(2))
-        ii = 1;
-    else
-        ii = 2;
+    arma::uword ii, ip;
+
+    ii = 0;
+    ip = 0;
+    peak_locations = uvec(length);
 
     //This loop finds the peak locations
     while (ii < length){
@@ -92,15 +93,18 @@ arma::vec Vespucci::Math::PeakFinding::FindPeaks(arma::vec X,
             peak_found = true;
             std::cout << "peak found!" << std::endl;
             left_min = X_extrema(ii);
-            peak_locations << temporary_location;
+            temporary_location = ii;
+            peak_locations[ip] = temporary_location;
+            ++ip;
         }
-        else if(X(ii) < left_min)
+        else if(X_extrema(ii) > left_min)
             left_min = X_extrema(ii);
 
         ++ii;
 
 
     }
+    peak_locations = peak_locations.head(ip);
     //remove all peaks below the threshold
 
     arma::vec results = arma::zeros(X.n_elem, 1);