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];
}
}
