← Back to branch summary

~ubuntu-branches/ubuntu/utopic/ardour3/utopic

« back to all changes in this revision

Viewing changes to libs/qm-dsp/dsp/segmentation/cluster_segmenter.c

  • Committer: Package Import Robot
  • Author(s): Felipe Sateler
  • Date: 2013-09-21 19:05:02 UTC
  • Revision ID: package-import@ubuntu.com-20130921190502-8gsftrku6jnzhd7v
Tags: upstream-3.4~dfsg
Import upstream version 3.4~dfsg

Show diffs side-by-side

added added

removed removed

Lines of Context:
libs/qm-dsp/dsp/segmentation/cluster_segmenter.c
 
1
/*
 
2
 *  cluster_segmenter.c
 
3
 *  soundbite
 
4
 *
 
5
 *  Created by Mark Levy on 06/04/2006.
 
6
 *  Copyright 2006 Centre for Digital Music, Queen Mary, University of London.
 
7
 
 
8
    This program is free software; you can redistribute it and/or
 
9
    modify it under the terms of the GNU General Public License as
 
10
    published by the Free Software Foundation; either version 2 of the
 
11
    License, or (at your option) any later version.  See the file
 
12
    COPYING included with this distribution for more information.
 
13
 *
 
14
 */
 
15
 
 
16
#include "cluster_segmenter.h"
 
17
 
 
18
extern int readmatarray_size(const char *filepath, int n_array, int* t, int* d);
 
19
extern int readmatarray(const char *filepath, int n_array, int t, int d, double** arr);
 
20
 
 
21
/* converts constant-Q features to normalised chroma */
 
22
void cq2chroma(double** cq, int nframes, int ncoeff, int bins, double** chroma)
 
23
{
 
24
        int noct = ncoeff / bins;       /* number of complete octaves in constant-Q */
 
25
        int t, b, oct, ix;
 
26
        //double maxchroma;     /* max chroma value at each time, for normalisation */
 
27
        //double sum;           /* for normalisation */
 
28
        
 
29
        for (t = 0; t < nframes; t++)
 
30
        {
 
31
                for (b = 0; b < bins; b++)
 
32
                        chroma[t][b] = 0;
 
33
                for (oct = 0; oct < noct; oct++)
 
34
                {
 
35
                        ix = oct * bins;
 
36
                        for (b = 0; b < bins; b++)
 
37
                                chroma[t][b] += fabs(cq[t][ix+b]);
 
38
                }
 
39
                /* normalise to unit sum
 
40
                sum = 0;
 
41
                for (b = 0; b < bins; b++)
 
42
                        sum += chroma[t][b];
 
43
                for (b = 0; b < bins; b++)
 
44
                        chroma[t][b] /= sum;
 
45
                */
 
46
                /* normalise to unit max - NO this made results much worse!
 
47
                maxchroma = 0;
 
48
                for (b = 0; b < bins; b++)
 
49
                        if (chroma[t][b] > maxchroma)
 
50
                                maxchroma = chroma[t][b];
 
51
                if (maxchroma > 0)
 
52
                        for (b = 0; b < bins; b++)
 
53
                                chroma[t][b] /= maxchroma;      
 
54
                */
 
55
        }
 
56
}
 
57
 
 
58
/* applies MPEG-7 normalisation to constant-Q features, storing normalised envelope (norm) in last feature dimension */
 
59
void mpeg7_constq(double** features, int nframes, int ncoeff)
 
60
{
 
61
        int i, j;
 
62
        double ss;
 
63
        double env;
 
64
        double maxenv = 0;
 
65
        
 
66
        /* convert const-Q features to dB scale */
 
67
        for (i = 0; i < nframes; i++)
 
68
                for (j = 0; j < ncoeff; j++)
 
69
                        features[i][j] = 10.0 * log10(features[i][j]+DBL_EPSILON);
 
70
        
 
71
        /* normalise each feature vector and add the norm as an extra feature dimension */      
 
72
        for (i = 0; i < nframes; i++)
 
73
        {
 
74
                ss = 0;
 
75
                for (j = 0; j < ncoeff; j++)
 
76
                        ss += features[i][j] * features[i][j];
 
77
                env = sqrt(ss);
 
78
                for (j = 0; j < ncoeff; j++)
 
79
                        features[i][j] /= env;
 
80
                features[i][ncoeff] = env;
 
81
                if (env > maxenv)
 
82
                        maxenv = env;
 
83
        } 
 
84
        /* normalise the envelopes */
 
85
        for (i = 0; i < nframes; i++)
 
86
                features[i][ncoeff] /= maxenv;  
 
87
}
 
88
 
 
89
/* return histograms h[nx*m] of data x[nx] into m bins using a sliding window of length h_len (MUST BE ODD) */
 
90
/* NB h is a vector in row major order, as required by cluster_melt() */
 
91
/* for historical reasons we normalise the histograms by their norm (not to sum to one) */
 
92
void create_histograms(int* x, int nx, int m, int hlen, double* h)
 
93
{
 
94
        int i, j, t;
 
95
        double norm;
 
96
 
 
97
        for (i = 0; i < nx*m; i++) 
 
98
                h[i] = 0;
 
99
 
 
100
        for (i = hlen/2; i < nx-hlen/2; i++)
 
101
        {
 
102
                for (j = 0; j < m; j++)
 
103
                        h[i*m+j] = 0;
 
104
                for (t = i-hlen/2; t <= i+hlen/2; t++)
 
105
                        ++h[i*m+x[t]];
 
106
                norm = 0;
 
107
                for (j = 0; j < m; j++)
 
108
                        norm += h[i*m+j] * h[i*m+j];
 
109
                for (j = 0; j < m; j++)
 
110
                        h[i*m+j] /= norm;
 
111
        }
 
112
        
 
113
        /* duplicate histograms at beginning and end to create one histogram for each data value supplied */
 
114
        for (i = 0; i < hlen/2; i++)
 
115
                for (j = 0; j < m; j++)
 
116
                        h[i*m+j] = h[hlen/2*m+j];
 
117
        for (i = nx-hlen/2; i < nx; i++)
 
118
                for (j = 0; j < m; j++)
 
119
                        h[i*m+j] = h[(nx-hlen/2-1)*m+j];
 
120
}
 
121
 
 
122
/* segment using HMM and then histogram clustering */
 
123
void cluster_segment(int* q, double** features, int frames_read, int feature_length, int nHMM_states, 
 
124
                                         int histogram_length, int nclusters, int neighbour_limit)
 
125
{
 
126
        int i, j;
 
127
        
 
128
        /*****************************/
 
129
        if (0) {
 
130
        /* try just using the predominant bin number as a 'decoded state' */
 
131
        nHMM_states = feature_length + 1;       /* allow a 'zero' state */
 
132
        double chroma_thresh = 0.05;
 
133
        double maxval;
 
134
        int maxbin;
 
135
        for (i = 0; i < frames_read; i++)
 
136
        {
 
137
                maxval = 0;
 
138
                for (j = 0; j < feature_length; j++)
 
139
                {
 
140
                        if (features[i][j] > maxval) 
 
141
                        {
 
142
                                maxval = features[i][j];
 
143
                                maxbin = j;
 
144
                        }                               
 
145
                }
 
146
                if (maxval > chroma_thresh)
 
147
                        q[i] = maxbin;
 
148
                else
 
149
                        q[i] = feature_length;
 
150
        }
 
151
        
 
152
        }
 
153
        if (1) {
 
154
        /*****************************/
 
155
                
 
156
        
 
157
        /* scale all the features to 'balance covariances' during HMM training */
 
158
        double scale = 10;
 
159
        for (i = 0; i < frames_read; i++)
 
160
                for (j = 0; j < feature_length; j++)
 
161
                        features[i][j] *= scale;
 
162
        
 
163
        /* train an HMM on the features */
 
164
        
 
165
        /* create a model */
 
166
        model_t* model = hmm_init(features, frames_read, feature_length, nHMM_states);
 
167
        
 
168
        /* train the model */
 
169
        hmm_train(features, frames_read, model);
 
170
/*      
 
171
        printf("\n\nafter training:\n");
 
172
        hmm_print(model);
 
173
*/      
 
174
        /* decode the hidden state sequence */
 
175
        viterbi_decode(features, frames_read, model, q);  
 
176
        hmm_close(model);
 
177
        
 
178
        /*****************************/
 
179
        }
 
180
        /*****************************/
 
181
        
 
182
    
 
183
/*
 
184
        fprintf(stderr, "HMM state sequence:\n");
 
185
        for (i = 0; i < frames_read; i++)
 
186
                fprintf(stderr, "%d ", q[i]);
 
187
        fprintf(stderr, "\n\n");
 
188
*/
 
189
        
 
190
        /* create histograms of states */
 
191
        double* h = (double*) malloc(frames_read*nHMM_states*sizeof(double));   /* vector in row major order */
 
192
        create_histograms(q, frames_read, nHMM_states, histogram_length, h);
 
193
        
 
194
        /* cluster the histograms */
 
195
        int nbsched = 20;       /* length of inverse temperature schedule */
 
196
        double* bsched = (double*) malloc(nbsched*sizeof(double));      /* inverse temperature schedule */
 
197
        double b0 = 100;
 
198
        double alpha = 0.7;
 
199
        bsched[0] = b0;
 
200
        for (i = 1; i < nbsched; i++)
 
201
                bsched[i] = alpha * bsched[i-1];
 
202
        cluster_melt(h, nHMM_states, frames_read, bsched, nbsched, nclusters, neighbour_limit, q);
 
203
        
 
204
        /* now q holds a sequence of cluster assignments */
 
205
        
 
206
        free(h);  
 
207
        free(bsched);
 
208
}
 
209
 
 
210
/* segment constant-Q or chroma features */
 
211
void constq_segment(int* q, double** features, int frames_read, int bins, int ncoeff, int feature_type, 
 
212
                         int nHMM_states, int histogram_length, int nclusters, int neighbour_limit)
 
213
{
 
214
        int feature_length;
 
215
        double** chroma;
 
216
        int i;
 
217
        
 
218
        if (feature_type == FEATURE_TYPE_CONSTQ)
 
219
        {
 
220
/*              fprintf(stderr, "Converting to dB and normalising...\n");
 
221
 */             
 
222
                mpeg7_constq(features, frames_read, ncoeff);
 
223
/*              
 
224
                fprintf(stderr, "Running PCA...\n");
 
225
*/              
 
226
                /* do PCA on the features (but not the envelope) */
 
227
                int ncomponents = 20;
 
228
                pca_project(features, frames_read, ncoeff, ncomponents);
 
229
                
 
230
                /* copy the envelope so that it immediatly follows the chosen components */
 
231
                for (i = 0; i < frames_read; i++)
 
232
                        features[i][ncomponents] = features[i][ncoeff]; 
 
233
                
 
234
                feature_length = ncomponents + 1;
 
235
                
 
236
                /**************************************
 
237
                //TEST
 
238
                // feature file name
 
239
                char* dir = "/Users/mark/documents/semma/audio/";
 
240
                char* file_name = (char*) malloc((strlen(dir) + strlen(trackname) + strlen("_features_c20r8h0.2f0.6.mat") + 1)*sizeof(char));
 
241
                strcpy(file_name, dir);
 
242
                strcat(file_name, trackname);
 
243
                strcat(file_name, "_features_c20r8h0.2f0.6.mat");
 
244
                
 
245
                // get the features from Matlab from mat-file
 
246
                int frames_in_file;
 
247
                readmatarray_size(file_name, 2, &frames_in_file, &feature_length);
 
248
                readmatarray(file_name, 2, frames_in_file, feature_length, features);
 
249
                // copy final frame to ensure that we get as many as we expected
 
250
                int missing_frames = frames_read - frames_in_file;
 
251
                while (missing_frames > 0)
 
252
                {
 
253
                        for (i = 0; i < feature_length; i++)
 
254
                                features[frames_read-missing_frames][i] = features[frames_read-missing_frames-1][i];
 
255
                        --missing_frames;
 
256
                }
 
257
                
 
258
                free(file_name);
 
259
                ******************************************/
 
260
        
 
261
                cluster_segment(q, features, frames_read, feature_length, nHMM_states, histogram_length, nclusters, neighbour_limit);
 
262
        }
 
263
        
 
264
        if (feature_type == FEATURE_TYPE_CHROMA)
 
265
        {
 
266
/*
 
267
                fprintf(stderr, "Converting to chroma features...\n");
 
268
*/              
 
269
                /* convert constant-Q to normalised chroma features */
 
270
                chroma = (double**) malloc(frames_read*sizeof(double*));
 
271
                for (i = 0; i < frames_read; i++)
 
272
                        chroma[i] = (double*) malloc(bins*sizeof(double));
 
273
                cq2chroma(features, frames_read, ncoeff, bins, chroma);
 
274
                feature_length = bins;
 
275
                
 
276
                cluster_segment(q, chroma, frames_read, feature_length, nHMM_states, histogram_length, nclusters, neighbour_limit);
 
277
        
 
278
                for (i = 0; i < frames_read; i++)
 
279
                        free(chroma[i]);
 
280
                free(chroma);
 
281
        }
 
282
}
 
283
 
 
284
 
 
285