Mimic Octave's Inv() function

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

Mimic Octave's Inv() function

Postby ionone » Tue Jul 09, 2013 5:54 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.

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
aa = [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++
Code: Select all
int nsmp = 1024;
   double* d = new double[nsmp+1];
   for (int i = 1; i <= nsmp; i++)
      d[i] = i/double(nsmp);
   float alpha = -0.2;

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

   if (p == -1)
     p = 12;

   if (h == -1)
     h = 128;

   if (w == -1)
     w = 2*h;

   if (ov == -1)
     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;
   }
   
   // Pre-emphasis
   //x = filter(pre,1,x);
   //y(n) = 1 * x(n) - 0.9 * x(n-1)
   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
    
     //free(xx);


     // Form autocorrelation (calculates *way* too many points)
     //rxx = xcorr(wxx);
     //début:
     // wxx =     [1 2 3  4]
     //     [1 2 3 4]
     // fin:
     // wxx = [1 2 3 4]
     //             [1 2 3 4]
    
     //int an = n*2-1;
    
     //example
     //wxx =     [1 2 3 4]           4
     //    [1 2 3 4]
     //wxx =     [1 2 3 4]           11
     //      [1 2 3 4]
     //wxx =     [1 2 3 4]           20
     //        [1 2 3 4]
     //wxx =     [1 2 3 4]           30
     //          [1 2 3 4]
     //wxx =     [1 2 3 4]           20
     //            [1 2 3 4]
     //wxx =     [1 2 3 4]           11
     //              [1 2 3 4]
     //wxx =     [1 2 3 4]           4
     //                [1 2 3 4]


    
     double total = 0;
     /*for (int i = 0; i < n; i++)
     {
        / n-1 à n-1 * x[0]
        | n-2 a n-1 * x[0 à 1]
        | ...
        \ 0    à n-1 * x[0 à n-1]
       
        / 0    à n-2 * x[1 à n-1]
        | ...
        \ 0    à 0    * x[n-1]
     }*/
   
     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++)
     {
    /*   / 0    à n-2 * x[1 à n-1]
        | ...
        \ 0    à 0    * x[n-1]         */
      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));
    
                a b c d e
                b a b c d
     T(a,b,c,d,e) = c b a b c
                d c b a b
                e d c b a*/
     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];
        }
     }


    //Solve for a (horribly inefficient to use full inv())
    // an = inv(R)*rxx(2:(p+1))';
     //invR (voir test inv matrice pour bon code)
     //exemple :
     // / 1 2 \ * [2, 5] = [8, 30]
     // \ 3 4 /
      
     //invDet ne marche qu'avec des arrays de 0 à <n et pas 1 à <=n
     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];
        }
      
     }
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 6 guests