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
