Mimic Octave's Inv() function

Open discussion regarding features, bugs, issues, vendors, etc.

Mimic Octave's Inv() function

Postby ionone » Tue Jul 09, 2013 6:11 am

Hi

I'm trying to replicate the Inv() function in Octave with no luck.

First i tried inverse matrix functions i found on the web but they weren't precise enough.

Then i heard Octave was using Lapack function so i tried those but i don't have the same numbers.
i tried dgetrf/dgetri and dsytrf/dsytri as it is a symetric matrix but the numbers don't match with Octave's.

for example in Octave the first value is : 1.7338e+010 and with dgetrf/dgetri i have : 17351976742.876453
with dsytrf/dsytri i have : 17352023489.128361 so you see they differ quite a lot...

So what Octave is using when calling "Inv()" ? thanks for any help

Jeff

here is the code i used in Visual Studio C++:
Code: Select all
#include <cstdio>
#include "stdafx.h"

extern "C" {
      // LU decomoposition of a general matrix
      void __cdecl dgetrf_(int* M, int *N, double* A, int* lda, int* IPIV, int* INFO);

      // generate inverse of a matrix given its LU decomposition
      void __cdecl dgetri_(int* N, double* A, int* lda, int* IPIV, double* WORK, int* lwork, int* INFO);


   
       int __cdecl  dsytrf_(char *uplo, int *n, double *a, int *lda, int *ipiv, double *work, int *lwork, int *info);
       int __cdecl  dsytri_(char *uplo, int *n, double *a, int *lda, int *ipiv, double *work, int *info);

   }



   

   void inverse(double* A, int N)
   {
      int *IPIV = new int[N+1];
      int LWORK = N*N;
      double *WORK = new double[LWORK];
      int INFO;

      dgetrf_(&N,&N,A,&N,IPIV,&INFO);
      dgetri_(&N,A,&N,IPIV,WORK,&LWORK,&INFO);

      delete IPIV;
      delete WORK;
   }
   void inverse2 (double * matrix, int matrix_rank)
   {
      int info = 0;
      int lwork = matrix_rank;
      int * ivpv = new   int [matrix_rank];
      double * work = new   double [lwork];

      dsytrf_("U", & matrix_rank, matrix, & matrix_rank, ivpv, work, & lwork, & info);
      dsytri_("U", & matrix_rank, matrix, & matrix_rank, ivpv, work, & info);

      for (int i = 0; i <matrix_rank; i++)
      {
         for (int j = i +1; j <matrix_rank; j++)
            matrix [i * matrix_rank + j] = matrix [j * matrix_rank + i];
      }

      delete [] ivpv;
      delete [] work;
   }



here is the Octave's code i used
Code: Select all
x = [1/1024:1/1024:1];

p = 20
h = 128;
w = 2*h;
ov = 1;

if (size(x,2) == 1)
  x = x';  % Convert X from column to row
end

npts = length(x);

nhops = floor(npts/h);

% Pad x with zeros so that we can extract complete w-length windows
% from it
x = [zeros(1,(w-h)/2),x,zeros(1,(w-h/2))];

a = zeros(nhops, p+1);
g = zeros(nhops, 1);
if ov == 0
  e = zeros(1, npts);
else
  e = zeros(1, (nhops-1)*h+w);
end

pre = [1 -0.9];
x = filter(pre,1,x);

nhops = 3
hop = 3

  % Extract segment of signal
  xx = x((hop - 1)*h + [1:w]);
  % Apply hanning window
  wxx = xx .* hanning(w)';
  % Form autocorrelation (calculates *way* too many points)
  rxx = xcorr(wxx);
  % extract just the points we need (middle p+1 points)
  rxx = rxx(w+[0:p]);% rxx ok
  % Setup the normal equations
  R = toeplitz(rxx(1:p));
inv(R)


and here is the C++ code, simply the same Octave code but in C++ (I kept the math way to count arrays (for 1 to n not from 0 to n-1 as in C)

Code: Select all
void fenetrage(double * t, int n)
{
   // hanning
   for (int i = 1; i <= n; i++)
   {
      t[i] *= (0.5 + 0.5 * cos( M_PI * (n/2 - ((i-1)*n/double(n-1))) / (n / 2) ));
   }
}

Code: Select all
              int nsmp = 1024;
   double* xx = new double[nsmp+1];
   for (int i = 1; i <= nsmp; i++)
      xx[i] = i/double(nsmp);

   float alpha = -0.2;

              double** a;
   double* g;
   double* e;
   double* x;

   int p = 12;

   int h = 128;

   int w = 2*h;

   int ov = 1;

   int sizeax;
   int sizeay;
   int sizee;

   int npts = nsmp;

   int nhops = floor(npts/double(h));

   //Pad x with zeros so that we can extract complete w-length windows
   // from it
   double* x2 = new double[nsmp+1];

   for (int i = 0; i <= nsmp; i++)
      x2[i] = xx[i];

   x = new double[(w-h)+nsmp+1];



   //x = [zeros(1,(w-h)/2),x,zeros(1,(w-h/2))];
   for (int i = 1; i <= (w-h)/2; i++)
      x[i] = 0;
   for (int i = (w-h)/2+1; i <= (w-h)/2+nsmp; i++)
   {
      x[i] = x2[i-(w-h)/2];
   }
   for (int i = (w-h)/2+nsmp+1; i <= (w-h)+nsmp; i++)
      x[i] = 0;

   
   
   sizeax = nhops;
              sizeay = p+1;
   //a = zeros(nhops, p+1);
   //g = zeros(nhops, 1);
   a = new double*[nhops+1];
   for (int i = 1; i <= nhops; i++)
   {
      a[i] = new double[(p+1)+1];
      for (int j = 1; j <= p+1; j++)
      {
         a[i][j] = 0;
      }
   }
   g = new double[nhops+1];
   for (int i = 1; i <= nhops; i++)
   {
      g[i] = 0;
   }

   int n1 = (nhops - 1)*h;
   if (ov == 0)
   {
      //e = zeros(1, npts);
      e = new double[npts+1];
      sizee = npts;
      for (int i = 1; i <= npts; i++)
         e[i] = 0;
   }
   else
   {
      //e = zeros(1, (nhops-1)*h+w);
      int n = (nhops-1)*h+w;
      e = new double[n+n1+1]; //on ajoute n1 smp car plus bas on en as besoin
      sizee = n;
      for (int i = 1; i <= n+n1; i++)
         e[i] = 0;
   }
   
   double * pre = new double[3];
   pre[1] = 1;
   pre[2] = -0.9;

   double prevX = 0;
   double xx2 = 0;
   for (int i = 1; i <= (w-h)+nsmp; i++)
   {
      xx2 = x[i];
      x[i] = pre[1] * xx2 + pre[2]*prevX;
      prevX = xx2;
   }

   free(pre);int n = w;
   double* rxx = new double[n*2+1];
   double* wxx = new double[w+1];double* rs = new double[w+1];
   //effacer
   nhops = 3;
   for (int hop = 3/*1*/; hop <= nhops/*nhops*/; hop++)
   {
     // Extract segment of signal
     //xx = x((hop - 1)*h + [1,2,3,4,...,w]);
     double* xx = new double[w+1];
     int nn = (hop - 1)*h;
     for (int i = 1; i <= w; i++)
     {
      xx[i] = x[nn + i];
     }

    

     // Apply hanning window
     //wxx = xx .* hanning(w)';
    
     for (int i = 1; i <= w; i++)
     {
      wxx[i] = xx[i];
     }
     fenetrage(wxx, w);//not CPU friendly
    
     double total = 0;
   
     for (int i = 1; i <= n; i++)
     {
        /*n-1 à n-1 * x[0]
        n-2 a n-1 * x[0 à 1]
        ...
        0    à n-1 * x[0 à n-1]*/
       
        total = 0;
        for (int j = n+1-i; j <= n; j++)
        {
            total = total + wxx[j]*wxx[j-(n+1-i)+1];
        }
       
        rxx[i] = total;
     }
 
    
     for (int i = 1; i < n; i++)
     {
      total = 0;
      for (int j = 1; j <= n - i; j++)
      {
         //i = 0    j = 0 à n-2    1 à n-1
         //i = n-2  j = 0 à 0          n-1
         total = total + wxx[j]* wxx[j + i];
      }
      rxx[i+n] = total;
     }   
         
     // extract just the points we need (middle p+1 points)
     //rxx = rxx(w+[0:p]);
     for (int i = 0; i <= p; i++)
     {
      rxx[i+1]= rxx[w+i];
     }
      
     //Setup the normal equations
     /*R = toeplitz(rxx(1:p));*/
     double** R = new double*[p+1];
     for (int i = 0; i <= p; i++)
     {
      R[i] = new double[p+1];
     }
     for (int i = p; i > 0; i--)
     {
        for (int j = 1; j <= i ; j++)
        {
           R[p-i+1][j+p-i] = rxx[j];
           R[j+p-i][p-i+1] = rxx[j];
        }
     }

     double** R2 = new double*[p];
     for (int i = 0; i < p; i++)
     {
      R2[i] = new double[p];
     }
     for (int i = 0; i < p; i++)
     {
        for (int j = 0; j < p ; j++)
        {
           R2[i][j] = R[i+1][j+1];
        }
     }
     double* R3 = new double[p*p];
     for (int i = 0; i < p; i++)
     {
        for (int j = 0; j < p ; j++)
        {
           R3[i*p+j] = R[i+1][j+1];
        }
      
     }
inverse(R3, p);
or : inverse2(R3,p);


thanks for any help

Jeffi
ionone
 
Posts: 2
Joined: Tue Jul 09, 2013 5:46 am

Return to User Discussion

Who is online

Users browsing this forum: No registered users and 1 guest

cron