mirror of
https://github.com/Ardour/ardour.git
synced 2025-12-06 14:54:56 +01:00
remove more unused qm-dsp code (fixes windows compile no LAPACK)
This commit is contained in:
Robin Gareus
parent
a543ae329c
commit
344728551d
14 changed files with 1 additions and 8117 deletions
|
|
@ -1,398 +0,0 @@
|
|||
/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */
|
||||
|
||||
/*
|
||||
* ClusterMeltSegmenter.cpp
|
||||
*
|
||||
* Created by Mark Levy on 23/03/2006.
|
||||
* Copyright 2006 Centre for Digital Music, Queen Mary, University of London.
|
||||
|
||||
This program is free software; you can redistribute it and/or
|
||||
modify it under the terms of the GNU General Public License as
|
||||
published by the Free Software Foundation; either version 2 of the
|
||||
License, or (at your option) any later version. See the file
|
||||
COPYING included with this distribution for more information.
|
||||
*/
|
||||
|
||||
#include <cfloat>
|
||||
#include <cmath>
|
||||
|
||||
#include "ClusterMeltSegmenter.h"
|
||||
#include "cluster_segmenter.h"
|
||||
#include "segment.h"
|
||||
|
||||
#include "dsp/transforms/FFT.h"
|
||||
#include "dsp/chromagram/ConstantQ.h"
|
||||
#include "dsp/rateconversion/Decimator.h"
|
||||
#include "dsp/mfcc/MFCC.h"
|
||||
|
||||
ClusterMeltSegmenter::ClusterMeltSegmenter(ClusterMeltSegmenterParams params) :
|
||||
window(NULL),
|
||||
fft(NULL),
|
||||
constq(NULL),
|
||||
mfcc(NULL),
|
||||
featureType(params.featureType),
|
||||
hopSize(params.hopSize),
|
||||
windowSize(params.windowSize),
|
||||
fmin(params.fmin),
|
||||
fmax(params.fmax),
|
||||
nbins(params.nbins),
|
||||
ncomponents(params.ncomponents), // NB currently not passed - no. of PCA components is set in cluser_segmenter.c
|
||||
nHMMStates(params.nHMMStates),
|
||||
nclusters(params.nclusters),
|
||||
histogramLength(params.histogramLength),
|
||||
neighbourhoodLimit(params.neighbourhoodLimit),
|
||||
decimator(NULL)
|
||||
{
|
||||
}
|
||||
|
||||
void ClusterMeltSegmenter::initialise(int fs)
|
||||
{
|
||||
samplerate = fs;
|
||||
|
||||
if (featureType == FEATURE_TYPE_CONSTQ ||
|
||||
featureType == FEATURE_TYPE_CHROMA) {
|
||||
|
||||
// run internal processing at 11025 or thereabouts
|
||||
int internalRate = 11025;
|
||||
int decimationFactor = samplerate / internalRate;
|
||||
if (decimationFactor < 1) decimationFactor = 1;
|
||||
|
||||
// must be a power of two
|
||||
while (decimationFactor & (decimationFactor - 1)) ++decimationFactor;
|
||||
|
||||
if (decimationFactor > Decimator::getHighestSupportedFactor()) {
|
||||
decimationFactor = Decimator::getHighestSupportedFactor();
|
||||
}
|
||||
|
||||
if (decimationFactor > 1) {
|
||||
decimator = new Decimator(getWindowsize(), decimationFactor);
|
||||
}
|
||||
|
||||
CQConfig config;
|
||||
config.FS = samplerate / decimationFactor;
|
||||
config.min = fmin;
|
||||
config.max = fmax;
|
||||
config.BPO = nbins;
|
||||
config.CQThresh = 0.0054;
|
||||
|
||||
constq = new ConstantQ(config);
|
||||
constq->sparsekernel();
|
||||
|
||||
ncoeff = constq->getK();
|
||||
|
||||
fft = new FFTReal(constq->getfftlength());
|
||||
|
||||
} else if (featureType == FEATURE_TYPE_MFCC) {
|
||||
|
||||
// run internal processing at 22050 or thereabouts
|
||||
int internalRate = 22050;
|
||||
int decimationFactor = samplerate / internalRate;
|
||||
if (decimationFactor < 1) decimationFactor = 1;
|
||||
|
||||
// must be a power of two
|
||||
while (decimationFactor & (decimationFactor - 1)) ++decimationFactor;
|
||||
|
||||
if (decimationFactor > Decimator::getHighestSupportedFactor()) {
|
||||
decimationFactor = Decimator::getHighestSupportedFactor();
|
||||
}
|
||||
|
||||
if (decimationFactor > 1) {
|
||||
decimator = new Decimator(getWindowsize(), decimationFactor);
|
||||
}
|
||||
|
||||
MFCCConfig config(samplerate / decimationFactor);
|
||||
config.fftsize = 2048;
|
||||
config.nceps = 19;
|
||||
config.want_c0 = true;
|
||||
|
||||
mfcc = new MFCC(config);
|
||||
ncoeff = config.nceps + 1;
|
||||
}
|
||||
}
|
||||
|
||||
ClusterMeltSegmenter::~ClusterMeltSegmenter()
|
||||
{
|
||||
delete window;
|
||||
delete constq;
|
||||
delete decimator;
|
||||
delete fft;
|
||||
}
|
||||
|
||||
int
|
||||
ClusterMeltSegmenter::getWindowsize()
|
||||
{
|
||||
return static_cast<int>(windowSize * samplerate + 0.001);
|
||||
}
|
||||
|
||||
int
|
||||
ClusterMeltSegmenter::getHopsize()
|
||||
{
|
||||
return static_cast<int>(hopSize * samplerate + 0.001);
|
||||
}
|
||||
|
||||
void ClusterMeltSegmenter::extractFeatures(const double* samples, int nsamples)
|
||||
{
|
||||
if (featureType == FEATURE_TYPE_CONSTQ ||
|
||||
featureType == FEATURE_TYPE_CHROMA) {
|
||||
extractFeaturesConstQ(samples, nsamples);
|
||||
} else if (featureType == FEATURE_TYPE_MFCC) {
|
||||
extractFeaturesMFCC(samples, nsamples);
|
||||
}
|
||||
}
|
||||
|
||||
void ClusterMeltSegmenter::extractFeaturesConstQ(const double* samples, int nsamples)
|
||||
{
|
||||
if (!constq) {
|
||||
std::cerr << "ERROR: ClusterMeltSegmenter::extractFeaturesConstQ: "
|
||||
<< "No const-q: initialise not called?"
|
||||
<< std::endl;
|
||||
return;
|
||||
}
|
||||
|
||||
if (nsamples < getWindowsize()) {
|
||||
std::cerr << "ERROR: ClusterMeltSegmenter::extractFeatures: nsamples < windowsize (" << nsamples << " < " << getWindowsize() << ")" << std::endl;
|
||||
return;
|
||||
}
|
||||
|
||||
int fftsize = constq->getfftlength();
|
||||
|
||||
if (!window || window->getSize() != fftsize) {
|
||||
delete window;
|
||||
window = new Window<double>(HammingWindow, fftsize);
|
||||
}
|
||||
|
||||
vector<double> cq(ncoeff);
|
||||
|
||||
for (int i = 0; i < ncoeff; ++i) cq[i] = 0.0;
|
||||
|
||||
const double *psource = samples;
|
||||
int pcount = nsamples;
|
||||
|
||||
if (decimator) {
|
||||
pcount = nsamples / decimator->getFactor();
|
||||
double *decout = new double[pcount];
|
||||
decimator->process(samples, decout);
|
||||
psource = decout;
|
||||
}
|
||||
|
||||
int origin = 0;
|
||||
|
||||
// std::cerr << "nsamples = " << nsamples << ", pcount = " << pcount << std::endl;
|
||||
|
||||
int frames = 0;
|
||||
|
||||
double *frame = new double[fftsize];
|
||||
double *real = new double[fftsize];
|
||||
double *imag = new double[fftsize];
|
||||
double *cqre = new double[ncoeff];
|
||||
double *cqim = new double[ncoeff];
|
||||
|
||||
while (origin <= pcount) {
|
||||
|
||||
// always need at least one fft window per block, but after
|
||||
// that we want to avoid having any incomplete ones
|
||||
if (origin > 0 && origin + fftsize >= pcount) break;
|
||||
|
||||
for (int i = 0; i < fftsize; ++i) {
|
||||
if (origin + i < pcount) {
|
||||
frame[i] = psource[origin + i];
|
||||
} else {
|
||||
frame[i] = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
for (int i = 0; i < fftsize/2; ++i) {
|
||||
double value = frame[i];
|
||||
frame[i] = frame[i + fftsize/2];
|
||||
frame[i + fftsize/2] = value;
|
||||
}
|
||||
|
||||
window->cut(frame);
|
||||
|
||||
fft->forward(frame, real, imag);
|
||||
|
||||
constq->process(real, imag, cqre, cqim);
|
||||
|
||||
for (int i = 0; i < ncoeff; ++i) {
|
||||
cq[i] += sqrt(cqre[i] * cqre[i] + cqim[i] * cqim[i]);
|
||||
}
|
||||
++frames;
|
||||
|
||||
origin += fftsize/2;
|
||||
}
|
||||
|
||||
delete [] cqre;
|
||||
delete [] cqim;
|
||||
delete [] real;
|
||||
delete [] imag;
|
||||
delete [] frame;
|
||||
|
||||
for (int i = 0; i < ncoeff; ++i) {
|
||||
cq[i] /= frames;
|
||||
}
|
||||
|
||||
if (decimator) delete[] psource;
|
||||
|
||||
features.push_back(cq);
|
||||
}
|
||||
|
||||
void ClusterMeltSegmenter::extractFeaturesMFCC(const double* samples, int nsamples)
|
||||
{
|
||||
if (!mfcc) {
|
||||
std::cerr << "ERROR: ClusterMeltSegmenter::extractFeaturesMFCC: "
|
||||
<< "No mfcc: initialise not called?"
|
||||
<< std::endl;
|
||||
return;
|
||||
}
|
||||
|
||||
if (nsamples < getWindowsize()) {
|
||||
std::cerr << "ERROR: ClusterMeltSegmenter::extractFeatures: nsamples < windowsize (" << nsamples << " < " << getWindowsize() << ")" << std::endl;
|
||||
return;
|
||||
}
|
||||
|
||||
int fftsize = mfcc->getfftlength();
|
||||
|
||||
vector<double> cc(ncoeff);
|
||||
|
||||
for (int i = 0; i < ncoeff; ++i) cc[i] = 0.0;
|
||||
|
||||
const double *psource = samples;
|
||||
int pcount = nsamples;
|
||||
|
||||
if (decimator) {
|
||||
pcount = nsamples / decimator->getFactor();
|
||||
double *decout = new double[pcount];
|
||||
decimator->process(samples, decout);
|
||||
psource = decout;
|
||||
}
|
||||
|
||||
int origin = 0;
|
||||
int frames = 0;
|
||||
|
||||
double *frame = new double[fftsize];
|
||||
double *ccout = new double[ncoeff];
|
||||
|
||||
while (origin <= pcount) {
|
||||
|
||||
// always need at least one fft window per block, but after
|
||||
// that we want to avoid having any incomplete ones
|
||||
if (origin > 0 && origin + fftsize >= pcount) break;
|
||||
|
||||
for (int i = 0; i < fftsize; ++i) {
|
||||
if (origin + i < pcount) {
|
||||
frame[i] = psource[origin + i];
|
||||
} else {
|
||||
frame[i] = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
mfcc->process(frame, ccout);
|
||||
|
||||
for (int i = 0; i < ncoeff; ++i) {
|
||||
cc[i] += ccout[i];
|
||||
}
|
||||
++frames;
|
||||
|
||||
origin += fftsize/2;
|
||||
}
|
||||
|
||||
delete [] ccout;
|
||||
delete [] frame;
|
||||
|
||||
for (int i = 0; i < ncoeff; ++i) {
|
||||
cc[i] /= frames;
|
||||
}
|
||||
|
||||
if (decimator) delete[] psource;
|
||||
|
||||
features.push_back(cc);
|
||||
}
|
||||
|
||||
void ClusterMeltSegmenter::segment(int m)
|
||||
{
|
||||
nclusters = m;
|
||||
segment();
|
||||
}
|
||||
|
||||
void ClusterMeltSegmenter::setFeatures(const vector<vector<double> >& f)
|
||||
{
|
||||
features = f;
|
||||
featureType = FEATURE_TYPE_UNKNOWN;
|
||||
}
|
||||
|
||||
void ClusterMeltSegmenter::segment()
|
||||
{
|
||||
delete constq;
|
||||
constq = 0;
|
||||
delete mfcc;
|
||||
mfcc = 0;
|
||||
delete decimator;
|
||||
decimator = 0;
|
||||
|
||||
if (features.size() < histogramLength) return;
|
||||
/*
|
||||
std::cerr << "ClusterMeltSegmenter::segment: have " << features.size()
|
||||
<< " features with " << features[0].size() << " coefficients (ncoeff = " << ncoeff << ", ncomponents = " << ncomponents << ")" << std::endl;
|
||||
*/
|
||||
// copy the features to a native array and use the existing C segmenter...
|
||||
double** arrFeatures = new double*[features.size()];
|
||||
for (int i = 0; i < features.size(); i++)
|
||||
{
|
||||
if (featureType == FEATURE_TYPE_UNKNOWN) {
|
||||
arrFeatures[i] = new double[features[0].size()];
|
||||
for (int j = 0; j < features[0].size(); j++)
|
||||
arrFeatures[i][j] = features[i][j];
|
||||
} else {
|
||||
arrFeatures[i] = new double[ncoeff+1]; // allow space for the normalised envelope
|
||||
for (int j = 0; j < ncoeff; j++)
|
||||
arrFeatures[i][j] = features[i][j];
|
||||
}
|
||||
}
|
||||
|
||||
q = new int[features.size()];
|
||||
|
||||
if (featureType == FEATURE_TYPE_UNKNOWN ||
|
||||
featureType == FEATURE_TYPE_MFCC)
|
||||
cluster_segment(q, arrFeatures, features.size(), features[0].size(), nHMMStates, histogramLength,
|
||||
nclusters, neighbourhoodLimit);
|
||||
else
|
||||
constq_segment(q, arrFeatures, features.size(), nbins, ncoeff, featureType,
|
||||
nHMMStates, histogramLength, nclusters, neighbourhoodLimit);
|
||||
|
||||
// convert the cluster assignment sequence to a segmentation
|
||||
makeSegmentation(q, features.size());
|
||||
|
||||
// de-allocate arrays
|
||||
delete [] q;
|
||||
for (int i = 0; i < features.size(); i++)
|
||||
delete [] arrFeatures[i];
|
||||
delete [] arrFeatures;
|
||||
|
||||
// clear the features
|
||||
clear();
|
||||
}
|
||||
|
||||
void ClusterMeltSegmenter::makeSegmentation(int* q, int len)
|
||||
{
|
||||
segmentation.segments.clear();
|
||||
segmentation.nsegtypes = nclusters;
|
||||
segmentation.samplerate = samplerate;
|
||||
|
||||
Segment segment;
|
||||
segment.start = 0;
|
||||
segment.type = q[0];
|
||||
|
||||
for (int i = 1; i < len; i++)
|
||||
{
|
||||
if (q[i] != q[i-1])
|
||||
{
|
||||
segment.end = i * getHopsize();
|
||||
segmentation.segments.push_back(segment);
|
||||
segment.type = q[i];
|
||||
segment.start = segment.end;
|
||||
}
|
||||
}
|
||||
segment.end = len * getHopsize();
|
||||
segmentation.segments.push_back(segment);
|
||||
}
|
||||
|
||||
|
|
@ -1,109 +0,0 @@
|
|||
/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */
|
||||
|
||||
/*
|
||||
* ClusterMeltSegmenter.h
|
||||
*
|
||||
* Created by Mark Levy on 23/03/2006.
|
||||
* Copyright 2006 Centre for Digital Music, Queen Mary, University of London.
|
||||
|
||||
This program is free software; you can redistribute it and/or
|
||||
modify it under the terms of the GNU General Public License as
|
||||
published by the Free Software Foundation; either version 2 of the
|
||||
License, or (at your option) any later version. See the file
|
||||
COPYING included with this distribution for more information.
|
||||
*/
|
||||
|
||||
#include <vector>
|
||||
|
||||
#include "segment.h"
|
||||
#include "Segmenter.h"
|
||||
#include "hmm/hmm.h"
|
||||
#include "base/Window.h"
|
||||
|
||||
using std::vector;
|
||||
|
||||
class Decimator;
|
||||
class ConstantQ;
|
||||
class MFCC;
|
||||
class FFTReal;
|
||||
|
||||
class ClusterMeltSegmenterParams
|
||||
// defaults are sensible for 11025Hz with 0.2 second hopsize
|
||||
{
|
||||
public:
|
||||
ClusterMeltSegmenterParams() :
|
||||
featureType(FEATURE_TYPE_CONSTQ),
|
||||
hopSize(0.2),
|
||||
windowSize(0.6),
|
||||
fmin(62),
|
||||
fmax(16000),
|
||||
nbins(8),
|
||||
ncomponents(20),
|
||||
nHMMStates(40),
|
||||
nclusters(10),
|
||||
histogramLength(15),
|
||||
neighbourhoodLimit(20) { }
|
||||
feature_types featureType;
|
||||
double hopSize; // in secs
|
||||
double windowSize; // in secs
|
||||
int fmin;
|
||||
int fmax;
|
||||
int nbins;
|
||||
int ncomponents;
|
||||
int nHMMStates;
|
||||
int nclusters;
|
||||
int histogramLength;
|
||||
int neighbourhoodLimit;
|
||||
};
|
||||
|
||||
class ClusterMeltSegmenter : public Segmenter
|
||||
{
|
||||
public:
|
||||
ClusterMeltSegmenter(ClusterMeltSegmenterParams params);
|
||||
virtual ~ClusterMeltSegmenter();
|
||||
virtual void initialise(int samplerate);
|
||||
virtual int getWindowsize();
|
||||
virtual int getHopsize();
|
||||
virtual void extractFeatures(const double* samples, int nsamples);
|
||||
void setFeatures(const vector<vector<double> >& f); // provide the features yourself
|
||||
virtual void segment(); // segment into default number of segment-types
|
||||
void segment(int m); // segment into m segment-types
|
||||
int getNSegmentTypes() { return nclusters; }
|
||||
|
||||
protected:
|
||||
void makeSegmentation(int* q, int len);
|
||||
|
||||
void extractFeaturesConstQ(const double *, int);
|
||||
void extractFeaturesMFCC(const double *, int);
|
||||
|
||||
Window<double> *window;
|
||||
FFTReal *fft;
|
||||
ConstantQ* constq;
|
||||
MFCC* mfcc;
|
||||
model_t* model; // the HMM
|
||||
int* q; // the decoded HMM state sequence
|
||||
vector<vector<double> > histograms;
|
||||
|
||||
feature_types featureType;
|
||||
double hopSize; // in seconds
|
||||
double windowSize; // in seconds
|
||||
|
||||
// constant-Q parameters
|
||||
int fmin;
|
||||
int fmax;
|
||||
int nbins;
|
||||
int ncoeff;
|
||||
|
||||
// PCA parameters
|
||||
int ncomponents;
|
||||
|
||||
// HMM parameters
|
||||
int nHMMStates;
|
||||
|
||||
// clustering parameters
|
||||
int nclusters;
|
||||
int histogramLength;
|
||||
int neighbourhoodLimit;
|
||||
|
||||
Decimator *decimator;
|
||||
};
|
||||
|
|
@ -1,31 +0,0 @@
|
|||
/*
|
||||
* Segmenter.cpp
|
||||
*
|
||||
* Created by Mark Levy on 04/04/2006.
|
||||
* Copyright 2006 Centre for Digital Music, Queen Mary, University of London.
|
||||
|
||||
This program is free software; you can redistribute it and/or
|
||||
modify it under the terms of the GNU General Public License as
|
||||
published by the Free Software Foundation; either version 2 of the
|
||||
License, or (at your option) any later version. See the file
|
||||
COPYING included with this distribution for more information.
|
||||
*
|
||||
*/
|
||||
|
||||
#include <iomanip>
|
||||
|
||||
#include "Segmenter.h"
|
||||
|
||||
ostream& operator<<(ostream& os, const Segmentation& s)
|
||||
{
|
||||
os << "structure_name : begin_time end_time\n";
|
||||
|
||||
for (int i = 0; i < s.segments.size(); i++)
|
||||
{
|
||||
Segment seg = s.segments[i];
|
||||
os << std::fixed << seg.type << ':' << '\t' << std::setprecision(6) << seg.start / static_cast<double>(s.samplerate)
|
||||
<< '\t' << std::setprecision(6) << seg.end / static_cast<double>(s.samplerate) << "\n";
|
||||
}
|
||||
|
||||
return os;
|
||||
}
|
||||
|
|
@ -1,62 +0,0 @@
|
|||
#ifndef _SEGMENTER_H
|
||||
#define _SEGMENTER_H
|
||||
|
||||
/*
|
||||
* Segmenter.h
|
||||
* soundbite
|
||||
*
|
||||
* Created by Mark Levy on 23/03/2006.
|
||||
* Copyright 2006 Centre for Digital Music, Queen Mary, University of London.
|
||||
|
||||
This program is free software; you can redistribute it and/or
|
||||
modify it under the terms of the GNU General Public License as
|
||||
published by the Free Software Foundation; either version 2 of the
|
||||
License, or (at your option) any later version. See the file
|
||||
COPYING included with this distribution for more information.
|
||||
*
|
||||
*/
|
||||
|
||||
#include <vector>
|
||||
#include <iostream>
|
||||
|
||||
using std::vector;
|
||||
using std::ostream;
|
||||
|
||||
class Segment
|
||||
{
|
||||
public:
|
||||
int start; // in samples
|
||||
int end;
|
||||
int type;
|
||||
};
|
||||
|
||||
class Segmentation
|
||||
{
|
||||
public:
|
||||
int nsegtypes; // number of segment types, so possible types are {0,1,...,nsegtypes-1}
|
||||
int samplerate;
|
||||
vector<Segment> segments;
|
||||
};
|
||||
|
||||
ostream& operator<<(ostream& os, const Segmentation& s);
|
||||
|
||||
class Segmenter
|
||||
{
|
||||
public:
|
||||
Segmenter() {}
|
||||
virtual ~Segmenter() {}
|
||||
virtual void initialise(int samplerate) = 0; // must be called before any other methods
|
||||
virtual int getWindowsize() = 0; // required window size for calls to extractFeatures()
|
||||
virtual int getHopsize() = 0; // required hop size for calls to extractFeatures()
|
||||
virtual void extractFeatures(const double* samples, int nsamples) = 0;
|
||||
virtual void segment() = 0; // call once all the features have been extracted
|
||||
virtual void segment(int m) = 0; // specify desired number of segment-types
|
||||
virtual void clear() { features.clear(); }
|
||||
const Segmentation& getSegmentation() const { return segmentation; }
|
||||
protected:
|
||||
vector<vector<double> > features;
|
||||
Segmentation segmentation;
|
||||
int samplerate;
|
||||
};
|
||||
|
||||
#endif
|
||||
|
|
@ -1,225 +0,0 @@
|
|||
/*
|
||||
* cluster.c
|
||||
* cluster_melt
|
||||
*
|
||||
* Created by Mark Levy on 21/02/2006.
|
||||
* Copyright 2006 Centre for Digital Music, Queen Mary, University of London.
|
||||
|
||||
This program is free software; you can redistribute it and/or
|
||||
modify it under the terms of the GNU General Public License as
|
||||
published by the Free Software Foundation; either version 2 of the
|
||||
License, or (at your option) any later version. See the file
|
||||
COPYING included with this distribution for more information.
|
||||
*
|
||||
*/
|
||||
|
||||
#include <stdlib.h>
|
||||
|
||||
#include "cluster_melt.h"
|
||||
|
||||
#define DEFAULT_LAMBDA 0.02;
|
||||
#define DEFAULT_LIMIT 20;
|
||||
|
||||
double kldist(double* a, double* b, int n) {
|
||||
/* NB assume that all a[i], b[i] are non-negative
|
||||
because a, b represent probability distributions */
|
||||
double q, d;
|
||||
int i;
|
||||
|
||||
d = 0;
|
||||
for (i = 0; i < n; i++)
|
||||
{
|
||||
q = (a[i] + b[i]) / 2.0;
|
||||
if (q > 0)
|
||||
{
|
||||
if (a[i] > 0)
|
||||
d += a[i] * log(a[i] / q);
|
||||
if (b[i] > 0)
|
||||
d += b[i] * log(b[i] / q);
|
||||
}
|
||||
}
|
||||
return d;
|
||||
}
|
||||
|
||||
void cluster_melt(double *h, int m, int n, double *Bsched, int t, int k, int l, int *c) {
|
||||
double lambda, sum, beta, logsumexp, maxlp;
|
||||
int i, j, a, b, b0, b1, limit, B, it, maxiter, maxiter0, maxiter1;
|
||||
double** cl; /* reference histograms for each cluster */
|
||||
int** nc; /* neighbour counts for each histogram */
|
||||
double** lp; /* soft assignment probs for each histogram */
|
||||
int* oldc; /* previous hard assignments (to check convergence) */
|
||||
|
||||
/* NB h is passed as a 1d row major array */
|
||||
|
||||
/* parameter values */
|
||||
lambda = DEFAULT_LAMBDA;
|
||||
if (l > 0)
|
||||
limit = l;
|
||||
else
|
||||
limit = DEFAULT_LIMIT; /* use default if no valid neighbourhood limit supplied */
|
||||
B = 2 * limit + 1;
|
||||
maxiter0 = 20; /* number of iterations at initial temperature */
|
||||
maxiter1 = 5; /* number of iterations at subsequent temperatures */
|
||||
|
||||
/* allocate memory */
|
||||
cl = (double**) malloc(k*sizeof(double*));
|
||||
for (i= 0; i < k; i++)
|
||||
cl[i] = (double*) malloc(m*sizeof(double));
|
||||
|
||||
nc = (int**) malloc(n*sizeof(int*));
|
||||
for (i= 0; i < n; i++)
|
||||
nc[i] = (int*) malloc(k*sizeof(int));
|
||||
|
||||
lp = (double**) malloc(n*sizeof(double*));
|
||||
for (i= 0; i < n; i++)
|
||||
lp[i] = (double*) malloc(k*sizeof(double));
|
||||
|
||||
oldc = (int*) malloc(n * sizeof(int));
|
||||
|
||||
/* initialise */
|
||||
for (i = 0; i < k; i++)
|
||||
{
|
||||
sum = 0;
|
||||
for (j = 0; j < m; j++)
|
||||
{
|
||||
cl[i][j] = rand(); /* random initial reference histograms */
|
||||
sum += cl[i][j] * cl[i][j];
|
||||
}
|
||||
sum = sqrt(sum);
|
||||
for (j = 0; j < m; j++)
|
||||
{
|
||||
cl[i][j] /= sum; /* normalise */
|
||||
}
|
||||
}
|
||||
//print_array(cl, k, m);
|
||||
|
||||
for (i = 0; i < n; i++)
|
||||
c[i] = 1; /* initially assign all histograms to cluster 1 */
|
||||
|
||||
for (a = 0; a < t; a++)
|
||||
{
|
||||
beta = Bsched[a];
|
||||
|
||||
if (a == 0)
|
||||
maxiter = maxiter0;
|
||||
else
|
||||
maxiter = maxiter1;
|
||||
|
||||
for (it = 0; it < maxiter; it++)
|
||||
{
|
||||
//if (it == maxiter - 1)
|
||||
// mexPrintf("hasn't converged after %d iterations\n", maxiter);
|
||||
|
||||
for (i = 0; i < n; i++)
|
||||
{
|
||||
/* save current hard assignments */
|
||||
oldc[i] = c[i];
|
||||
|
||||
/* calculate soft assignment logprobs for each cluster */
|
||||
sum = 0;
|
||||
for (j = 0; j < k; j++)
|
||||
{
|
||||
lp[i][ j] = -beta * kldist(cl[j], &h[i*m], m);
|
||||
|
||||
/* update matching neighbour counts for this histogram, based on current hard assignments */
|
||||
/* old version:
|
||||
nc[i][j] = 0;
|
||||
if (i >= limit && i <= n - 1 - limit)
|
||||
{
|
||||
for (b = i - limit; b <= i + limit; b++)
|
||||
{
|
||||
if (c[b] == j+1)
|
||||
nc[i][j]++;
|
||||
}
|
||||
nc[i][j] = B - nc[i][j];
|
||||
}
|
||||
*/
|
||||
b0 = i - limit;
|
||||
if (b0 < 0)
|
||||
b0 = 0;
|
||||
b1 = i + limit;
|
||||
if (b1 >= n)
|
||||
b1 = n - 1;
|
||||
nc[i][j] = b1 - b0 + 1; /* = B except at edges */
|
||||
for (b = b0; b <= b1; b++)
|
||||
if (c[b] == j+1)
|
||||
nc[i][j]--;
|
||||
|
||||
sum += exp(lp[i][j]);
|
||||
}
|
||||
|
||||
/* normalise responsibilities and add duration logprior */
|
||||
logsumexp = log(sum);
|
||||
for (j = 0; j < k; j++)
|
||||
lp[i][j] -= logsumexp + lambda * nc[i][j];
|
||||
}
|
||||
//print_array(lp, n, k);
|
||||
/*
|
||||
for (i = 0; i < n; i++)
|
||||
{
|
||||
for (j = 0; j < k; j++)
|
||||
mexPrintf("%d ", nc[i][j]);
|
||||
mexPrintf("\n");
|
||||
}
|
||||
*/
|
||||
|
||||
|
||||
/* update the assignments now that we know the duration priors
|
||||
based on the current assignments */
|
||||
for (i = 0; i < n; i++)
|
||||
{
|
||||
maxlp = lp[i][0];
|
||||
c[i] = 1;
|
||||
for (j = 1; j < k; j++)
|
||||
if (lp[i][j] > maxlp)
|
||||
{
|
||||
maxlp = lp[i][j];
|
||||
c[i] = j+1;
|
||||
}
|
||||
}
|
||||
|
||||
/* break if assignments haven't changed */
|
||||
i = 0;
|
||||
while (i < n && oldc[i] == c[i])
|
||||
i++;
|
||||
if (i == n)
|
||||
break;
|
||||
|
||||
/* update reference histograms now we know new responsibilities */
|
||||
for (j = 0; j < k; j++)
|
||||
{
|
||||
for (b = 0; b < m; b++)
|
||||
{
|
||||
cl[j][b] = 0;
|
||||
for (i = 0; i < n; i++)
|
||||
{
|
||||
cl[j][b] += exp(lp[i][j]) * h[i*m+b];
|
||||
}
|
||||
}
|
||||
|
||||
sum = 0;
|
||||
for (i = 0; i < n; i++)
|
||||
sum += exp(lp[i][j]);
|
||||
for (b = 0; b < m; b++)
|
||||
cl[j][b] /= sum; /* normalise */
|
||||
}
|
||||
|
||||
//print_array(cl, k, m);
|
||||
//mexPrintf("\n\n");
|
||||
}
|
||||
}
|
||||
|
||||
/* free memory */
|
||||
for (i = 0; i < k; i++)
|
||||
free(cl[i]);
|
||||
free(cl);
|
||||
for (i = 0; i < n; i++)
|
||||
free(nc[i]);
|
||||
free(nc);
|
||||
for (i = 0; i < n; i++)
|
||||
free(lp[i]);
|
||||
free(lp);
|
||||
free(oldc);
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -1,39 +0,0 @@
|
|||
#ifndef _CLUSTER_MELT_H
|
||||
#define _CLUSTER_MELT_H
|
||||
/*
|
||||
* cluster_melt.h
|
||||
* cluster_melt
|
||||
*
|
||||
* Created by Mark Levy on 21/02/2006.
|
||||
* Copyright 2006 Centre for Digital Music, Queen Mary, University of London.
|
||||
|
||||
This program is free software; you can redistribute it and/or
|
||||
modify it under the terms of the GNU General Public License as
|
||||
published by the Free Software Foundation; either version 2 of the
|
||||
License, or (at your option) any later version. See the file
|
||||
COPYING included with this distribution for more information.
|
||||
*
|
||||
*/
|
||||
|
||||
#include <stdlib.h>
|
||||
#include <math.h>
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
|
||||
void cluster_melt(double *h, /* normalised histograms, as a vector in row major order */
|
||||
int m, /* number of dimensions (i.e. histogram bins) */
|
||||
int n, /* number of histograms */
|
||||
double *Bsched, /* inverse temperature schedule */
|
||||
int t, /* length of schedule */
|
||||
int k, /* number of clusters */
|
||||
int l, /* neighbourhood limit (supply zero to use default value) */
|
||||
int *c /* sequence of cluster assignments */
|
||||
);
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
||||
#endif
|
||||
|
|
@ -1,285 +0,0 @@
|
|||
/*
|
||||
* cluster_segmenter.c
|
||||
* soundbite
|
||||
*
|
||||
* Created by Mark Levy on 06/04/2006.
|
||||
* Copyright 2006 Centre for Digital Music, Queen Mary, University of London.
|
||||
|
||||
This program is free software; you can redistribute it and/or
|
||||
modify it under the terms of the GNU General Public License as
|
||||
published by the Free Software Foundation; either version 2 of the
|
||||
License, or (at your option) any later version. See the file
|
||||
COPYING included with this distribution for more information.
|
||||
*
|
||||
*/
|
||||
|
||||
#include "cluster_segmenter.h"
|
||||
|
||||
extern int readmatarray_size(const char *filepath, int n_array, int* t, int* d);
|
||||
extern int readmatarray(const char *filepath, int n_array, int t, int d, double** arr);
|
||||
|
||||
/* converts constant-Q features to normalised chroma */
|
||||
void cq2chroma(double** cq, int nframes, int ncoeff, int bins, double** chroma)
|
||||
{
|
||||
int noct = ncoeff / bins; /* number of complete octaves in constant-Q */
|
||||
int t, b, oct, ix;
|
||||
//double maxchroma; /* max chroma value at each time, for normalisation */
|
||||
//double sum; /* for normalisation */
|
||||
|
||||
for (t = 0; t < nframes; t++)
|
||||
{
|
||||
for (b = 0; b < bins; b++)
|
||||
chroma[t][b] = 0;
|
||||
for (oct = 0; oct < noct; oct++)
|
||||
{
|
||||
ix = oct * bins;
|
||||
for (b = 0; b < bins; b++)
|
||||
chroma[t][b] += fabs(cq[t][ix+b]);
|
||||
}
|
||||
/* normalise to unit sum
|
||||
sum = 0;
|
||||
for (b = 0; b < bins; b++)
|
||||
sum += chroma[t][b];
|
||||
for (b = 0; b < bins; b++)
|
||||
chroma[t][b] /= sum;
|
||||
*/
|
||||
/* normalise to unit max - NO this made results much worse!
|
||||
maxchroma = 0;
|
||||
for (b = 0; b < bins; b++)
|
||||
if (chroma[t][b] > maxchroma)
|
||||
maxchroma = chroma[t][b];
|
||||
if (maxchroma > 0)
|
||||
for (b = 0; b < bins; b++)
|
||||
chroma[t][b] /= maxchroma;
|
||||
*/
|
||||
}
|
||||
}
|
||||
|
||||
/* applies MPEG-7 normalisation to constant-Q features, storing normalised envelope (norm) in last feature dimension */
|
||||
void mpeg7_constq(double** features, int nframes, int ncoeff)
|
||||
{
|
||||
int i, j;
|
||||
double ss;
|
||||
double env;
|
||||
double maxenv = 0;
|
||||
|
||||
/* convert const-Q features to dB scale */
|
||||
for (i = 0; i < nframes; i++)
|
||||
for (j = 0; j < ncoeff; j++)
|
||||
features[i][j] = 10.0 * log10(features[i][j]+DBL_EPSILON);
|
||||
|
||||
/* normalise each feature vector and add the norm as an extra feature dimension */
|
||||
for (i = 0; i < nframes; i++)
|
||||
{
|
||||
ss = 0;
|
||||
for (j = 0; j < ncoeff; j++)
|
||||
ss += features[i][j] * features[i][j];
|
||||
env = sqrt(ss);
|
||||
for (j = 0; j < ncoeff; j++)
|
||||
features[i][j] /= env;
|
||||
features[i][ncoeff] = env;
|
||||
if (env > maxenv)
|
||||
maxenv = env;
|
||||
}
|
||||
/* normalise the envelopes */
|
||||
for (i = 0; i < nframes; i++)
|
||||
features[i][ncoeff] /= maxenv;
|
||||
}
|
||||
|
||||
/* return histograms h[nx*m] of data x[nx] into m bins using a sliding window of length h_len (MUST BE ODD) */
|
||||
/* NB h is a vector in row major order, as required by cluster_melt() */
|
||||
/* for historical reasons we normalise the histograms by their norm (not to sum to one) */
|
||||
void create_histograms(int* x, int nx, int m, int hlen, double* h)
|
||||
{
|
||||
int i, j, t;
|
||||
double norm;
|
||||
|
||||
for (i = 0; i < nx*m; i++)
|
||||
h[i] = 0;
|
||||
|
||||
for (i = hlen/2; i < nx-hlen/2; i++)
|
||||
{
|
||||
for (j = 0; j < m; j++)
|
||||
h[i*m+j] = 0;
|
||||
for (t = i-hlen/2; t <= i+hlen/2; t++)
|
||||
++h[i*m+x[t]];
|
||||
norm = 0;
|
||||
for (j = 0; j < m; j++)
|
||||
norm += h[i*m+j] * h[i*m+j];
|
||||
for (j = 0; j < m; j++)
|
||||
h[i*m+j] /= norm;
|
||||
}
|
||||
|
||||
/* duplicate histograms at beginning and end to create one histogram for each data value supplied */
|
||||
for (i = 0; i < hlen/2; i++)
|
||||
for (j = 0; j < m; j++)
|
||||
h[i*m+j] = h[hlen/2*m+j];
|
||||
for (i = nx-hlen/2; i < nx; i++)
|
||||
for (j = 0; j < m; j++)
|
||||
h[i*m+j] = h[(nx-hlen/2-1)*m+j];
|
||||
}
|
||||
|
||||
/* segment using HMM and then histogram clustering */
|
||||
void cluster_segment(int* q, double** features, int frames_read, int feature_length, int nHMM_states,
|
||||
int histogram_length, int nclusters, int neighbour_limit)
|
||||
{
|
||||
int i, j;
|
||||
|
||||
/*****************************/
|
||||
if (0) {
|
||||
/* try just using the predominant bin number as a 'decoded state' */
|
||||
nHMM_states = feature_length + 1; /* allow a 'zero' state */
|
||||
double chroma_thresh = 0.05;
|
||||
double maxval;
|
||||
int maxbin;
|
||||
for (i = 0; i < frames_read; i++)
|
||||
{
|
||||
maxval = 0;
|
||||
for (j = 0; j < feature_length; j++)
|
||||
{
|
||||
if (features[i][j] > maxval)
|
||||
{
|
||||
maxval = features[i][j];
|
||||
maxbin = j;
|
||||
}
|
||||
}
|
||||
if (maxval > chroma_thresh)
|
||||
q[i] = maxbin;
|
||||
else
|
||||
q[i] = feature_length;
|
||||
}
|
||||
|
||||
}
|
||||
if (1) {
|
||||
/*****************************/
|
||||
|
||||
|
||||
/* scale all the features to 'balance covariances' during HMM training */
|
||||
double scale = 10;
|
||||
for (i = 0; i < frames_read; i++)
|
||||
for (j = 0; j < feature_length; j++)
|
||||
features[i][j] *= scale;
|
||||
|
||||
/* train an HMM on the features */
|
||||
|
||||
/* create a model */
|
||||
model_t* model = hmm_init(features, frames_read, feature_length, nHMM_states);
|
||||
|
||||
/* train the model */
|
||||
hmm_train(features, frames_read, model);
|
||||
/*
|
||||
printf("\n\nafter training:\n");
|
||||
hmm_print(model);
|
||||
*/
|
||||
/* decode the hidden state sequence */
|
||||
viterbi_decode(features, frames_read, model, q);
|
||||
hmm_close(model);
|
||||
|
||||
/*****************************/
|
||||
}
|
||||
/*****************************/
|
||||
|
||||
|
||||
/*
|
||||
fprintf(stderr, "HMM state sequence:\n");
|
||||
for (i = 0; i < frames_read; i++)
|
||||
fprintf(stderr, "%d ", q[i]);
|
||||
fprintf(stderr, "\n\n");
|
||||
*/
|
||||
|
||||
/* create histograms of states */
|
||||
double* h = (double*) malloc(frames_read*nHMM_states*sizeof(double)); /* vector in row major order */
|
||||
create_histograms(q, frames_read, nHMM_states, histogram_length, h);
|
||||
|
||||
/* cluster the histograms */
|
||||
int nbsched = 20; /* length of inverse temperature schedule */
|
||||
double* bsched = (double*) malloc(nbsched*sizeof(double)); /* inverse temperature schedule */
|
||||
double b0 = 100;
|
||||
double alpha = 0.7;
|
||||
bsched[0] = b0;
|
||||
for (i = 1; i < nbsched; i++)
|
||||
bsched[i] = alpha * bsched[i-1];
|
||||
cluster_melt(h, nHMM_states, frames_read, bsched, nbsched, nclusters, neighbour_limit, q);
|
||||
|
||||
/* now q holds a sequence of cluster assignments */
|
||||
|
||||
free(h);
|
||||
free(bsched);
|
||||
}
|
||||
|
||||
/* segment constant-Q or chroma features */
|
||||
void constq_segment(int* q, double** features, int frames_read, int bins, int ncoeff, int feature_type,
|
||||
int nHMM_states, int histogram_length, int nclusters, int neighbour_limit)
|
||||
{
|
||||
int feature_length;
|
||||
double** chroma;
|
||||
int i;
|
||||
|
||||
if (feature_type == FEATURE_TYPE_CONSTQ)
|
||||
{
|
||||
/* fprintf(stderr, "Converting to dB and normalising...\n");
|
||||
*/
|
||||
mpeg7_constq(features, frames_read, ncoeff);
|
||||
/*
|
||||
fprintf(stderr, "Running PCA...\n");
|
||||
*/
|
||||
/* do PCA on the features (but not the envelope) */
|
||||
int ncomponents = 20;
|
||||
pca_project(features, frames_read, ncoeff, ncomponents);
|
||||
|
||||
/* copy the envelope so that it immediatly follows the chosen components */
|
||||
for (i = 0; i < frames_read; i++)
|
||||
features[i][ncomponents] = features[i][ncoeff];
|
||||
|
||||
feature_length = ncomponents + 1;
|
||||
|
||||
/**************************************
|
||||
//TEST
|
||||
// feature file name
|
||||
char* dir = "/Users/mark/documents/semma/audio/";
|
||||
char* file_name = (char*) malloc((strlen(dir) + strlen(trackname) + strlen("_features_c20r8h0.2f0.6.mat") + 1)*sizeof(char));
|
||||
strcpy(file_name, dir);
|
||||
strcat(file_name, trackname);
|
||||
strcat(file_name, "_features_c20r8h0.2f0.6.mat");
|
||||
|
||||
// get the features from Matlab from mat-file
|
||||
int frames_in_file;
|
||||
readmatarray_size(file_name, 2, &frames_in_file, &feature_length);
|
||||
readmatarray(file_name, 2, frames_in_file, feature_length, features);
|
||||
// copy final frame to ensure that we get as many as we expected
|
||||
int missing_frames = frames_read - frames_in_file;
|
||||
while (missing_frames > 0)
|
||||
{
|
||||
for (i = 0; i < feature_length; i++)
|
||||
features[frames_read-missing_frames][i] = features[frames_read-missing_frames-1][i];
|
||||
--missing_frames;
|
||||
}
|
||||
|
||||
free(file_name);
|
||||
******************************************/
|
||||
|
||||
cluster_segment(q, features, frames_read, feature_length, nHMM_states, histogram_length, nclusters, neighbour_limit);
|
||||
}
|
||||
|
||||
if (feature_type == FEATURE_TYPE_CHROMA)
|
||||
{
|
||||
/*
|
||||
fprintf(stderr, "Converting to chroma features...\n");
|
||||
*/
|
||||
/* convert constant-Q to normalised chroma features */
|
||||
chroma = (double**) malloc(frames_read*sizeof(double*));
|
||||
for (i = 0; i < frames_read; i++)
|
||||
chroma[i] = (double*) malloc(bins*sizeof(double));
|
||||
cq2chroma(features, frames_read, ncoeff, bins, chroma);
|
||||
feature_length = bins;
|
||||
|
||||
cluster_segment(q, chroma, frames_read, feature_length, nHMM_states, histogram_length, nclusters, neighbour_limit);
|
||||
|
||||
for (i = 0; i < frames_read; i++)
|
||||
free(chroma[i]);
|
||||
free(chroma);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
|
@ -1,51 +0,0 @@
|
|||
#ifndef _CLUSTER_SEGMENTER_H
|
||||
#define _CLUSTER_SEGMENTER_H
|
||||
|
||||
/*
|
||||
* cluster_segmenter.h
|
||||
* soundbite
|
||||
*
|
||||
* Created by Mark Levy on 06/04/2006.
|
||||
* Copyright 2006 Centre for Digital Music, Queen Mary, University of London.
|
||||
|
||||
This program is free software; you can redistribute it and/or
|
||||
modify it under the terms of the GNU General Public License as
|
||||
published by the Free Software Foundation; either version 2 of the
|
||||
License, or (at your option) any later version. See the file
|
||||
COPYING included with this distribution for more information.
|
||||
*
|
||||
*/
|
||||
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <math.h>
|
||||
#include <float.h>
|
||||
|
||||
#include "segment.h"
|
||||
#include "cluster_melt.h"
|
||||
#include "hmm/hmm.h"
|
||||
#include "maths/pca/pca.h"
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
|
||||
/* applies MPEG-7 normalisation to constant-Q features, storing normalised envelope (norm) in last feature dimension */
|
||||
void mpeg7_constq(double** features, int nframes, int ncoeff);
|
||||
|
||||
/* converts constant-Q features to normalised chroma */
|
||||
void cq2chroma(double** cq, int nframes, int ncoeff, int bins, double** chroma);
|
||||
|
||||
void create_histograms(int* x, int nx, int m, int hlen, double* h);
|
||||
|
||||
void cluster_segment(int* q, double** features, int frames_read, int feature_length, int nHMM_states,
|
||||
int histogram_length, int nclusters, int neighbour_limit);
|
||||
|
||||
void constq_segment(int* q, double** features, int frames_read, int bins, int ncoeff, int feature_type,
|
||||
int nHMM_states, int histogram_length, int nclusters, int neighbour_limit);
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
||||
#endif
|
||||
|
|
@ -1,50 +0,0 @@
|
|||
#ifndef _SEGMENT_H
|
||||
#define _SEGMENT_H
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
|
||||
/*
|
||||
* segment.h
|
||||
*
|
||||
* Created by Mark Levy on 06/04/2006.
|
||||
* Copyright 2006 Centre for Digital Music, Queen Mary, University of London.
|
||||
|
||||
This program is free software; you can redistribute it and/or
|
||||
modify it under the terms of the GNU General Public License as
|
||||
published by the Free Software Foundation; either version 2 of the
|
||||
License, or (at your option) any later version. See the file
|
||||
COPYING included with this distribution for more information.
|
||||
*
|
||||
*/
|
||||
|
||||
typedef struct segment_t
|
||||
{
|
||||
long start; /* in samples */
|
||||
long end;
|
||||
int type;
|
||||
} segment_t;
|
||||
|
||||
typedef struct segmentation_t
|
||||
{
|
||||
int nsegs; /* number of segments */
|
||||
int nsegtypes; /* number of segment types, so possible types are {0,1,...,nsegtypes-1} */
|
||||
int samplerate;
|
||||
segment_t* segments;
|
||||
} segmentation_t;
|
||||
|
||||
typedef enum
|
||||
{
|
||||
FEATURE_TYPE_UNKNOWN = 0,
|
||||
FEATURE_TYPE_CONSTQ = 1,
|
||||
FEATURE_TYPE_CHROMA = 2,
|
||||
FEATURE_TYPE_MFCC = 3
|
||||
} feature_types;
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
||||
#endif
|
||||
|
||||
|
|
@ -1,837 +0,0 @@
|
|||
/*
|
||||
* hmm.c
|
||||
*
|
||||
* Created by Mark Levy on 12/02/2006.
|
||||
* Copyright 2006 Centre for Digital Music, Queen Mary, University of London.
|
||||
|
||||
This program is free software; you can redistribute it and/or
|
||||
modify it under the terms of the GNU General Public License as
|
||||
published by the Free Software Foundation; either version 2 of the
|
||||
License, or (at your option) any later version. See the file
|
||||
COPYING included with this distribution for more information.
|
||||
*
|
||||
*/
|
||||
|
||||
#include <stdio.h>
|
||||
#include <math.h>
|
||||
#include <stdlib.h>
|
||||
#include <float.h>
|
||||
#include <time.h> /* to seed random number generator */
|
||||
|
||||
#include <clapack.h> /* LAPACK for matrix inversion */
|
||||
|
||||
#include "maths/nan-inf.h"
|
||||
|
||||
#ifdef ATLAS_ORDER
|
||||
#define HAVE_ATLAS 1
|
||||
#endif
|
||||
|
||||
#ifdef HAVE_ATLAS
|
||||
// Using ATLAS C interface to LAPACK
|
||||
#define dgetrf_(m, n, a, lda, ipiv, info) \
|
||||
clapack_dgetrf(CblasColMajor, *m, *n, a, *lda, ipiv)
|
||||
#define dgetri_(n, a, lda, ipiv, work, lwork, info) \
|
||||
clapack_dgetri(CblasColMajor, *n, a, *lda, ipiv)
|
||||
#endif
|
||||
|
||||
#ifdef _MAC_OS_X
|
||||
#include <vecLib/cblas.h>
|
||||
#else
|
||||
#include <cblas.h> /* BLAS for matrix multiplication */
|
||||
#endif
|
||||
|
||||
#include "hmm.h"
|
||||
|
||||
model_t* hmm_init(double** x, int T, int L, int N)
|
||||
{
|
||||
int i, j, d, e, t;
|
||||
double s, ss;
|
||||
|
||||
model_t* model;
|
||||
model = (model_t*) malloc(sizeof(model_t));
|
||||
model->N = N;
|
||||
model->L = L;
|
||||
model->p0 = (double*) malloc(N*sizeof(double));
|
||||
model->a = (double**) malloc(N*sizeof(double*));
|
||||
model->mu = (double**) malloc(N*sizeof(double*));
|
||||
for (i = 0; i < N; i++)
|
||||
{
|
||||
model->a[i] = (double*) malloc(N*sizeof(double));
|
||||
model->mu[i] = (double*) malloc(L*sizeof(double));
|
||||
}
|
||||
model->cov = (double**) malloc(L*sizeof(double*));
|
||||
for (i = 0; i < L; i++)
|
||||
model->cov[i] = (double*) malloc(L*sizeof(double));
|
||||
|
||||
srand(time(0));
|
||||
double* global_mean = (double*) malloc(L*sizeof(double));
|
||||
|
||||
/* find global mean */
|
||||
for (d = 0; d < L; d++)
|
||||
{
|
||||
global_mean[d] = 0;
|
||||
for (t = 0; t < T; t++)
|
||||
global_mean[d] += x[t][d];
|
||||
global_mean[d] /= T;
|
||||
}
|
||||
|
||||
/* calculate global diagonal covariance */
|
||||
for (d = 0; d < L; d++)
|
||||
{
|
||||
for (e = 0; e < L; e++)
|
||||
model->cov[d][e] = 0;
|
||||
for (t = 0; t < T; t++)
|
||||
model->cov[d][d] += (x[t][d] - global_mean[d]) * (x[t][d] - global_mean[d]);
|
||||
model->cov[d][d] /= T-1;
|
||||
}
|
||||
|
||||
/* set all means close to global mean */
|
||||
for (i = 0; i < N; i++)
|
||||
{
|
||||
for (d = 0; d < L; d++)
|
||||
{
|
||||
/* add some random noise related to covariance */
|
||||
/* ideally the random number would be Gaussian(0,1), as a hack we make it uniform on [-0.25,0.25] */
|
||||
model->mu[i][d] = global_mean[d] + (0.5 * rand() / (double) RAND_MAX - 0.25) * sqrt(model->cov[d][d]);
|
||||
}
|
||||
}
|
||||
|
||||
/* random intial and transition probs */
|
||||
s = 0;
|
||||
for (i = 0; i < N; i++)
|
||||
{
|
||||
model->p0[i] = 1 + rand() / (double) RAND_MAX;
|
||||
s += model->p0[i];
|
||||
ss = 0;
|
||||
for (j = 0; j < N; j++)
|
||||
{
|
||||
model->a[i][j] = 1 + rand() / (double) RAND_MAX;
|
||||
ss += model->a[i][j];
|
||||
}
|
||||
for (j = 0; j < N; j++)
|
||||
{
|
||||
model->a[i][j] /= ss;
|
||||
}
|
||||
}
|
||||
for (i = 0; i < N; i++)
|
||||
model->p0[i] /= s;
|
||||
|
||||
free(global_mean);
|
||||
|
||||
return model;
|
||||
}
|
||||
|
||||
void hmm_close(model_t* model)
|
||||
{
|
||||
int i;
|
||||
|
||||
for (i = 0; i < model->N; i++)
|
||||
{
|
||||
free(model->a[i]);
|
||||
free(model->mu[i]);
|
||||
}
|
||||
free(model->a);
|
||||
free(model->mu);
|
||||
for (i = 0; i < model->L; i++)
|
||||
free(model->cov[i]);
|
||||
free(model->cov);
|
||||
free(model);
|
||||
}
|
||||
|
||||
void hmm_train(double** x, int T, model_t* model)
|
||||
{
|
||||
int i, t;
|
||||
double loglik; /* overall log-likelihood at each iteration */
|
||||
|
||||
int N = model->N;
|
||||
int L = model->L;
|
||||
double* p0 = model->p0;
|
||||
double** a = model->a;
|
||||
double** mu = model->mu;
|
||||
double** cov = model->cov;
|
||||
|
||||
/* allocate memory */
|
||||
double** gamma = (double**) malloc(T*sizeof(double*));
|
||||
double*** xi = (double***) malloc(T*sizeof(double**));
|
||||
for (t = 0; t < T; t++)
|
||||
{
|
||||
gamma[t] = (double*) malloc(N*sizeof(double));
|
||||
xi[t] = (double**) malloc(N*sizeof(double*));
|
||||
for (i = 0; i < N; i++)
|
||||
xi[t][i] = (double*) malloc(N*sizeof(double));
|
||||
}
|
||||
|
||||
/* temporary memory */
|
||||
double* gauss_y = (double*) malloc(L*sizeof(double));
|
||||
double* gauss_z = (double*) malloc(L*sizeof(double));
|
||||
|
||||
/* obs probs P(j|{x}) */
|
||||
double** b = (double**) malloc(T*sizeof(double*));
|
||||
for (t = 0; t < T; t++)
|
||||
b[t] = (double*) malloc(N*sizeof(double));
|
||||
|
||||
/* inverse covariance and its determinant */
|
||||
double** icov = (double**) malloc(L*sizeof(double*));
|
||||
for (i = 0; i < L; i++)
|
||||
icov[i] = (double*) malloc(L*sizeof(double));
|
||||
double detcov;
|
||||
|
||||
double thresh = 0.0001;
|
||||
int niter = 50;
|
||||
int iter = 0;
|
||||
double loglik1, loglik2;
|
||||
int foundnan = 0;
|
||||
|
||||
while (iter < niter && !foundnan && !(iter > 1 && (loglik - loglik1) < thresh * (loglik1 - loglik2)))
|
||||
{
|
||||
++iter;
|
||||
/*
|
||||
fprintf(stderr, "calculating obsprobs...\n");
|
||||
fflush(stderr);
|
||||
*/
|
||||
/* precalculate obs probs */
|
||||
invert(cov, L, icov, &detcov);
|
||||
|
||||
for (t = 0; t < T; t++)
|
||||
{
|
||||
//int allzero = 1;
|
||||
for (i = 0; i < N; i++)
|
||||
{
|
||||
b[t][i] = exp(loggauss(x[t], L, mu[i], icov, detcov, gauss_y, gauss_z));
|
||||
|
||||
//if (b[t][i] != 0)
|
||||
// allzero = 0;
|
||||
}
|
||||
/*
|
||||
if (allzero)
|
||||
{
|
||||
printf("all the b[t][i] were zero for t = %d, correcting...\n", t);
|
||||
for (i = 0; i < N; i++)
|
||||
{
|
||||
b[t][i] = 0.00001;
|
||||
}
|
||||
}
|
||||
*/
|
||||
}
|
||||
/*
|
||||
fprintf(stderr, "forwards-backwards...\n");
|
||||
fflush(stderr);
|
||||
*/
|
||||
forward_backwards(xi, gamma, &loglik, &loglik1, &loglik2, iter, N, T, p0, a, b);
|
||||
/*
|
||||
fprintf(stderr, "iteration %d: loglik = %f\n", iter, loglik);
|
||||
fprintf(stderr, "re-estimation...\n");
|
||||
fflush(stderr);
|
||||
*/
|
||||
if (ISNAN(loglik)) {
|
||||
foundnan = 1;
|
||||
continue;
|
||||
}
|
||||
|
||||
baum_welch(p0, a, mu, cov, N, T, L, x, xi, gamma);
|
||||
|
||||
/*
|
||||
printf("a:\n");
|
||||
for (i = 0; i < model->N; i++)
|
||||
{
|
||||
for (j = 0; j < model->N; j++)
|
||||
printf("%f ", model->a[i][j]);
|
||||
printf("\n");
|
||||
}
|
||||
printf("\n\n");
|
||||
*/
|
||||
//hmm_print(model);
|
||||
}
|
||||
|
||||
/* deallocate memory */
|
||||
for (t = 0; t < T; t++)
|
||||
{
|
||||
free(gamma[t]);
|
||||
free(b[t]);
|
||||
for (i = 0; i < N; i++)
|
||||
free(xi[t][i]);
|
||||
free(xi[t]);
|
||||
}
|
||||
free(gamma);
|
||||
free(xi);
|
||||
free(b);
|
||||
|
||||
for (i = 0; i < L; i++)
|
||||
free(icov[i]);
|
||||
free(icov);
|
||||
|
||||
free(gauss_y);
|
||||
free(gauss_z);
|
||||
}
|
||||
|
||||
void mlss_reestimate(double* p0, double** a, double** mu, double** cov, int N, int T, int L, int* q, double** x)
|
||||
{
|
||||
/* fit a single Gaussian to observations in each state */
|
||||
|
||||
/* calculate the mean observation in each state */
|
||||
|
||||
/* calculate the overall covariance */
|
||||
|
||||
/* count transitions */
|
||||
|
||||
/* estimate initial probs from transitions (???) */
|
||||
}
|
||||
|
||||
void baum_welch(double* p0, double** a, double** mu, double** cov, int N, int T, int L, double** x, double*** xi, double** gamma)
|
||||
{
|
||||
int i, j, t;
|
||||
|
||||
double* sum_gamma = (double*) malloc(N*sizeof(double));
|
||||
|
||||
/* temporary memory */
|
||||
double* u = (double*) malloc(L*L*sizeof(double));
|
||||
double* yy = (double*) malloc(T*L*sizeof(double));
|
||||
double* yy2 = (double*) malloc(T*L*sizeof(double));
|
||||
|
||||
/* re-estimate transition probs */
|
||||
for (i = 0; i < N; i++)
|
||||
{
|
||||
sum_gamma[i] = 0;
|
||||
for (t = 0; t < T-1; t++)
|
||||
sum_gamma[i] += gamma[t][i];
|
||||
}
|
||||
|
||||
for (i = 0; i < N; i++)
|
||||
{
|
||||
if (sum_gamma[i] == 0)
|
||||
{
|
||||
/* fprintf(stderr, "sum_gamma[%d] was zero...\n", i); */
|
||||
}
|
||||
//double s = 0;
|
||||
for (j = 0; j < N; j++)
|
||||
{
|
||||
a[i][j] = 0;
|
||||
if (sum_gamma[i] == 0.) continue;
|
||||
for (t = 0; t < T-1; t++)
|
||||
a[i][j] += xi[t][i][j];
|
||||
//s += a[i][j];
|
||||
a[i][j] /= sum_gamma[i];
|
||||
}
|
||||
/*
|
||||
for (j = 0; j < N; j++)
|
||||
{
|
||||
a[i][j] /= s;
|
||||
}
|
||||
*/
|
||||
}
|
||||
|
||||
/* NB: now we need to sum gamma over all t */
|
||||
for (i = 0; i < N; i++)
|
||||
sum_gamma[i] += gamma[T-1][i];
|
||||
|
||||
/* re-estimate initial probs */
|
||||
for (i = 0; i < N; i++)
|
||||
p0[i] = gamma[0][i];
|
||||
|
||||
/* re-estimate covariance */
|
||||
int d, e;
|
||||
double sum_sum_gamma = 0;
|
||||
for (i = 0; i < N; i++)
|
||||
sum_sum_gamma += sum_gamma[i];
|
||||
|
||||
/*
|
||||
for (d = 0; d < L; d++)
|
||||
{
|
||||
for (e = d; e < L; e++)
|
||||
{
|
||||
cov[d][e] = 0;
|
||||
for (t = 0; t < T; t++)
|
||||
for (j = 0; j < N; j++)
|
||||
cov[d][e] += gamma[t][j] * (x[t][d] - mu[j][d]) * (x[t][e] - mu[j][e]);
|
||||
|
||||
cov[d][e] /= sum_sum_gamma;
|
||||
|
||||
if (ISNAN(cov[d][e]))
|
||||
{
|
||||
printf("cov[%d][%d] was nan\n", d, e);
|
||||
for (j = 0; j < N; j++)
|
||||
for (i = 0; i < L; i++)
|
||||
if (ISNAN(mu[j][i]))
|
||||
printf("mu[%d][%d] was nan\n", j, i);
|
||||
for (t = 0; t < T; t++)
|
||||
for (j = 0; j < N; j++)
|
||||
if (ISNAN(gamma[t][j]))
|
||||
printf("gamma[%d][%d] was nan\n", t, j);
|
||||
exit(-1);
|
||||
}
|
||||
}
|
||||
}
|
||||
for (d = 0; d < L; d++)
|
||||
for (e = 0; e < d; e++)
|
||||
cov[d][e] = cov[e][d];
|
||||
*/
|
||||
|
||||
/* using BLAS */
|
||||
for (d = 0; d < L; d++)
|
||||
for (e = 0; e < L; e++)
|
||||
cov[d][e] = 0;
|
||||
|
||||
for (j = 0; j < N; j++)
|
||||
{
|
||||
for (d = 0; d < L; d++)
|
||||
for (t = 0; t < T; t++)
|
||||
{
|
||||
yy[d*T+t] = x[t][d] - mu[j][d];
|
||||
yy2[d*T+t] = gamma[t][j] * (x[t][d] - mu[j][d]);
|
||||
}
|
||||
|
||||
cblas_dgemm(CblasColMajor, CblasTrans, CblasNoTrans, L, L, T, 1.0, yy, T, yy2, T, 0, u, L);
|
||||
|
||||
for (e = 0; e < L; e++)
|
||||
for (d = 0; d < L; d++)
|
||||
cov[d][e] += u[e*L+d];
|
||||
}
|
||||
|
||||
for (d = 0; d < L; d++)
|
||||
for (e = 0; e < L; e++)
|
||||
cov[d][e] /= T; /* sum_sum_gamma; */
|
||||
|
||||
//printf("sum_sum_gamma = %f\n", sum_sum_gamma); /* fine, = T IS THIS ALWAYS TRUE with pooled cov?? */
|
||||
|
||||
/* re-estimate means */
|
||||
for (j = 0; j < N; j++)
|
||||
{
|
||||
for (d = 0; d < L; d++)
|
||||
{
|
||||
mu[j][d] = 0;
|
||||
for (t = 0; t < T; t++)
|
||||
mu[j][d] += gamma[t][j] * x[t][d];
|
||||
mu[j][d] /= sum_gamma[j];
|
||||
}
|
||||
}
|
||||
|
||||
/* deallocate memory */
|
||||
free(sum_gamma);
|
||||
free(yy);
|
||||
free(yy2);
|
||||
free(u);
|
||||
}
|
||||
|
||||
void forward_backwards(double*** xi, double** gamma, double* loglik, double* loglik1, double* loglik2, int iter, int N, int T, double* p0, double** a, double** b)
|
||||
{
|
||||
/* forwards-backwards with scaling */
|
||||
int i, j, t;
|
||||
|
||||
double** alpha = (double**) malloc(T*sizeof(double*));
|
||||
double** beta = (double**) malloc(T*sizeof(double*));
|
||||
for (t = 0; t < T; t++)
|
||||
{
|
||||
alpha[t] = (double*) malloc(N*sizeof(double));
|
||||
beta[t] = (double*) malloc(N*sizeof(double));
|
||||
}
|
||||
|
||||
/* scaling coefficients */
|
||||
double* c = (double*) malloc(T*sizeof(double));
|
||||
|
||||
/* calculate forward probs and scale coefficients */
|
||||
c[0] = 0;
|
||||
for (i = 0; i < N; i++)
|
||||
{
|
||||
alpha[0][i] = p0[i] * b[0][i];
|
||||
c[0] += alpha[0][i];
|
||||
|
||||
//printf("p0[%d] = %f, b[0][%d] = %f\n", i, p0[i], i, b[0][i]);
|
||||
}
|
||||
c[0] = 1 / c[0];
|
||||
for (i = 0; i < N; i++)
|
||||
{
|
||||
alpha[0][i] *= c[0];
|
||||
|
||||
//printf("alpha[0][%d] = %f\n", i, alpha[0][i]); /* OK agrees with Matlab */
|
||||
}
|
||||
|
||||
*loglik1 = *loglik;
|
||||
*loglik = -log(c[0]);
|
||||
if (iter == 2)
|
||||
*loglik2 = *loglik;
|
||||
|
||||
for (t = 1; t < T; t++)
|
||||
{
|
||||
c[t] = 0;
|
||||
for (j = 0; j < N; j++)
|
||||
{
|
||||
alpha[t][j] = 0;
|
||||
for (i = 0; i < N; i++)
|
||||
alpha[t][j] += alpha[t-1][i] * a[i][j];
|
||||
alpha[t][j] *= b[t][j];
|
||||
|
||||
c[t] += alpha[t][j];
|
||||
}
|
||||
|
||||
/*
|
||||
if (c[t] == 0)
|
||||
{
|
||||
printf("c[%d] = 0, going to blow up so exiting\n", t);
|
||||
for (i = 0; i < N; i++)
|
||||
if (b[t][i] == 0)
|
||||
fprintf(stderr, "b[%d][%d] was zero\n", t, i);
|
||||
fprintf(stderr, "x[t] was \n");
|
||||
for (i = 0; i < L; i++)
|
||||
fprintf(stderr, "%f ", x[t][i]);
|
||||
fprintf(stderr, "\n\n");
|
||||
exit(-1);
|
||||
}
|
||||
*/
|
||||
|
||||
c[t] = 1 / c[t];
|
||||
for (j = 0; j < N; j++)
|
||||
alpha[t][j] *= c[t];
|
||||
|
||||
//printf("c[%d] = %e\n", t, c[t]);
|
||||
|
||||
*loglik -= log(c[t]);
|
||||
}
|
||||
|
||||
/* calculate backwards probs using same coefficients */
|
||||
for (i = 0; i < N; i++)
|
||||
beta[T-1][i] = 1;
|
||||
t = T - 1;
|
||||
while (1)
|
||||
{
|
||||
for (i = 0; i < N; i++)
|
||||
beta[t][i] *= c[t];
|
||||
|
||||
if (t == 0)
|
||||
break;
|
||||
|
||||
for (i = 0; i < N; i++)
|
||||
{
|
||||
beta[t-1][i] = 0;
|
||||
for (j = 0; j < N; j++)
|
||||
beta[t-1][i] += a[i][j] * b[t][j] * beta[t][j];
|
||||
}
|
||||
|
||||
t--;
|
||||
}
|
||||
|
||||
/*
|
||||
printf("alpha:\n");
|
||||
for (t = 0; t < T; t++)
|
||||
{
|
||||
for (i = 0; i < N; i++)
|
||||
printf("%4.4e\t\t", alpha[t][i]);
|
||||
printf("\n");
|
||||
}
|
||||
printf("\n\n");printf("beta:\n");
|
||||
for (t = 0; t < T; t++)
|
||||
{
|
||||
for (i = 0; i < N; i++)
|
||||
printf("%4.4e\t\t", beta[t][i]);
|
||||
printf("\n");
|
||||
}
|
||||
printf("\n\n");
|
||||
*/
|
||||
|
||||
/* calculate posterior probs */
|
||||
double tot;
|
||||
for (t = 0; t < T; t++)
|
||||
{
|
||||
tot = 0;
|
||||
for (i = 0; i < N; i++)
|
||||
{
|
||||
gamma[t][i] = alpha[t][i] * beta[t][i];
|
||||
tot += gamma[t][i];
|
||||
}
|
||||
for (i = 0; i < N; i++)
|
||||
{
|
||||
gamma[t][i] /= tot;
|
||||
|
||||
//printf("gamma[%d][%d] = %f\n", t, i, gamma[t][i]);
|
||||
}
|
||||
}
|
||||
|
||||
for (t = 0; t < T-1; t++)
|
||||
{
|
||||
tot = 0;
|
||||
for (i = 0; i < N; i++)
|
||||
{
|
||||
for (j = 0; j < N; j++)
|
||||
{
|
||||
xi[t][i][j] = alpha[t][i] * a[i][j] * b[t+1][j] * beta[t+1][j];
|
||||
tot += xi[t][i][j];
|
||||
}
|
||||
}
|
||||
for (i = 0; i < N; i++)
|
||||
for (j = 0; j < N; j++)
|
||||
xi[t][i][j] /= tot;
|
||||
}
|
||||
|
||||
/*
|
||||
// CHECK - fine
|
||||
// gamma[t][i] = \sum_j{xi[t][i][j]}
|
||||
tot = 0;
|
||||
for (j = 0; j < N; j++)
|
||||
tot += xi[3][1][j];
|
||||
printf("gamma[3][1] = %f, sum_j(xi[3][1][j]) = %f\n", gamma[3][1], tot);
|
||||
*/
|
||||
|
||||
for (t = 0; t < T; t++)
|
||||
{
|
||||
free(alpha[t]);
|
||||
free(beta[t]);
|
||||
}
|
||||
free(alpha);
|
||||
free(beta);
|
||||
free(c);
|
||||
}
|
||||
|
||||
void viterbi_decode(double** x, int T, model_t* model, int* q)
|
||||
{
|
||||
int i, j, t;
|
||||
double max;
|
||||
int havemax;
|
||||
|
||||
int N = model->N;
|
||||
int L = model->L;
|
||||
double* p0 = model->p0;
|
||||
double** a = model->a;
|
||||
double** mu = model->mu;
|
||||
double** cov = model->cov;
|
||||
|
||||
/* inverse covariance and its determinant */
|
||||
double** icov = (double**) malloc(L*sizeof(double*));
|
||||
for (i = 0; i < L; i++)
|
||||
icov[i] = (double*) malloc(L*sizeof(double));
|
||||
double detcov;
|
||||
|
||||
double** logb = (double**) malloc(T*sizeof(double*));
|
||||
double** phi = (double**) malloc(T*sizeof(double*));
|
||||
int** psi = (int**) malloc(T*sizeof(int*));
|
||||
for (t = 0; t < T; t++)
|
||||
{
|
||||
logb[t] = (double*) malloc(N*sizeof(double));
|
||||
phi[t] = (double*) malloc(N*sizeof(double));
|
||||
psi[t] = (int*) malloc(N*sizeof(int));
|
||||
}
|
||||
|
||||
/* temporary memory */
|
||||
double* gauss_y = (double*) malloc(L*sizeof(double));
|
||||
double* gauss_z = (double*) malloc(L*sizeof(double));
|
||||
|
||||
/* calculate observation logprobs */
|
||||
invert(cov, L, icov, &detcov);
|
||||
for (t = 0; t < T; t++)
|
||||
for (i = 0; i < N; i++)
|
||||
logb[t][i] = loggauss(x[t], L, mu[i], icov, detcov, gauss_y, gauss_z);
|
||||
|
||||
/* initialise */
|
||||
for (i = 0; i < N; i++)
|
||||
{
|
||||
phi[0][i] = log(p0[i]) + logb[0][i];
|
||||
psi[0][i] = 0;
|
||||
}
|
||||
|
||||
for (t = 1; t < T; t++)
|
||||
{
|
||||
for (j = 0; j < N; j++)
|
||||
{
|
||||
max = -1000000;
|
||||
havemax = 0;
|
||||
|
||||
psi[t][j] = 0;
|
||||
for (i = 0; i < N; i++)
|
||||
{
|
||||
if (phi[t-1][i] + log(a[i][j]) > max || !havemax)
|
||||
{
|
||||
max = phi[t-1][i] + log(a[i][j]);
|
||||
phi[t][j] = max + logb[t][j];
|
||||
psi[t][j] = i;
|
||||
havemax = 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* find maximising state at time T-1 */
|
||||
max = phi[T-1][0];
|
||||
q[T-1] = 0;
|
||||
for (i = 1; i < N; i++)
|
||||
{
|
||||
if (phi[T-1][i] > max)
|
||||
{
|
||||
max = phi[T-1][i];
|
||||
q[T-1] = i;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* track back */
|
||||
t = T - 2;
|
||||
while (t >= 0)
|
||||
{
|
||||
q[t] = psi[t+1][q[t+1]];
|
||||
t--;
|
||||
}
|
||||
|
||||
/* de-allocate memory */
|
||||
for (i = 0; i < L; i++)
|
||||
free(icov[i]);
|
||||
free(icov);
|
||||
for (t = 0; t < T; t++)
|
||||
{
|
||||
free(logb[t]);
|
||||
free(phi[t]);
|
||||
free(psi[t]);
|
||||
}
|
||||
free(logb);
|
||||
free(phi);
|
||||
free(psi);
|
||||
|
||||
free(gauss_y);
|
||||
free(gauss_z);
|
||||
}
|
||||
|
||||
/* invert matrix and calculate determinant using LAPACK */
|
||||
void invert(double** cov, int L, double** icov, double* detcov)
|
||||
{
|
||||
/* copy square matrix into a vector in column-major order */
|
||||
double* a = (double*) malloc(L*L*sizeof(double));
|
||||
int i, j;
|
||||
for(j=0; j < L; j++)
|
||||
for (i=0; i < L; i++)
|
||||
a[j*L+i] = cov[i][j];
|
||||
|
||||
int M = (int) L;
|
||||
int* ipiv = (int *) malloc(L*L*sizeof(int));
|
||||
int ret;
|
||||
|
||||
/* LU decomposition */
|
||||
ret = dgetrf_(&M, &M, a, &M, ipiv, &ret); /* ret should be zero, negative if cov is singular */
|
||||
if (ret < 0)
|
||||
{
|
||||
fprintf(stderr, "Covariance matrix was singular, couldn't invert\n");
|
||||
exit(-1);
|
||||
}
|
||||
|
||||
/* find determinant */
|
||||
double det = 1;
|
||||
for(i = 0; i < L; i++)
|
||||
det *= a[i*L+i];
|
||||
// TODO: get this to work!!! If detcov < 0 then cov is bad anyway...
|
||||
/*
|
||||
int sign = 1;
|
||||
for (i = 0; i < L; i++)
|
||||
if (ipiv[i] != i)
|
||||
sign = -sign;
|
||||
det *= sign;
|
||||
*/
|
||||
if (det < 0)
|
||||
det = -det;
|
||||
*detcov = det;
|
||||
|
||||
/* allocate required working storage */
|
||||
#ifndef HAVE_ATLAS
|
||||
int lwork = -1;
|
||||
double lwbest = 0.0;
|
||||
dgetri_(&M, a, &M, ipiv, &lwbest, &lwork, &ret);
|
||||
lwork = (int) lwbest;
|
||||
double* work = (double*) malloc(lwork*sizeof(double));
|
||||
#endif
|
||||
|
||||
/* find inverse */
|
||||
dgetri_(&M, a, &M, ipiv, work, &lwork, &ret);
|
||||
|
||||
for(j=0; j < L; j++)
|
||||
for (i=0; i < L; i++)
|
||||
icov[i][j] = a[j*L+i];
|
||||
|
||||
#ifndef HAVE_ATLAS
|
||||
free(work);
|
||||
#endif
|
||||
free(a);
|
||||
}
|
||||
|
||||
/* probability of multivariate Gaussian given mean, inverse and determinant of covariance */
|
||||
double gauss(double* x, int L, double* mu, double** icov, double detcov, double* y, double* z)
|
||||
{
|
||||
int i, j;
|
||||
double s = 0;
|
||||
for (i = 0; i < L; i++)
|
||||
y[i] = x[i] - mu[i];
|
||||
for (i = 0; i < L; i++)
|
||||
{
|
||||
//z[i] = 0;
|
||||
//for (j = 0; j < L; j++)
|
||||
// z[i] += icov[i][j] * y[j];
|
||||
z[i] = cblas_ddot(L, &icov[i][0], 1, y, 1);
|
||||
}
|
||||
s = cblas_ddot(L, z, 1, y, 1);
|
||||
//for (i = 0; i < L; i++)
|
||||
// s += z[i] * y[i];
|
||||
|
||||
return exp(-s/2.0) / (pow(2*PI, L/2.0) * sqrt(detcov));
|
||||
}
|
||||
|
||||
/* log probability of multivariate Gaussian given mean, inverse and determinant of covariance */
|
||||
double loggauss(double* x, int L, double* mu, double** icov, double detcov, double* y, double* z)
|
||||
{
|
||||
int i, j;
|
||||
double s = 0;
|
||||
double ret;
|
||||
for (i = 0; i < L; i++)
|
||||
y[i] = x[i] - mu[i];
|
||||
for (i = 0; i < L; i++)
|
||||
{
|
||||
//z[i] = 0;
|
||||
//for (j = 0; j < L; j++)
|
||||
// z[i] += icov[i][j] * y[j];
|
||||
z[i] = cblas_ddot(L, &icov[i][0], 1, y, 1);
|
||||
}
|
||||
s = cblas_ddot(L, z, 1, y, 1);
|
||||
//for (i = 0; i < L; i++)
|
||||
// s += z[i] * y[i];
|
||||
|
||||
ret = -0.5 * (s + L * log(2*PI) + log(detcov));
|
||||
|
||||
/*
|
||||
// TEST
|
||||
if (ISINF(ret) > 0)
|
||||
printf("loggauss returning infinity\n");
|
||||
if (ISINF(ret) < 0)
|
||||
printf("loggauss returning -infinity\n");
|
||||
if (ISNAN(ret))
|
||||
printf("loggauss returning nan\n");
|
||||
*/
|
||||
|
||||
return ret;
|
||||
}
|
||||
|
||||
void hmm_print(model_t* model)
|
||||
{
|
||||
int i, j;
|
||||
printf("p0:\n");
|
||||
for (i = 0; i < model->N; i++)
|
||||
printf("%f ", model->p0[i]);
|
||||
printf("\n\n");
|
||||
printf("a:\n");
|
||||
for (i = 0; i < model->N; i++)
|
||||
{
|
||||
for (j = 0; j < model->N; j++)
|
||||
printf("%f ", model->a[i][j]);
|
||||
printf("\n");
|
||||
}
|
||||
printf("\n\n");
|
||||
printf("mu:\n");
|
||||
for (i = 0; i < model->N; i++)
|
||||
{
|
||||
for (j = 0; j < model->L; j++)
|
||||
printf("%f ", model->mu[i][j]);
|
||||
printf("\n");
|
||||
}
|
||||
printf("\n\n");
|
||||
printf("cov:\n");
|
||||
for (i = 0; i < model->L; i++)
|
||||
{
|
||||
for (j = 0; j < model->L; j++)
|
||||
printf("%f ", model->cov[i][j]);
|
||||
printf("\n");
|
||||
}
|
||||
printf("\n\n");
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -1,52 +0,0 @@
|
|||
#ifndef _HMM_H
|
||||
#define _HMM_H
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
|
||||
/*
|
||||
* hmm.h
|
||||
*
|
||||
* Created by Mark Levy on 12/02/2006.
|
||||
* Copyright 2006 Centre for Digital Music, Queen Mary, University of London.
|
||||
|
||||
This program is free software; you can redistribute it and/or
|
||||
modify it under the terms of the GNU General Public License as
|
||||
published by the Free Software Foundation; either version 2 of the
|
||||
License, or (at your option) any later version. See the file
|
||||
COPYING included with this distribution for more information.
|
||||
*
|
||||
*/
|
||||
|
||||
#ifndef PI
|
||||
#define PI 3.14159265358979323846264338327950288
|
||||
#endif
|
||||
|
||||
typedef struct _model_t {
|
||||
int N; /* number of states */
|
||||
double* p0; /* initial probs */
|
||||
double** a; /* transition probs */
|
||||
int L; /* dimensionality of data */
|
||||
double** mu; /* state means */
|
||||
double** cov; /* covariance, tied between all states */
|
||||
} model_t;
|
||||
|
||||
void hmm_train(double** x, int T, model_t* model); /* with scaling */
|
||||
void forward_backwards(double*** xi, double** gamma, double* loglik, double* loglik1, double* loglik2, int iter,
|
||||
int N, int T, double* p0, double** a, double** b);
|
||||
void baum_welch(double* p0, double** a, double** mu, double** cov, int N, int T, int L, double** x, double*** xi, double** gamma);
|
||||
void viterbi_decode(double** x, int T, model_t* model, int* q); /* using logs */
|
||||
model_t* hmm_init(double** x, int T, int L, int N);
|
||||
void hmm_close(model_t* model);
|
||||
void invert(double** cov, int L, double** icov, double* detcov); /* uses LAPACK (included with Mac OSX) */
|
||||
double gauss(double* x, int L, double* mu, double** icov, double detcov, double* y, double* z);
|
||||
double loggauss(double* x, int L, double* mu, double** icov, double detcov, double* y, double* z);
|
||||
void hmm_print(model_t* model);
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
||||
#endif
|
||||
|
||||
|
|
@ -1,578 +0,0 @@
|
|||
#ifndef CBLAS_H
|
||||
#define CBLAS_H
|
||||
#include <stddef.h>
|
||||
|
||||
/* Allow the use in C++ code. */
|
||||
#ifdef __cplusplus
|
||||
extern "C"
|
||||
{
|
||||
#endif
|
||||
|
||||
/*
|
||||
* Enumerated and derived types
|
||||
*/
|
||||
#define CBLAS_INDEX size_t /* this may vary between platforms */
|
||||
|
||||
enum CBLAS_ORDER {CblasRowMajor=101, CblasColMajor=102};
|
||||
enum CBLAS_TRANSPOSE {CblasNoTrans=111, CblasTrans=112, CblasConjTrans=113};
|
||||
enum CBLAS_UPLO {CblasUpper=121, CblasLower=122};
|
||||
enum CBLAS_DIAG {CblasNonUnit=131, CblasUnit=132};
|
||||
enum CBLAS_SIDE {CblasLeft=141, CblasRight=142};
|
||||
|
||||
/*
|
||||
* ===========================================================================
|
||||
* Prototypes for level 1 BLAS functions (complex are recast as routines)
|
||||
* ===========================================================================
|
||||
*/
|
||||
float cblas_sdsdot(const int N, const float alpha, const float *X,
|
||||
const int incX, const float *Y, const int incY);
|
||||
double cblas_dsdot(const int N, const float *X, const int incX, const float *Y,
|
||||
const int incY);
|
||||
float cblas_sdot(const int N, const float *X, const int incX,
|
||||
const float *Y, const int incY);
|
||||
double cblas_ddot(const int N, const double *X, const int incX,
|
||||
const double *Y, const int incY);
|
||||
|
||||
/*
|
||||
* Functions having prefixes Z and C only
|
||||
*/
|
||||
void cblas_cdotu_sub(const int N, const void *X, const int incX,
|
||||
const void *Y, const int incY, void *dotu);
|
||||
void cblas_cdotc_sub(const int N, const void *X, const int incX,
|
||||
const void *Y, const int incY, void *dotc);
|
||||
|
||||
void cblas_zdotu_sub(const int N, const void *X, const int incX,
|
||||
const void *Y, const int incY, void *dotu);
|
||||
void cblas_zdotc_sub(const int N, const void *X, const int incX,
|
||||
const void *Y, const int incY, void *dotc);
|
||||
|
||||
|
||||
/*
|
||||
* Functions having prefixes S D SC DZ
|
||||
*/
|
||||
float cblas_snrm2(const int N, const float *X, const int incX);
|
||||
float cblas_sasum(const int N, const float *X, const int incX);
|
||||
|
||||
double cblas_dnrm2(const int N, const double *X, const int incX);
|
||||
double cblas_dasum(const int N, const double *X, const int incX);
|
||||
|
||||
float cblas_scnrm2(const int N, const void *X, const int incX);
|
||||
float cblas_scasum(const int N, const void *X, const int incX);
|
||||
|
||||
double cblas_dznrm2(const int N, const void *X, const int incX);
|
||||
double cblas_dzasum(const int N, const void *X, const int incX);
|
||||
|
||||
|
||||
/*
|
||||
* Functions having standard 4 prefixes (S D C Z)
|
||||
*/
|
||||
CBLAS_INDEX cblas_isamax(const int N, const float *X, const int incX);
|
||||
CBLAS_INDEX cblas_idamax(const int N, const double *X, const int incX);
|
||||
CBLAS_INDEX cblas_icamax(const int N, const void *X, const int incX);
|
||||
CBLAS_INDEX cblas_izamax(const int N, const void *X, const int incX);
|
||||
|
||||
/*
|
||||
* ===========================================================================
|
||||
* Prototypes for level 1 BLAS routines
|
||||
* ===========================================================================
|
||||
*/
|
||||
|
||||
/*
|
||||
* Routines with standard 4 prefixes (s, d, c, z)
|
||||
*/
|
||||
void cblas_sswap(const int N, float *X, const int incX,
|
||||
float *Y, const int incY);
|
||||
void cblas_scopy(const int N, const float *X, const int incX,
|
||||
float *Y, const int incY);
|
||||
void cblas_saxpy(const int N, const float alpha, const float *X,
|
||||
const int incX, float *Y, const int incY);
|
||||
|
||||
void cblas_dswap(const int N, double *X, const int incX,
|
||||
double *Y, const int incY);
|
||||
void cblas_dcopy(const int N, const double *X, const int incX,
|
||||
double *Y, const int incY);
|
||||
void cblas_daxpy(const int N, const double alpha, const double *X,
|
||||
const int incX, double *Y, const int incY);
|
||||
|
||||
void cblas_cswap(const int N, void *X, const int incX,
|
||||
void *Y, const int incY);
|
||||
void cblas_ccopy(const int N, const void *X, const int incX,
|
||||
void *Y, const int incY);
|
||||
void cblas_caxpy(const int N, const void *alpha, const void *X,
|
||||
const int incX, void *Y, const int incY);
|
||||
|
||||
void cblas_zswap(const int N, void *X, const int incX,
|
||||
void *Y, const int incY);
|
||||
void cblas_zcopy(const int N, const void *X, const int incX,
|
||||
void *Y, const int incY);
|
||||
void cblas_zaxpy(const int N, const void *alpha, const void *X,
|
||||
const int incX, void *Y, const int incY);
|
||||
|
||||
|
||||
/*
|
||||
* Routines with S and D prefix only
|
||||
*/
|
||||
void cblas_srotg(float *a, float *b, float *c, float *s);
|
||||
void cblas_srotmg(float *d1, float *d2, float *b1, const float b2, float *P);
|
||||
void cblas_srot(const int N, float *X, const int incX,
|
||||
float *Y, const int incY, const float c, const float s);
|
||||
void cblas_srotm(const int N, float *X, const int incX,
|
||||
float *Y, const int incY, const float *P);
|
||||
|
||||
void cblas_drotg(double *a, double *b, double *c, double *s);
|
||||
void cblas_drotmg(double *d1, double *d2, double *b1, const double b2, double *P);
|
||||
void cblas_drot(const int N, double *X, const int incX,
|
||||
double *Y, const int incY, const double c, const double s);
|
||||
void cblas_drotm(const int N, double *X, const int incX,
|
||||
double *Y, const int incY, const double *P);
|
||||
|
||||
|
||||
/*
|
||||
* Routines with S D C Z CS and ZD prefixes
|
||||
*/
|
||||
void cblas_sscal(const int N, const float alpha, float *X, const int incX);
|
||||
void cblas_dscal(const int N, const double alpha, double *X, const int incX);
|
||||
void cblas_cscal(const int N, const void *alpha, void *X, const int incX);
|
||||
void cblas_zscal(const int N, const void *alpha, void *X, const int incX);
|
||||
void cblas_csscal(const int N, const float alpha, void *X, const int incX);
|
||||
void cblas_zdscal(const int N, const double alpha, void *X, const int incX);
|
||||
|
||||
/*
|
||||
* ===========================================================================
|
||||
* Prototypes for level 2 BLAS
|
||||
* ===========================================================================
|
||||
*/
|
||||
|
||||
/*
|
||||
* Routines with standard 4 prefixes (S, D, C, Z)
|
||||
*/
|
||||
void cblas_sgemv(const enum CBLAS_ORDER order,
|
||||
const enum CBLAS_TRANSPOSE TransA, const int M, const int N,
|
||||
const float alpha, const float *A, const int lda,
|
||||
const float *X, const int incX, const float beta,
|
||||
float *Y, const int incY);
|
||||
void cblas_sgbmv(const enum CBLAS_ORDER order,
|
||||
const enum CBLAS_TRANSPOSE TransA, const int M, const int N,
|
||||
const int KL, const int KU, const float alpha,
|
||||
const float *A, const int lda, const float *X,
|
||||
const int incX, const float beta, float *Y, const int incY);
|
||||
void cblas_strmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
|
||||
const int N, const float *A, const int lda,
|
||||
float *X, const int incX);
|
||||
void cblas_stbmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
|
||||
const int N, const int K, const float *A, const int lda,
|
||||
float *X, const int incX);
|
||||
void cblas_stpmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
|
||||
const int N, const float *Ap, float *X, const int incX);
|
||||
void cblas_strsv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
|
||||
const int N, const float *A, const int lda, float *X,
|
||||
const int incX);
|
||||
void cblas_stbsv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
|
||||
const int N, const int K, const float *A, const int lda,
|
||||
float *X, const int incX);
|
||||
void cblas_stpsv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
|
||||
const int N, const float *Ap, float *X, const int incX);
|
||||
|
||||
void cblas_dgemv(const enum CBLAS_ORDER order,
|
||||
const enum CBLAS_TRANSPOSE TransA, const int M, const int N,
|
||||
const double alpha, const double *A, const int lda,
|
||||
const double *X, const int incX, const double beta,
|
||||
double *Y, const int incY);
|
||||
void cblas_dgbmv(const enum CBLAS_ORDER order,
|
||||
const enum CBLAS_TRANSPOSE TransA, const int M, const int N,
|
||||
const int KL, const int KU, const double alpha,
|
||||
const double *A, const int lda, const double *X,
|
||||
const int incX, const double beta, double *Y, const int incY);
|
||||
void cblas_dtrmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
|
||||
const int N, const double *A, const int lda,
|
||||
double *X, const int incX);
|
||||
void cblas_dtbmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
|
||||
const int N, const int K, const double *A, const int lda,
|
||||
double *X, const int incX);
|
||||
void cblas_dtpmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
|
||||
const int N, const double *Ap, double *X, const int incX);
|
||||
void cblas_dtrsv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
|
||||
const int N, const double *A, const int lda, double *X,
|
||||
const int incX);
|
||||
void cblas_dtbsv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
|
||||
const int N, const int K, const double *A, const int lda,
|
||||
double *X, const int incX);
|
||||
void cblas_dtpsv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
|
||||
const int N, const double *Ap, double *X, const int incX);
|
||||
|
||||
void cblas_cgemv(const enum CBLAS_ORDER order,
|
||||
const enum CBLAS_TRANSPOSE TransA, const int M, const int N,
|
||||
const void *alpha, const void *A, const int lda,
|
||||
const void *X, const int incX, const void *beta,
|
||||
void *Y, const int incY);
|
||||
void cblas_cgbmv(const enum CBLAS_ORDER order,
|
||||
const enum CBLAS_TRANSPOSE TransA, const int M, const int N,
|
||||
const int KL, const int KU, const void *alpha,
|
||||
const void *A, const int lda, const void *X,
|
||||
const int incX, const void *beta, void *Y, const int incY);
|
||||
void cblas_ctrmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
|
||||
const int N, const void *A, const int lda,
|
||||
void *X, const int incX);
|
||||
void cblas_ctbmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
|
||||
const int N, const int K, const void *A, const int lda,
|
||||
void *X, const int incX);
|
||||
void cblas_ctpmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
|
||||
const int N, const void *Ap, void *X, const int incX);
|
||||
void cblas_ctrsv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
|
||||
const int N, const void *A, const int lda, void *X,
|
||||
const int incX);
|
||||
void cblas_ctbsv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
|
||||
const int N, const int K, const void *A, const int lda,
|
||||
void *X, const int incX);
|
||||
void cblas_ctpsv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
|
||||
const int N, const void *Ap, void *X, const int incX);
|
||||
|
||||
void cblas_zgemv(const enum CBLAS_ORDER order,
|
||||
const enum CBLAS_TRANSPOSE TransA, const int M, const int N,
|
||||
const void *alpha, const void *A, const int lda,
|
||||
const void *X, const int incX, const void *beta,
|
||||
void *Y, const int incY);
|
||||
void cblas_zgbmv(const enum CBLAS_ORDER order,
|
||||
const enum CBLAS_TRANSPOSE TransA, const int M, const int N,
|
||||
const int KL, const int KU, const void *alpha,
|
||||
const void *A, const int lda, const void *X,
|
||||
const int incX, const void *beta, void *Y, const int incY);
|
||||
void cblas_ztrmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
|
||||
const int N, const void *A, const int lda,
|
||||
void *X, const int incX);
|
||||
void cblas_ztbmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
|
||||
const int N, const int K, const void *A, const int lda,
|
||||
void *X, const int incX);
|
||||
void cblas_ztpmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
|
||||
const int N, const void *Ap, void *X, const int incX);
|
||||
void cblas_ztrsv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
|
||||
const int N, const void *A, const int lda, void *X,
|
||||
const int incX);
|
||||
void cblas_ztbsv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
|
||||
const int N, const int K, const void *A, const int lda,
|
||||
void *X, const int incX);
|
||||
void cblas_ztpsv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
|
||||
const int N, const void *Ap, void *X, const int incX);
|
||||
|
||||
|
||||
/*
|
||||
* Routines with S and D prefixes only
|
||||
*/
|
||||
void cblas_ssymv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const int N, const float alpha, const float *A,
|
||||
const int lda, const float *X, const int incX,
|
||||
const float beta, float *Y, const int incY);
|
||||
void cblas_ssbmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const int N, const int K, const float alpha, const float *A,
|
||||
const int lda, const float *X, const int incX,
|
||||
const float beta, float *Y, const int incY);
|
||||
void cblas_sspmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const int N, const float alpha, const float *Ap,
|
||||
const float *X, const int incX,
|
||||
const float beta, float *Y, const int incY);
|
||||
void cblas_sger(const enum CBLAS_ORDER order, const int M, const int N,
|
||||
const float alpha, const float *X, const int incX,
|
||||
const float *Y, const int incY, float *A, const int lda);
|
||||
void cblas_ssyr(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const int N, const float alpha, const float *X,
|
||||
const int incX, float *A, const int lda);
|
||||
void cblas_sspr(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const int N, const float alpha, const float *X,
|
||||
const int incX, float *Ap);
|
||||
void cblas_ssyr2(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const int N, const float alpha, const float *X,
|
||||
const int incX, const float *Y, const int incY, float *A,
|
||||
const int lda);
|
||||
void cblas_sspr2(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const int N, const float alpha, const float *X,
|
||||
const int incX, const float *Y, const int incY, float *A);
|
||||
|
||||
void cblas_dsymv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const int N, const double alpha, const double *A,
|
||||
const int lda, const double *X, const int incX,
|
||||
const double beta, double *Y, const int incY);
|
||||
void cblas_dsbmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const int N, const int K, const double alpha, const double *A,
|
||||
const int lda, const double *X, const int incX,
|
||||
const double beta, double *Y, const int incY);
|
||||
void cblas_dspmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const int N, const double alpha, const double *Ap,
|
||||
const double *X, const int incX,
|
||||
const double beta, double *Y, const int incY);
|
||||
void cblas_dger(const enum CBLAS_ORDER order, const int M, const int N,
|
||||
const double alpha, const double *X, const int incX,
|
||||
const double *Y, const int incY, double *A, const int lda);
|
||||
void cblas_dsyr(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const int N, const double alpha, const double *X,
|
||||
const int incX, double *A, const int lda);
|
||||
void cblas_dspr(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const int N, const double alpha, const double *X,
|
||||
const int incX, double *Ap);
|
||||
void cblas_dsyr2(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const int N, const double alpha, const double *X,
|
||||
const int incX, const double *Y, const int incY, double *A,
|
||||
const int lda);
|
||||
void cblas_dspr2(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const int N, const double alpha, const double *X,
|
||||
const int incX, const double *Y, const int incY, double *A);
|
||||
|
||||
|
||||
/*
|
||||
* Routines with C and Z prefixes only
|
||||
*/
|
||||
void cblas_chemv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const int N, const void *alpha, const void *A,
|
||||
const int lda, const void *X, const int incX,
|
||||
const void *beta, void *Y, const int incY);
|
||||
void cblas_chbmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const int N, const int K, const void *alpha, const void *A,
|
||||
const int lda, const void *X, const int incX,
|
||||
const void *beta, void *Y, const int incY);
|
||||
void cblas_chpmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const int N, const void *alpha, const void *Ap,
|
||||
const void *X, const int incX,
|
||||
const void *beta, void *Y, const int incY);
|
||||
void cblas_cgeru(const enum CBLAS_ORDER order, const int M, const int N,
|
||||
const void *alpha, const void *X, const int incX,
|
||||
const void *Y, const int incY, void *A, const int lda);
|
||||
void cblas_cgerc(const enum CBLAS_ORDER order, const int M, const int N,
|
||||
const void *alpha, const void *X, const int incX,
|
||||
const void *Y, const int incY, void *A, const int lda);
|
||||
void cblas_cher(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const int N, const float alpha, const void *X, const int incX,
|
||||
void *A, const int lda);
|
||||
void cblas_chpr(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const int N, const float alpha, const void *X,
|
||||
const int incX, void *A);
|
||||
void cblas_cher2(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo, const int N,
|
||||
const void *alpha, const void *X, const int incX,
|
||||
const void *Y, const int incY, void *A, const int lda);
|
||||
void cblas_chpr2(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo, const int N,
|
||||
const void *alpha, const void *X, const int incX,
|
||||
const void *Y, const int incY, void *Ap);
|
||||
|
||||
void cblas_zhemv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const int N, const void *alpha, const void *A,
|
||||
const int lda, const void *X, const int incX,
|
||||
const void *beta, void *Y, const int incY);
|
||||
void cblas_zhbmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const int N, const int K, const void *alpha, const void *A,
|
||||
const int lda, const void *X, const int incX,
|
||||
const void *beta, void *Y, const int incY);
|
||||
void cblas_zhpmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const int N, const void *alpha, const void *Ap,
|
||||
const void *X, const int incX,
|
||||
const void *beta, void *Y, const int incY);
|
||||
void cblas_zgeru(const enum CBLAS_ORDER order, const int M, const int N,
|
||||
const void *alpha, const void *X, const int incX,
|
||||
const void *Y, const int incY, void *A, const int lda);
|
||||
void cblas_zgerc(const enum CBLAS_ORDER order, const int M, const int N,
|
||||
const void *alpha, const void *X, const int incX,
|
||||
const void *Y, const int incY, void *A, const int lda);
|
||||
void cblas_zher(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const int N, const double alpha, const void *X, const int incX,
|
||||
void *A, const int lda);
|
||||
void cblas_zhpr(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
|
||||
const int N, const double alpha, const void *X,
|
||||
const int incX, void *A);
|
||||
void cblas_zher2(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo, const int N,
|
||||
const void *alpha, const void *X, const int incX,
|
||||
const void *Y, const int incY, void *A, const int lda);
|
||||
void cblas_zhpr2(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo, const int N,
|
||||
const void *alpha, const void *X, const int incX,
|
||||
const void *Y, const int incY, void *Ap);
|
||||
|
||||
/*
|
||||
* ===========================================================================
|
||||
* Prototypes for level 3 BLAS
|
||||
* ===========================================================================
|
||||
*/
|
||||
|
||||
/*
|
||||
* Routines with standard 4 prefixes (S, D, C, Z)
|
||||
*/
|
||||
void cblas_sgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA,
|
||||
const enum CBLAS_TRANSPOSE TransB, const int M, const int N,
|
||||
const int K, const float alpha, const float *A,
|
||||
const int lda, const float *B, const int ldb,
|
||||
const float beta, float *C, const int ldc);
|
||||
void cblas_ssymm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side,
|
||||
const enum CBLAS_UPLO Uplo, const int M, const int N,
|
||||
const float alpha, const float *A, const int lda,
|
||||
const float *B, const int ldb, const float beta,
|
||||
float *C, const int ldc);
|
||||
void cblas_ssyrk(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE Trans, const int N, const int K,
|
||||
const float alpha, const float *A, const int lda,
|
||||
const float beta, float *C, const int ldc);
|
||||
void cblas_ssyr2k(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE Trans, const int N, const int K,
|
||||
const float alpha, const float *A, const int lda,
|
||||
const float *B, const int ldb, const float beta,
|
||||
float *C, const int ldc);
|
||||
void cblas_strmm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side,
|
||||
const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE TransA,
|
||||
const enum CBLAS_DIAG Diag, const int M, const int N,
|
||||
const float alpha, const float *A, const int lda,
|
||||
float *B, const int ldb);
|
||||
void cblas_strsm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side,
|
||||
const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE TransA,
|
||||
const enum CBLAS_DIAG Diag, const int M, const int N,
|
||||
const float alpha, const float *A, const int lda,
|
||||
float *B, const int ldb);
|
||||
|
||||
void cblas_dgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA,
|
||||
const enum CBLAS_TRANSPOSE TransB, const int M, const int N,
|
||||
const int K, const double alpha, const double *A,
|
||||
const int lda, const double *B, const int ldb,
|
||||
const double beta, double *C, const int ldc);
|
||||
void cblas_dsymm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side,
|
||||
const enum CBLAS_UPLO Uplo, const int M, const int N,
|
||||
const double alpha, const double *A, const int lda,
|
||||
const double *B, const int ldb, const double beta,
|
||||
double *C, const int ldc);
|
||||
void cblas_dsyrk(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE Trans, const int N, const int K,
|
||||
const double alpha, const double *A, const int lda,
|
||||
const double beta, double *C, const int ldc);
|
||||
void cblas_dsyr2k(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE Trans, const int N, const int K,
|
||||
const double alpha, const double *A, const int lda,
|
||||
const double *B, const int ldb, const double beta,
|
||||
double *C, const int ldc);
|
||||
void cblas_dtrmm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side,
|
||||
const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE TransA,
|
||||
const enum CBLAS_DIAG Diag, const int M, const int N,
|
||||
const double alpha, const double *A, const int lda,
|
||||
double *B, const int ldb);
|
||||
void cblas_dtrsm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side,
|
||||
const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE TransA,
|
||||
const enum CBLAS_DIAG Diag, const int M, const int N,
|
||||
const double alpha, const double *A, const int lda,
|
||||
double *B, const int ldb);
|
||||
|
||||
void cblas_cgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA,
|
||||
const enum CBLAS_TRANSPOSE TransB, const int M, const int N,
|
||||
const int K, const void *alpha, const void *A,
|
||||
const int lda, const void *B, const int ldb,
|
||||
const void *beta, void *C, const int ldc);
|
||||
void cblas_csymm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side,
|
||||
const enum CBLAS_UPLO Uplo, const int M, const int N,
|
||||
const void *alpha, const void *A, const int lda,
|
||||
const void *B, const int ldb, const void *beta,
|
||||
void *C, const int ldc);
|
||||
void cblas_csyrk(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE Trans, const int N, const int K,
|
||||
const void *alpha, const void *A, const int lda,
|
||||
const void *beta, void *C, const int ldc);
|
||||
void cblas_csyr2k(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE Trans, const int N, const int K,
|
||||
const void *alpha, const void *A, const int lda,
|
||||
const void *B, const int ldb, const void *beta,
|
||||
void *C, const int ldc);
|
||||
void cblas_ctrmm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side,
|
||||
const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE TransA,
|
||||
const enum CBLAS_DIAG Diag, const int M, const int N,
|
||||
const void *alpha, const void *A, const int lda,
|
||||
void *B, const int ldb);
|
||||
void cblas_ctrsm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side,
|
||||
const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE TransA,
|
||||
const enum CBLAS_DIAG Diag, const int M, const int N,
|
||||
const void *alpha, const void *A, const int lda,
|
||||
void *B, const int ldb);
|
||||
|
||||
void cblas_zgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA,
|
||||
const enum CBLAS_TRANSPOSE TransB, const int M, const int N,
|
||||
const int K, const void *alpha, const void *A,
|
||||
const int lda, const void *B, const int ldb,
|
||||
const void *beta, void *C, const int ldc);
|
||||
void cblas_zsymm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side,
|
||||
const enum CBLAS_UPLO Uplo, const int M, const int N,
|
||||
const void *alpha, const void *A, const int lda,
|
||||
const void *B, const int ldb, const void *beta,
|
||||
void *C, const int ldc);
|
||||
void cblas_zsyrk(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE Trans, const int N, const int K,
|
||||
const void *alpha, const void *A, const int lda,
|
||||
const void *beta, void *C, const int ldc);
|
||||
void cblas_zsyr2k(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE Trans, const int N, const int K,
|
||||
const void *alpha, const void *A, const int lda,
|
||||
const void *B, const int ldb, const void *beta,
|
||||
void *C, const int ldc);
|
||||
void cblas_ztrmm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side,
|
||||
const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE TransA,
|
||||
const enum CBLAS_DIAG Diag, const int M, const int N,
|
||||
const void *alpha, const void *A, const int lda,
|
||||
void *B, const int ldb);
|
||||
void cblas_ztrsm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side,
|
||||
const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE TransA,
|
||||
const enum CBLAS_DIAG Diag, const int M, const int N,
|
||||
const void *alpha, const void *A, const int lda,
|
||||
void *B, const int ldb);
|
||||
|
||||
|
||||
/*
|
||||
* Routines with prefixes C and Z only
|
||||
*/
|
||||
void cblas_chemm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side,
|
||||
const enum CBLAS_UPLO Uplo, const int M, const int N,
|
||||
const void *alpha, const void *A, const int lda,
|
||||
const void *B, const int ldb, const void *beta,
|
||||
void *C, const int ldc);
|
||||
void cblas_cherk(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE Trans, const int N, const int K,
|
||||
const float alpha, const void *A, const int lda,
|
||||
const float beta, void *C, const int ldc);
|
||||
void cblas_cher2k(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE Trans, const int N, const int K,
|
||||
const void *alpha, const void *A, const int lda,
|
||||
const void *B, const int ldb, const float beta,
|
||||
void *C, const int ldc);
|
||||
|
||||
void cblas_zhemm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side,
|
||||
const enum CBLAS_UPLO Uplo, const int M, const int N,
|
||||
const void *alpha, const void *A, const int lda,
|
||||
const void *B, const int ldb, const void *beta,
|
||||
void *C, const int ldc);
|
||||
void cblas_zherk(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE Trans, const int N, const int K,
|
||||
const double alpha, const void *A, const int lda,
|
||||
const double beta, void *C, const int ldc);
|
||||
void cblas_zher2k(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo,
|
||||
const enum CBLAS_TRANSPOSE Trans, const int N, const int K,
|
||||
const void *alpha, const void *A, const int lda,
|
||||
const void *B, const int ldb, const double beta,
|
||||
void *C, const int ldc);
|
||||
|
||||
void cblas_xerbla(int p, const char *rout, const char *form, ...);
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
||||
#endif
|
||||
File diff suppressed because it is too large
Load diff
|
|
@ -50,10 +50,6 @@ def build(bld):
|
|||
dsp/rateconversion/Decimator.cpp
|
||||
dsp/rateconversion/DecimatorB.cpp
|
||||
dsp/rhythm/BeatSpectrum.cpp
|
||||
dsp/segmentation/cluster_melt.c
|
||||
dsp/segmentation/ClusterMeltSegmenter.cpp
|
||||
dsp/segmentation/cluster_segmenter.c
|
||||
dsp/segmentation/Segmenter.cpp
|
||||
dsp/signalconditioning/DFProcess.cpp
|
||||
dsp/signalconditioning/Filter.cpp
|
||||
dsp/signalconditioning/FiltFilt.cpp
|
||||
|
|
@ -66,7 +62,6 @@ def build(bld):
|
|||
dsp/tonal/TonalEstimator.cpp
|
||||
dsp/transforms/FFT.cpp
|
||||
dsp/wavelet/Wavelet.cpp
|
||||
hmm/hmm.c
|
||||
maths/Correlation.cpp
|
||||
maths/CosineDistance.cpp
|
||||
maths/KLDivergence.cpp
|
||||
|
|
@ -77,7 +72,7 @@ def build(bld):
|
|||
'''
|
||||
autowaf.ensure_visible_symbols (obj, True)
|
||||
obj.export_includes = ['.']
|
||||
obj.includes = ['.', 'include/', 'ext/kissfft', 'ext/kissfft/tools/']
|
||||
obj.includes = ['.', 'ext/kissfft', 'ext/kissfft/tools/']
|
||||
obj.defines = ['kiss_fft_scalar=double']
|
||||
obj.name = 'libqm-dsp'
|
||||
obj.target = 'qm-dsp'
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue