diff --git a/matlab.h b/matlab.h deleted file mode 100644 index 27d163c6a99d83a3f11d77eae03367db11824da2..0000000000000000000000000000000000000000 --- a/matlab.h +++ /dev/null @@ -1,241 +0,0 @@ -#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; -} - -//-------------------------------------------------------------------------------------------------