diff --git a/matlab.h b/matlab.h
new file mode 100644
index 0000000000000000000000000000000000000000..27d163c6a99d83a3f11d77eae03367db11824da2
--- /dev/null
+++ b/matlab.h
@@ -0,0 +1,241 @@
+#pragma once
+
+#include <random>
+#include <iostream>
+#include <fstream>
+#include "../types.h"
+#include <complex>
+/*
+ * Code translated from Matlab base functions
+ * matlab translation guide : https://eigen.tuxfamily.org/dox/AsciiQuickReference.txt
+ */
+
+//-------------------------------------------------------------------------------------------------
+// VARIOUS
+
+const number pi = M_PI;
+
+// squares a number
+template<typename T> inline T squared(const T& x)
+{
+ return x*x;
+}
+
+//-------------------------------------------------------------------------------------------------
+// ALGEBRA
+
+// build a vector with values in the given interval
+Vector make_vector(const int fromValue, const int toValue, const int by = 1)
+{
+ int size = 1 + (toValue - fromValue) / by;
+ Vector result(size);
+
+ int currentValue = fromValue;
+ for(int i = 0; i < size; i++)
+ {
+ result(i) = currentValue;
+ currentValue += by;
+ }
+
+ return result;
+}
+
+// repeats a column vector n times in order to build a matrix
+Matrix repmat(const Vector& v, const int colNumber)
+{
+ Matrix result(v.size(), colNumber);
+
+ for(int col = 0; col < colNumber; col++)
+ {
+ result.col(col) = v;
+ }
+
+ return result;
+}
+
+// solves m = m1 \ m2
+inline ComplexMatrix mldivide(const ComplexMatrix& m1, const ComplexMatrix& m2)
+{
+ //return m1.colPivHouseholderQr().solve(m2);
+ return m1.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(m2); // you can use JacobiSVD for improved precision
+}
+
+// returns an hermitian matrix where v is the first row and conj(v) the first column
+Matrix toeplitz(const Vector& v)
+{
+ auto size = v.size();
+ Matrix result(size, size);
+
+ for(int row=0; row<size; row++)
+ {
+ for(int col=0; col<size; col++)
+ {
+ // complex
+ /*number value = v(std::abs(row - col));
+ if(row > col)
+ {
+ value = Sstd::conj(value);
+ }
+ result(row,col) = value;*/
+
+ // real
+ result(row,col) = v(std::abs(row - col));
+ }
+ }
+
+ return result;
+}
+
+// flip the column and rows
+inline Matrix flip(Matrix mat)
+{
+ Matrix matflipr = mat.rowwise().reverse(); // fliplr
+ Matrix matflipc = matflipr.colwise().reverse(); // flipud
+ return matflipc;
+}
+
+//-------------------------------------------------------------------------------------------------
+// POLYNOMIAL
+
+/*
+ * evaluates a polynomial on the given data
+ * use Horner's scheme
+ * NOTE: one could use a more numerically stable algorithm
+ */
+Vector polyval(const Polynomial& pol, const Vector& x)
+{
+ Vector result = Vector::Zero(x.size());
+
+ for(const auto& coef : pol)
+ {
+ result = result.binaryExpr(x, [coef](const number& result, const number& x){return result*x + coef;});
+ }
+
+ return result;
+}
+
+//-------------------------------------------------------------------------------------------------
+
+// RANDOM HELPING FUNCTION
+
+const unsigned int n = 65535;//0xFFFF;
+const unsigned int m = 16807;
+const unsigned int b = 0;
+static long seed = 43690;//0xAAAA;
+// generate an uniformaly distributed random number sampled from [0;0x7FFFFFFF]
+long uniform() {
+ long unsigned int hi = m * (seed & n);
+ long unsigned int lo = m * (seed >> 16);
+
+ lo+= (hi & 32767) << 16;
+ lo+= hi >> 15;
+
+ if (lo > 2147483647) //0x7FFFFFFF)
+ lo -= 2147483647; //0x7FFFFFFF;
+ seed=(long) lo;
+ return seed ;
+}
+
+
+//--------------------------------------------------------------------------------------------------
+
+// RANDOM
+
+// generate an uniformaly distributed random number sampled from [-1;1]
+inline number randu()
+{
+ //on le veut entre -1 et 1 / pour l'instant on l'a entre 0 et 1
+ number result = uniform()/2147483647; //(double)0x7FFFFFFF;
+ // number result = (uniform() - (number)3xFFFFFFF.8)/(number)3xFFFFFFF.8);
+ return result;
+ // throw std::logic_error("Function not yet implemented.");
+}
+
+// generates a vector of the given size, the entries of the vector are all complex number
+// their real and imaginary parts are uniformaly distributed random number sampled from [-1;1]
+ComplexVector randu(int size)
+{
+ ComplexVector complexVector(size);
+ std::complex<number> mycomplex;
+ for(int i = 0; i < size; i++)
+ {
+ std::complex<number> mycomplex(randu(), randu());
+ complexVector(i) = mycomplex;
+ }
+ return complexVector;
+
+ //throw std::logic_error("Function not yet implemented.");
+}
+
+// generate a normally distributed random number
+// sampled from a distribution with mean 0 and std 1
+inline number randn()
+{
+
+ number res;
+ number U0 = uniform()/2147483647;//(double)0x7FFFFFFF;
+ number U1 = uniform()/2147483647;//(double)0x7FFFFFFF;
+ res = sqrt(-2 * log(U0)) * cos(2 * M_PI * U1);
+
+ return res;
+
+ // throw std::logic_error("Function not yet implemented.");
+}
+
+// generates a vector of the given size, the entries of the vector are all complex number
+// their real and imaginary parts are sampled from a normal distribution with mean 0 and std 1
+ComplexVector randn(int size)
+{
+ // NOTE: this can be built on the previous function
+ ComplexVector complexVector(size);
+ std::complex<number> mycomplex;
+ for(int i = 0; i < size; i++)
+ {
+ std::complex<number> mycomplex(randn(), randn());
+ complexVector(i) = mycomplex;
+ }
+ return complexVector;
+ //throw std::logic_error("Function not yet implemented.");
+}
+
+
+
+//-------------------------------------------------------------------------------------------------
+// IMPORT
+
+/// import a vector of reals from a file
+Vector loadRealVector(const std::string& file)
+{
+ std::vector<number> result;
+ std::ifstream input(file);
+
+ number x;
+ while (input >> x)
+ {
+ result.push_back(x);
+ }
+
+ Vector output = Eigen::Map<Vector, Eigen::Unaligned>(result.data(), result.size());
+ return output;
+}
+
+/// import a vector of complex from a file
+ComplexVector loadComplexVector(const std::string& file)
+{
+ std::vector<complexNumber> result;
+ std::ifstream input(file);
+
+ number r;
+ number i;
+ while (input >> r)
+ {
+ input >> i;
+ complexNumber x(r,i);
+ result.push_back(x);
+ }
+
+ ComplexVector output = Eigen::Map<ComplexVector, Eigen::Unaligned>(result.data(), result.size());
+ return output;
+}
+
+//-------------------------------------------------------------------------------------------------