view maths/pca/pca.c @ 321:f1e6be2de9a5

A threshold (delta) is added in the peak picking parameters structure (PPickParams). It is used as an offset when computing the smoothed detection function. A constructor for the structure PPickParams is also added to set the parameters to 0 when a structure instance is created. Hence programmes using the peak picking parameter structure and which do not set the delta parameter (e.g. QM Vamp note onset detector) won't be affected by the modifications.

Functions modified:
- dsp/onsets/PeakPicking.cpp
- dsp/onsets/PeakPicking.h
- dsp/signalconditioning/DFProcess.cpp
- dsp/signalconditioning/DFProcess.h

/*********************************/
/* Principal Components Analysis */
/*********************************/

/*********************************************************************/
/* Principal Components Analysis or the Karhunen-Loeve expansion is a
   classical method for dimensionality reduction or exploratory data
   analysis.  One reference among many is: F. Murtagh and A. Heck,
   Multivariate Data Analysis, Kluwer Academic, Dordrecht, 1987.

   Author:
   F. Murtagh
   Phone:        + 49 89 32006298 (work)
                 + 49 89 965307 (home)
   Earn/Bitnet:  fionn@dgaeso51,  fim@dgaipp1s,  murtagh@stsci
   Span:         esomc1::fionn
   Internet:     [e-mail]
   
   F. Murtagh, Munich, 6 June 1989                                   */   
/*********************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "pca.h"

#define SIGN(a, b) ( (b) < 0 ? -fabs(a) : fabs(a) )

/**  Variance-covariance matrix: creation  *****************************/

/* Create m * m covariance matrix from given n * m data matrix. */
void covcol(double** data, int n, int m, double** symmat)
{
double *mean;
int i, j, j1, j2;

/* Allocate storage for mean vector */

mean = (double*) malloc(m*sizeof(double));

/* Determine mean of column vectors of input data matrix */

for (j = 0; j < m; j++)
    {
    mean[j] = 0.0;
    for (i = 0; i < n; i++)
        {
        mean[j] += data[i][j];
        }
    mean[j] /= (double)n;
    }

/*
printf("\nMeans of column vectors:\n");
for (j = 0; j < m; j++)  {
    printf("%12.1f",mean[j]);  }   printf("\n");
 */

/* Center the column vectors. */

for (i = 0; i < n; i++)
    {
    for (j = 0; j < m; j++)
        {
        data[i][j] -= mean[j];
        }
    }

/* Calculate the m * m covariance matrix. */
for (j1 = 0; j1 < m; j1++)
    {
    for (j2 = j1; j2 < m; j2++)
        {
        symmat[j1][j2] = 0.0;
        for (i = 0; i < n; i++)
            {
            symmat[j1][j2] += data[i][j1] * data[i][j2];
            }
        symmat[j2][j1] = symmat[j1][j2];
        }
    }

free(mean);

return;

}

/**  Error handler  **************************************************/

void erhand(char* err_msg)
{
    fprintf(stderr,"Run-time error:\n");
    fprintf(stderr,"%s\n", err_msg);
    fprintf(stderr,"Exiting to system.\n");
    exit(1);
}


/**  Reduce a real, symmetric matrix to a symmetric, tridiag. matrix. */

/* Householder reduction of matrix a to tridiagonal form.
Algorithm: Martin et al., Num. Math. 11, 181-195, 1968.
Ref: Smith et al., Matrix Eigensystem Routines -- EISPACK Guide
Springer-Verlag, 1976, pp. 489-494.
W H Press et al., Numerical Recipes in C, Cambridge U P,
1988, pp. 373-374.  */
void tred2(double** a, int n, double* d, double* e)
{
	int l, k, j, i;
	double scale, hh, h, g, f;
	
	for (i = n-1; i >= 1; i--)
    {
		l = i - 1;
		h = scale = 0.0;
		if (l > 0)
		{
			for (k = 0; k <= l; k++)
				scale += fabs(a[i][k]);
			if (scale == 0.0)
				e[i] = a[i][l];
			else
			{
				for (k = 0; k <= l; k++)
				{
					a[i][k] /= scale;
					h += a[i][k] * a[i][k];
				}
				f = a[i][l];
				g = f>0 ? -sqrt(h) : sqrt(h);
				e[i] = scale * g;
				h -= f * g;
				a[i][l] = f - g;
				f = 0.0;
				for (j = 0; j <= l; j++)
				{
					a[j][i] = a[i][j]/h;
					g = 0.0;
					for (k = 0; k <= j; k++)
						g += a[j][k] * a[i][k];
					for (k = j+1; k <= l; k++)
						g += a[k][j] * a[i][k];
					e[j] = g / h;
					f += e[j] * a[i][j];
				}
				hh = f / (h + h);
				for (j = 0; j <= l; j++)
				{
					f = a[i][j];
					e[j] = g = e[j] - hh * f;
					for (k = 0; k <= j; k++)
						a[j][k] -= (f * e[k] + g * a[i][k]);
				}
			}
		}
		else
			e[i] = a[i][l];
		d[i] = h;
    }
	d[0] = 0.0;
	e[0] = 0.0;
	for (i = 0; i < n; i++)
    {
		l = i - 1;
		if (d[i])
		{
			for (j = 0; j <= l; j++)
			{
				g = 0.0;
				for (k = 0; k <= l; k++)
					g += a[i][k] * a[k][j];
				for (k = 0; k <= l; k++)
					a[k][j] -= g * a[k][i];
			}
		}
		d[i] = a[i][i];
		a[i][i] = 1.0;
		for (j = 0; j <= l; j++)
			a[j][i] = a[i][j] = 0.0;
    }
}

/**  Tridiagonal QL algorithm -- Implicit  **********************/

void tqli(double* d, double* e, int n, double** z)
{
	int m, l, iter, i, k;
	double s, r, p, g, f, dd, c, b;
	
	for (i = 1; i < n; i++)
		e[i-1] = e[i];
	e[n-1] = 0.0;
	for (l = 0; l < n; l++)
    {
		iter = 0;
		do
		{
			for (m = l; m < n-1; m++)
			{
				dd = fabs(d[m]) + fabs(d[m+1]);
				if (fabs(e[m]) + dd == dd) break;
			}
			if (m != l)
			{
				if (iter++ == 30) erhand("No convergence in TLQI.");
				g = (d[l+1] - d[l]) / (2.0 * e[l]);
				r = sqrt((g * g) + 1.0);
				g = d[m] - d[l] + e[l] / (g + SIGN(r, g));
				s = c = 1.0;
				p = 0.0;
				for (i = m-1; i >= l; i--)
				{
					f = s * e[i];
					b = c * e[i];
					if (fabs(f) >= fabs(g))
                    {
						c = g / f;
						r = sqrt((c * c) + 1.0);
						e[i+1] = f * r;
						c *= (s = 1.0/r);
                    }
					else
                    {
						s = f / g;
						r = sqrt((s * s) + 1.0);
						e[i+1] = g * r;
						s *= (c = 1.0/r);
                    }
					g = d[i+1] - p;
					r = (d[i] - g) * s + 2.0 * c * b;
					p = s * r;
					d[i+1] = g + p;
					g = c * r - b;
					for (k = 0; k < n; k++)
					{
						f = z[k][i+1];
						z[k][i+1] = s * z[k][i] + c * f;
						z[k][i] = c * z[k][i] - s * f;
					}
				}
				d[l] = d[l] - p;
				e[l] = g;
				e[m] = 0.0;
			}
		}  while (m != l);
	}
}

/* In place projection onto basis vectors */
void pca_project(double** data, int n, int m, int ncomponents)
{
	int  i, j, k, k2;
	double  **symmat, **symmat2, *evals, *interm;
	
	//TODO: assert ncomponents < m
	
	symmat = (double**) malloc(m*sizeof(double*));
	for (i = 0; i < m; i++)
		symmat[i] = (double*) malloc(m*sizeof(double));
		
	covcol(data, n, m, symmat);
	
	/*********************************************************************
		Eigen-reduction
		**********************************************************************/
	
    /* Allocate storage for dummy and new vectors. */
    evals = (double*) malloc(m*sizeof(double));     /* Storage alloc. for vector of eigenvalues */
    interm = (double*) malloc(m*sizeof(double));    /* Storage alloc. for 'intermediate' vector */
    //MALLOC_ARRAY(symmat2,m,m,double);    
	//for (i = 0; i < m; i++) {
	//	for (j = 0; j < m; j++) {
	//		symmat2[i][j] = symmat[i][j]; /* Needed below for col. projections */
	//	}
	//}
    tred2(symmat, m, evals, interm);  /* Triangular decomposition */
tqli(evals, interm, m, symmat);   /* Reduction of sym. trid. matrix */
/* evals now contains the eigenvalues,
columns of symmat now contain the associated eigenvectors. */	

/*
	printf("\nEigenvalues:\n");
	for (j = m-1; j >= 0; j--) {
		printf("%18.5f\n", evals[j]); }
	printf("\n(Eigenvalues should be strictly positive; limited\n");
	printf("precision machine arithmetic may affect this.\n");
	printf("Eigenvalues are often expressed as cumulative\n");
	printf("percentages, representing the 'percentage variance\n");
	printf("explained' by the associated axis or principal component.)\n");
	
	printf("\nEigenvectors:\n");
	printf("(First three; their definition in terms of original vbes.)\n");
	for (j = 0; j < m; j++) {
		for (i = 1; i <= 3; i++)  {
			printf("%12.4f", symmat[j][m-i]);  }
		printf("\n");  }
 */

/* Form projections of row-points on prin. components. */
/* Store in 'data', overwriting original data. */
for (i = 0; i < n; i++) {
	for (j = 0; j < m; j++) {
		interm[j] = data[i][j]; }   /* data[i][j] will be overwritten */
        for (k = 0; k < ncomponents; k++) {
			data[i][k] = 0.0;
			for (k2 = 0; k2 < m; k2++) {
				data[i][k] += interm[k2] * symmat[k2][m-k-1]; }
        }
}

/*	
printf("\nProjections of row-points on first 3 prin. comps.:\n");
 for (i = 0; i < n; i++) {
	 for (j = 0; j < 3; j++)  {
		 printf("%12.4f", data[i][j]);  }
	 printf("\n");  }
 */

/* Form projections of col.-points on first three prin. components. */
/* Store in 'symmat2', overwriting what was stored in this. */
//for (j = 0; j < m; j++) {
//	 for (k = 0; k < m; k++) {
//		 interm[k] = symmat2[j][k]; }  /*symmat2[j][k] will be overwritten*/
//  for (i = 0; i < 3; i++) {
//	symmat2[j][i] = 0.0;
//		for (k2 = 0; k2 < m; k2++) {
//			symmat2[j][i] += interm[k2] * symmat[k2][m-i-1]; }
//		if (evals[m-i-1] > 0.0005)   /* Guard against zero eigenvalue */
//			symmat2[j][i] /= sqrt(evals[m-i-1]);   /* Rescale */
//		else
//			symmat2[j][i] = 0.0;    /* Standard kludge */
//    }
// }

/*
 printf("\nProjections of column-points on first 3 prin. comps.:\n");
 for (j = 0; j < m; j++) {
	 for (k = 0; k < 3; k++)  {
		 printf("%12.4f", symmat2[j][k]);  }
	 printf("\n");  }
	*/


for (i = 0; i < m; i++)
	free(symmat[i]);
free(symmat);
//FREE_ARRAY(symmat2,m);
free(evals);
free(interm);

}