Entropy Piano Tuner  1.1.3 (documentation not yet complete)
An open-source experimental software for piano tuning by entropy minimization
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
auditorypreprocessing.cpp
Go to the documentation of this file.
1 /*****************************************************************************
2  * Copyright 2015 Haye Hinrichsen, Christoph Wick
3  *
4  * This file is part of Entropy Piano Tuner.
5  *
6  * Entropy Piano Tuner is free software: you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by the
8  * Free Software Foundation, either version 3 of the License, or (at your
9  * option) any later version.
10  *
11  * Entropy Piano Tuner is distributed in the hope that it will be useful, but
12  * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
13  * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
14  * more details.
15  *
16  * You should have received a copy of the GNU General Public License along with
17  * Entropy Piano Tuner. If not, see http://www.gnu.org/licenses/.
18  *****************************************************************************/
19 
20 //=============================================================================
21 // Auditory preprocessing
22 //=============================================================================
23 
24 #include "auditorypreprocessing.h"
25 
26 #include <iostream>
27 
31 #include "core/piano/keyboard.h"
32 #include "core/piano/key.h"
33 #include "core/math/mathtools.h"
34 #include "core/config.h"
35 
37 {
38 
39 //-----------------------------------------------------------------------------
40 // Constructor
41 //-----------------------------------------------------------------------------
42 
48 
50  mPiano(piano),
51  mKeyboard(mPiano.getKeyboard()),
52  mKeys(mKeyboard.getKeys()),
53  mNumberOfKeys(mKeyboard.getNumberOfKeys()),
54  mKeyNumberOfA4(mKeyboard.getKeyNumberOfA4()),
55  mdBA()
56 {
57 }
58 
59 
60 //-----------------------------------------------------------------------------
61 // Check consistency of dataset
62 //-----------------------------------------------------------------------------
63 
72 {
73  // First check whether the number of keys is consistent.
74  EptAssert(mKeys.size() > 0, "Piano should have at least one key");
75  EptAssert(mKeys.size() == static_cast<size_t>(mNumberOfKeys),
76  "Key vector length mismatch");
77 
78  // Next let us see whether all keys have been recorded.
79  bool allkeysrecorded = true;
80  bool somekeysrecorded = false;
81  for (Key key : mKeys)
82  {
83  if (key.isRecorded()) somekeysrecorded = true;
84  else allkeysrecorded = false;
85  }
86  if (not somekeysrecorded)
87  {
88  MessageHandler::send<MessageCaluclationProgress>
91  LogW("Calculation started without data");
92  return false;
93  }
94  if (not allkeysrecorded)
95  {
96  MessageHandler::send<MessageCaluclationProgress>
99  LogW("Not all keys have been recorded");
100  return false;
101  }
102 
103  // Thirdly, let us see whether the frequencies and inharmonicities
104  // are within reasonable bounds
105  bool consistent = true;
106  for (int keynumber=0; keynumber<mNumberOfKeys; keynumber++)
107  {
108  Key &key = mKeys[keynumber];
109  double f = key.getRecordedFrequency();
110  if (f<20 or f>20000)
111  {
112  LogW("Key %d: Frequency f=%f out of range.",keynumber,f);
113  consistent = false;
114  }
115  double B = key.getMeasuredInharmonicity();
116  if (B<0 or B>1)
117  {
118  LogW("Key %d: Inharmonicity B=%f out of range.",keynumber,B);
119  consistent=false;
120  }
121  Key::SpectrumType &spectrum = key.getSpectrum();
122  if (spectrum.size() != static_cast<size_t>(Key::NumberOfBins))
123  {
124  LogW("Key %d: Logspec size is %d, expected %d.",
125  keynumber,
126  static_cast<int>(spectrum.size()),
127  static_cast<int>(Key::NumberOfBins));
128  consistent=false;
129  }
130  else if (MathTools::computeNorm(spectrum)==0)
131  {
132  LogW("Key %d: Logspec norm = %f.",keynumber,
133  MathTools::computeNorm(spectrum));
134  consistent=false;
135  }
136  }
137 
138  // If data is found to be inconsistent send a corresponding message.
139  if (not consistent)
140  {
141  MessageHandler::send<MessageCaluclationProgress>
144  }
145  return consistent;
146 }
147 
148 
149 //-----------------------------------------------------------------------------
150 // Normalize all spectra
151 //-----------------------------------------------------------------------------
152 
159 
161 {
162  Key::SpectrumType &spectrum = key.getSpectrum();
163  MathTools::normalize(spectrum);
164 }
165 
166 
167 //-----------------------------------------------------------------------------
168 // Compute frequency of inharmonic partial
169 //-----------------------------------------------------------------------------
170 
181 
182 double AuditoryPreprocessing::getInharmonicPartial (double n, double f1, double B)
183 {
184  return f1 * n * sqrt((1+B*n*n)/(1+B));
185 }
186 
187 
188 //-----------------------------------------------------------------------------
189 // Compute the inharmonic index for a given frequency
190 //-----------------------------------------------------------------------------
191 
202 
203 double AuditoryPreprocessing::getInharmonicPartialIndex (double f, double f1, double B)
204 {
205  double x=f/f1;
206  if (B==0) return x;
207  else return 0.7071067811865475*sqrt((-1 + sqrt(1 + 4*B*(1+B)*x*x))/B);
208 }
209 
210 
211 //-----------------------------------------------------------------------------
212 // Clean spectrum
213 //-----------------------------------------------------------------------------
214 
229 
231 {
232  SpectrumType &spectrum = key.getSpectrum();
233  int M = spectrum.size();
234 
235  double f = key.getRecordedFrequency();
236  double B = key.getMeasuredInharmonicity();
237  auto wave = [f,B,this] (int m)
238  { return cos(MathTools::PI*getInharmonicPartialIndex(mtof(m),f,B)); };
239  auto envelope = [wave,f,this] (int m)
240  { return pow(fabs(wave(m)),200.0/pow(mtof(m)/f,1.5)); };
241  for (int m=0; m<M; m++) spectrum[m] *= envelope(m);
242 }
243 
244 
245 
246 //-----------------------------------------------------------------------------
247 // Cut frequencies below f1
248 //-----------------------------------------------------------------------------
249 
254 
256 {
257  SpectrumType &spectrum = key.getSpectrum();
258  double f = key.getRecordedFrequency();
259  int lowcutindex = std::min(static_cast<size_t>((5*ftom(f)))/6,NumberOfBins);
260  for (int m=0; m<lowcutindex; m++) spectrum[m]=0;
261 }
262 
263 
264 
265 //-----------------------------------------------------------------------------
266 // Initialize the filter function
267 //-----------------------------------------------------------------------------
268 
278 
280 {
281  mdBA.clear();
282  mdBA.resize(NumberOfBins);
283  for (uint m=0; m<NumberOfBins; ++m)
284  {
285  double f = mtof(m);
286  double Ra = 12200.0*12200.0*f*f*f*f / (f*f+20.6*20.6) /
287  sqrt((f*f+107.7*107.7)*(f*f+737.9*737.9)) / (f*f+12200.0*12200.0);
288  mdBA[m] = 2.0+20*log10(Ra);
289  //std::cout << f << "\t" << mdBA[m] << std::endl;
290  }
291 }
292 
293 
294 //-----------------------------------------------------------------------------
295 // Compute A-weighted sound pressure level SPLA
296 //-----------------------------------------------------------------------------
297 
309 
311 {
312  if (mdBA.size()==0) initializeSPLAFilter();
313  EptAssert(mdBA.size()==NumberOfBins,"mdBA should be initialized.");
314  const double I0 = 1E-7; // auditory threshold intensity
315  for (uint m=0; m<NumberOfBins; ++m)
316  {
317  double SPLA = 10 * log10(spectrum[m]/I0) + mdBA[m];
318  if (SPLA<0) spectrum[m]=0;
319  else spectrum[m] = I0 * pow(10.0, SPLA/10.0);
320  }
321 }
322 
323 
324 //-----------------------------------------------------------------------------
325 // extrapolate missing inharmonicity estimates
326 //-----------------------------------------------------------------------------
327 
339 
341 {
342  const int firstkey = mKeyNumberOfA4-8;
343  double K=0, Y=0, KK=0, KY=0, N=0, BE=0;
344  for (int k = firstkey; k < mNumberOfKeys; ++k)
345  {
346  if (N>1)
347  {
348  double a = (N*KY-K*Y)/(N*KK-K*K);
349  double b = (KK*Y-K*KY)/(N*KK-K*K);
350  BE = exp(a*k+b);
351  }
352  double B = mKeys[k].getMeasuredInharmonicity();
353  bool B_is_valid = (B>0);
354  if (B>0 and BE>0 and N>5) if (fabs(log(B/BE))>0.2) B_is_valid=false;
355  if (B_is_valid)
356  {
357  double y = log(B);
358  K+=k; Y+=y; KK+=k*k; KY+=k*y; N++;
359  }
360  else
361  {
362  if (BE==0)
363  {
364  double f1 = mPiano.getEqualTempFrequency(k,0,440);
366  }
367  mKeys[k].setMeasuredInharmonicity(BE);
368  }
369  }
370 }
371 
372 
373 //-----------------------------------------------------------------------------
374 // Improve high-frequency spectral lines
375 //-----------------------------------------------------------------------------
376 
383 
385 {
386  int start = mKeyNumberOfA4;
387  for (int k = start; k < mNumberOfKeys; k++)
388  {
389  Key &key = mKeys[k];
390  SpectrumType &spectrum = key.getSpectrum();
391  double f = key.getRecordedFrequency();
392  double B = key.getMeasuredInharmonicity();
393  if (f<=0 or B<=0) continue;
394  int m = MathTools::roundToInteger(ftom(f));
395  if (m<0 or m>static_cast<int>(NumberOfBins)) continue;
396  double factor = (1.0*(k-start))/(mNumberOfKeys-start);
397  double intensity = spectrum[m] * factor;
398  for (int n=2; n<=6; ++n)
399  {
400  double fn = getInharmonicPartial(n,f,B);
401  if (fn<20 or fn>10000) continue;
402  int mn = MathTools::roundToInteger(ftom(fn));
403  for (int i=mn-10; i<=mn+10; ++i)
404  if (i>=0 and i < static_cast<int>(NumberOfBins))
405  spectrum[i] = intensity * pow(4,-n) * exp(-0.1*(i-mn)*(i-mn));
406  }
407  }
408 }
409 
410 
411 
412 //-----------------------------------------------------------------------------
413 // Mollify
414 //-----------------------------------------------------------------------------
415 
417 {
418  SpectrumType &spectrum = key.getSpectrum();
419  SpectrumType copy = spectrum;
420  //double f1 = key.getRecordedFrequency();
421 
422 
423  //SpectrumType copy = spectrum;
424  int M = static_cast<int>(NumberOfBins);
425  for (int m=0; m<M; ++m)
426  {
427  double f = mtof(m);
428  const double df=55.0/f+f/2000.0;
429  int dm = MathTools::roundToInteger(ftom(f+df)) - m;
430  double sum=0,norm=0;
431  for (int ms = std::max(1,m-3*dm); ms<=std::min(m+3*dm,M-1); ms++)
432  {
433  double weight = exp(-1.0*(ms-m)*(ms-m)/dm/dm);
434  norm+=weight;
435  sum+=copy[ms]*weight;
436  }
437  if (norm>0) spectrum[m]=sum/norm;
438  }
439 
440 }
441 
442 
443 
444 
445 //-----------------------------------------------------------------------------
446 } // namespace entropyminimizer
double getRecordedFrequency() const
Get recorded frequency.
Definition: key.cpp:119
void extrapolateInharmonicity()
Extrapolate inharmonicity and estimate missing values.
#define LogW(...)
Definition: log.h:56
Class describing a single piano key.
Definition: key.h:45
static const int NumberOfBins
Total number of slots: 9 octaves.
Definition: key.h:50
void initializeSPLAFilter()
Initialize the filter function in the vector mdBA.
int roundToInteger(T x)
Round a floating point number to an integer.
Definition: mathtools.h:43
Definition: piano.h:40
double getEqualTempFrequency(int keynumber, double cents=0, double A4=0) const
Function returning the equal temperament.
Definition: piano.cpp:123
double getInharmonicPartial(double n, double f1, double B)
Compute frequency of inharmonic partial.
const double PI
Definition: mathtools.h:37
AuditoryPreprocessing(Piano &piano)
Constructor.
double getMeasuredInharmonicity() const
Get estimated inharmonicity.
Definition: key.cpp:141
void convertToSPLA(SpectrumType &s)
Compute A-weighted sound pressure level SPLA.
Namespace for all entropy minimizer components.
void normalize(std::vector< double > &vec)
Normalize a distribution stored in a vector.
Definition: mathtools.cpp:92
bool checkDataConsistency()
Check consistency of the piano dataset.
std::vector< double > SpectrumType
Type of a log-binned spectrum.
Definition: key.h:54
const SpectrumType & getSpectrum() const
Get a read-only reference to mSpectrum.
Definition: key.cpp:237
double getInharmonicPartialIndex(double f, double f1, double B)
#define EptAssert(a, b)
Definition: eptexception.h:47
void improveHighFrequencyPeaks()
Improve high-frequency spectral lines.
void normalizeSpectrum(Key &key)
Normalize spectrum.
void cutLowFrequencies(Key &key)
Cut all frequencies below 5/6*f1 in order to reduce noise.
void cleanSpectrum(Key &key)
Clean logarithmically binned spectrum.
double getExpectedInharmonicity(double f) const
Compute expected approximative inharmonicity.
Definition: piano.cpp:102
double computeNorm(std::vector< double > &vec)
Compute the norm of a vector.
Definition: mathtools.cpp:82