comparison src/EM.cpp @ 45:e92376d450b0 preshift

Get notes out using pre-shifted templates (not working properly)
author Chris Cannam
date Mon, 07 Apr 2014 13:01:08 +0100
parents c0c4a945577a
children ce1d88759557
comparison
equal deleted inserted replaced
44:eec530c4300d 45:e92376d450b0
29 using std::endl; 29 using std::endl;
30 30
31 static double epsilon = 1e-16; 31 static double epsilon = 1e-16;
32 32
33 EM::EM() : 33 EM::EM() :
34 m_notes(SILVET_TEMPLATE_NOTE_COUNT), 34 m_noteCount(SILVET_TEMPLATE_NOTE_COUNT),
35 m_bins(SILVET_TEMPLATE_HEIGHT), 35 m_shiftCount(SILVET_TEMPLATE_MAX_SHIFT * 2 + 1),
36 m_instruments(SILVET_TEMPLATE_COUNT), 36 m_pitchCount(m_noteCount * m_shiftCount),
37 m_binCount(SILVET_TEMPLATE_HEIGHT),
38 m_instrumentCount(SILVET_TEMPLATE_COUNT),
37 m_pitchSparsity(1.1), 39 m_pitchSparsity(1.1),
38 m_sourceSparsity(1.3) 40 m_sourceSparsity(1.3)
39 { 41 {
40 m_lowest = 0; 42 m_lowestPitch =
41 m_highest = m_notes - 1; 43 silvet_templates_lowest_note * m_shiftCount;
42 44 m_highestPitch =
43 for (int i = 0; i < m_instruments; ++i) { 45 silvet_templates_highest_note * m_shiftCount + m_shiftCount - 1;
44 if (i == 0 || silvet_templates[i].lowest < m_lowest) { 46
45 m_lowest = silvet_templates[i].lowest; 47 m_pitches = V(m_pitchCount);
46 } 48
47 if (i == 0 || silvet_templates[i].highest > m_highest) { 49 for (int n = 0; n < m_pitchCount; ++n) {
48 m_highest = silvet_templates[i].highest;
49 }
50 }
51
52 m_pitches = V(m_notes);
53
54 for (int n = 0; n < m_notes; ++n) {
55 m_pitches[n] = drand48(); 50 m_pitches[n] = drand48();
56 } 51 }
57 52
58 m_sources = Grid(m_instruments); 53 m_sources = Grid(m_instrumentCount);
59 54
60 for (int i = 0; i < m_instruments; ++i) { 55 for (int i = 0; i < m_instrumentCount; ++i) {
61 m_sources[i] = V(m_notes); 56 m_sources[i] = V(m_pitchCount);
62 for (int n = 0; n < m_notes; ++n) { 57 for (int n = 0; n < m_pitchCount; ++n) {
63 m_sources[i][n] = (inRange(i, n) ? 1.0 : 0.0); 58 m_sources[i][n] = (inRange(i, n) ? 1.0 : 0.0);
64 } 59 }
65 } 60 }
66 61
67 m_estimate = V(m_bins); 62 m_estimate = V(m_binCount);
68 m_q = V(m_bins); 63 m_q = V(m_binCount);
69 } 64 }
70 65
71 EM::~EM() 66 EM::~EM()
72 { 67 {
73 } 68 }
74 69
70 void
71 EM::rangeFor(int instrument, int &minPitch, int &maxPitch)
72 {
73 minPitch = silvet_templates[instrument].lowest * m_shiftCount;
74 maxPitch = silvet_templates[instrument].highest * m_shiftCount
75 + m_shiftCount - 1;
76 }
77
75 bool 78 bool
76 EM::inRange(int instrument, int note) 79 EM::inRange(int instrument, int pitch)
77 { 80 {
78 return (note >= silvet_templates[instrument].lowest && 81 int minPitch, maxPitch;
79 note <= silvet_templates[instrument].highest); 82 rangeFor(instrument, minPitch, maxPitch);
83 return (pitch >= minPitch && pitch <= maxPitch);
80 } 84 }
81 85
82 void 86 void
83 EM::normalise(V &column) 87 EM::normalise(V &column)
84 { 88 {
97 normalise(column); 101 normalise(column);
98 expectation(column); 102 expectation(column);
99 maximisation(column); 103 maximisation(column);
100 } 104 }
101 105
106 const float *
107 EM::templateFor(int instrument, int pitch)
108 {
109 int note = pitch / m_shiftCount;
110 int shift = pitch % m_shiftCount;
111 return silvet_templates[instrument].data[note] + shift;
112 }
113
102 void 114 void
103 EM::expectation(const V &column) 115 EM::expectation(const V &column)
104 { 116 {
105 cerr << "."; 117 cerr << ".";
106 118
107 for (int i = 0; i < m_bins; ++i) { 119 for (int i = 0; i < m_binCount; ++i) {
108 m_estimate[i] = epsilon; 120 m_estimate[i] = epsilon;
109 } 121 }
110 122
111 for (int i = 0; i < m_instruments; ++i) { 123 for (int i = 0; i < m_instrumentCount; ++i) {
112 for (int n = 0; n < m_notes; ++n) { 124 for (int n = 0; n < m_pitchCount; ++n) {
113 float *w = silvet_templates[i].data[n]; 125 const float *w = templateFor(i, n);
114 double pitch = m_pitches[n]; 126 double pitch = m_pitches[n];
115 double source = m_sources[i][n]; 127 double source = m_sources[i][n];
116 for (int j = 0; j < m_bins; ++j) { 128 for (int j = 0; j < m_binCount; ++j) {
117 m_estimate[j] += w[j] * pitch * source; 129 m_estimate[j] += w[j] * pitch * source;
118 } 130 }
119 } 131 }
120 } 132 }
121 133
122 for (int i = 0; i < m_bins; ++i) { 134 for (int i = 0; i < m_binCount; ++i) {
123 m_q[i] = column[i] / m_estimate[i]; 135 m_q[i] = column[i] / m_estimate[i];
124 } 136 }
125 } 137 }
126 138
127 void 139 void
128 EM::maximisation(const V &column) 140 EM::maximisation(const V &column)
129 { 141 {
130 V newPitches = m_pitches; 142 V newPitches = m_pitches;
131 143
132 for (int n = 0; n < m_notes; ++n) { 144 for (int n = 0; n < m_pitchCount; ++n) {
133 newPitches[n] = epsilon; 145 newPitches[n] = epsilon;
134 if (n >= m_lowest && n <= m_highest) { 146 if (n >= m_lowestPitch && n <= m_highestPitch) {
135 for (int i = 0; i < m_instruments; ++i) { 147 for (int i = 0; i < m_instrumentCount; ++i) {
136 float *w = silvet_templates[i].data[n]; 148 const float *w = templateFor(i, n);
137 double pitch = m_pitches[n]; 149 double pitch = m_pitches[n];
138 double source = m_sources[i][n]; 150 double source = m_sources[i][n];
139 for (int j = 0; j < m_bins; ++j) { 151 for (int j = 0; j < m_binCount; ++j) {
140 newPitches[n] += w[j] * m_q[j] * pitch * source; 152 newPitches[n] += w[j] * m_q[j] * pitch * source;
141 } 153 }
142 } 154 }
143 } 155 }
144 if (m_pitchSparsity != 1.0) { 156 if (m_pitchSparsity != 1.0) {
147 } 159 }
148 normalise(newPitches); 160 normalise(newPitches);
149 161
150 Grid newSources = m_sources; 162 Grid newSources = m_sources;
151 163
152 for (int i = 0; i < m_instruments; ++i) { 164 for (int i = 0; i < m_instrumentCount; ++i) {
153 for (int n = 0; n < m_notes; ++n) { 165 for (int n = 0; n < m_pitchCount; ++n) {
154 newSources[i][n] = epsilon; 166 newSources[i][n] = epsilon;
155 if (inRange(i, n)) { 167 if (inRange(i, n)) {
156 float *w = silvet_templates[i].data[n]; 168 const float *w = templateFor(i, n);
157 double pitch = m_pitches[n]; 169 double pitch = m_pitches[n];
158 double source = m_sources[i][n]; 170 double source = m_sources[i][n];
159 for (int j = 0; j < m_bins; ++j) { 171 for (int j = 0; j < m_binCount; ++j) {
160 newSources[i][n] += w[j] * m_q[j] * pitch * source; 172 newSources[i][n] += w[j] * m_q[j] * pitch * source;
161 } 173 }
162 } 174 }
163 if (m_sourceSparsity != 1.0) { 175 if (m_sourceSparsity != 1.0) {
164 newSources[i][n] = pow(newSources[i][n], m_sourceSparsity); 176 newSources[i][n] = pow(newSources[i][n], m_sourceSparsity);
173 185
174 void 186 void
175 EM::report() 187 EM::report()
176 { 188 {
177 vector<int> sounding; 189 vector<int> sounding;
178 for (int n = 0; n < m_notes; ++n) { 190 for (int n = 0; n < m_pitchCount; ++n) {
179 if (m_pitches[n] > 0.05) { 191 if (m_pitches[n] > 0.05) {
180 sounding.push_back(n); 192 sounding.push_back(n);
181 } 193 }
182 } 194 }
183 cerr << " sounding: "; 195 cerr << " sounding: ";
184 for (int i = 0; i < (int)sounding.size(); ++i) { 196 for (int i = 0; i < (int)sounding.size(); ++i) {
185 cerr << sounding[i] << " "; 197 cerr << sounding[i] << " ";
186 int maxj = -1; 198 int maxj = -1;
187 double maxs = 0.0; 199 double maxs = 0.0;
188 for (int j = 0; j < m_instruments; ++j) { 200 for (int j = 0; j < m_instrumentCount; ++j) {
189 if (j == 0 || m_sources[j][sounding[i]] > maxs) { 201 if (j == 0 || m_sources[j][sounding[i]] > maxs) {
190 maxj = j; 202 maxj = j;
191 maxs = m_sources[j][sounding[i]]; 203 maxs = m_sources[j][sounding[i]];
192 } 204 }
193 } 205 }