chris@52
|
1 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */
|
chris@52
|
2
|
chris@52
|
3 /*
|
chris@52
|
4 pYIN - A fundamental frequency estimator for monophonic audio
|
chris@52
|
5 Centre for Digital Music, Queen Mary, University of London.
|
chris@52
|
6
|
chris@52
|
7 This program is free software; you can redistribute it and/or
|
chris@52
|
8 modify it under the terms of the GNU General Public License as
|
chris@52
|
9 published by the Free Software Foundation; either version 2 of the
|
chris@52
|
10 License, or (at your option) any later version. See the file
|
chris@52
|
11 COPYING included with this distribution for more information.
|
chris@52
|
12 */
|
chris@52
|
13
|
chris@52
|
14 #include "Yin.h"
|
chris@52
|
15
|
chris@52
|
16 #include "vamp-sdk/FFT.h"
|
chris@52
|
17 #include "MeanFilter.h"
|
chris@52
|
18 #include "YinUtil.h"
|
chris@52
|
19
|
chris@52
|
20 #include <vector>
|
chris@52
|
21 #include <cstdlib>
|
chris@52
|
22 #include <cstdio>
|
chris@52
|
23 #include <cmath>
|
chris@52
|
24 #include <complex>
|
chris@52
|
25
|
chris@52
|
26 using std::vector;
|
chris@52
|
27
|
chris@52
|
28 Yin::Yin(size_t frameSize, size_t inputSampleRate, double thresh) :
|
chris@52
|
29 m_frameSize(frameSize),
|
chris@52
|
30 m_inputSampleRate(inputSampleRate),
|
chris@52
|
31 m_thresh(thresh),
|
chris@52
|
32 m_threshDistr(2),
|
chris@52
|
33 m_yinBufferSize(frameSize/2)
|
chris@52
|
34 {
|
chris@52
|
35 if (frameSize & (frameSize-1)) {
|
chris@52
|
36 // throw "N must be a power of two";
|
chris@52
|
37 }
|
chris@52
|
38 }
|
chris@52
|
39
|
chris@52
|
40 Yin::~Yin()
|
chris@52
|
41 {
|
chris@52
|
42 }
|
chris@52
|
43
|
chris@52
|
44 Yin::YinOutput
|
chris@52
|
45 Yin::process(const double *in) const {
|
chris@52
|
46
|
chris@52
|
47 double* yinBuffer = new double[m_yinBufferSize];
|
chris@52
|
48
|
chris@52
|
49 // calculate aperiodicity function for all periods
|
chris@52
|
50 YinUtil::fastDifference(in, yinBuffer, m_yinBufferSize);
|
chris@52
|
51 YinUtil::cumulativeDifference(yinBuffer, m_yinBufferSize);
|
chris@52
|
52
|
chris@52
|
53 int tau = 0;
|
chris@52
|
54 tau = YinUtil::absoluteThreshold(yinBuffer, m_yinBufferSize, m_thresh);
|
chris@52
|
55
|
chris@52
|
56 double interpolatedTau;
|
chris@52
|
57 double aperiodicity;
|
chris@52
|
58 double f0;
|
chris@52
|
59
|
chris@54
|
60 if (tau!=0 && tau!=m_yinBufferSize-1)
|
chris@52
|
61 {
|
chris@52
|
62 interpolatedTau = YinUtil::parabolicInterpolation(yinBuffer, abs(tau), m_yinBufferSize);
|
chris@52
|
63 f0 = m_inputSampleRate * (1.0 / interpolatedTau);
|
chris@52
|
64 } else {
|
chris@52
|
65 interpolatedTau = 0;
|
chris@52
|
66 f0 = 0;
|
chris@52
|
67 }
|
chris@52
|
68 double rms = std::sqrt(YinUtil::sumSquare(in, 0, m_yinBufferSize)/m_yinBufferSize);
|
chris@52
|
69 aperiodicity = yinBuffer[abs(tau)];
|
chris@52
|
70 // std::cerr << aperiodicity << std::endl;
|
chris@52
|
71 if (tau < 0) f0 = -f0;
|
chris@52
|
72
|
chris@52
|
73 Yin::YinOutput yo(f0, 1-aperiodicity, rms);
|
chris@52
|
74 for (size_t iBuf = 0; iBuf < m_yinBufferSize; ++iBuf)
|
chris@52
|
75 {
|
chris@52
|
76 yo.salience.push_back(yinBuffer[iBuf] < 1 ? 1-yinBuffer[iBuf] : 0); // why are the values sometimes < 0 if I don't check?
|
chris@52
|
77 }
|
chris@52
|
78
|
chris@52
|
79 delete [] yinBuffer;
|
chris@52
|
80 return yo;
|
chris@52
|
81 }
|
chris@52
|
82
|
chris@52
|
83 Yin::YinOutput
|
chris@52
|
84 Yin::processProbabilisticYin(const double *in) const {
|
chris@54
|
85
|
chris@52
|
86 double* yinBuffer = new double[m_yinBufferSize];
|
chris@52
|
87
|
chris@52
|
88 // calculate aperiodicity function for all periods
|
chris@52
|
89 YinUtil::fastDifference(in, yinBuffer, m_yinBufferSize);
|
chris@52
|
90 YinUtil::cumulativeDifference(yinBuffer, m_yinBufferSize);
|
chris@52
|
91
|
chris@52
|
92 vector<double> peakProbability = YinUtil::yinProb(yinBuffer, m_threshDistr, m_yinBufferSize);
|
chris@54
|
93
|
chris@54
|
94 // basic yin output
|
chris@54
|
95 Yin::YinOutput yo(0,0,0);
|
chris@54
|
96 for (size_t iBuf = 1; iBuf < m_yinBufferSize-1; ++iBuf)
|
chris@52
|
97 {
|
chris@52
|
98 if (peakProbability[iBuf] > 0)
|
chris@52
|
99 {
|
chris@52
|
100 double currentF0 =
|
chris@52
|
101 m_inputSampleRate * (1.0 /
|
chris@52
|
102 YinUtil::parabolicInterpolation(yinBuffer, iBuf, m_yinBufferSize));
|
chris@52
|
103 yo.freqProb.push_back(pair<double, double>(currentF0, peakProbability[iBuf]));
|
chris@52
|
104 }
|
chris@52
|
105 }
|
chris@52
|
106
|
chris@54
|
107 // add salience
|
chris@54
|
108 for (size_t iBuf = 0; iBuf < m_yinBufferSize; ++iBuf) {
|
chris@54
|
109 yo.salience.push_back(peakProbability[iBuf]);
|
chris@54
|
110 }
|
chris@54
|
111
|
chris@52
|
112 // std::cerr << yo.freqProb.size() << std::endl;
|
chris@52
|
113
|
chris@52
|
114 delete [] yinBuffer;
|
chris@52
|
115 return yo;
|
chris@52
|
116 }
|
chris@52
|
117
|
chris@52
|
118
|
chris@52
|
119 int
|
chris@52
|
120 Yin::setThreshold(double parameter)
|
chris@52
|
121 {
|
chris@52
|
122 m_thresh = static_cast<float>(parameter);
|
chris@52
|
123 return 0;
|
chris@52
|
124 }
|
chris@52
|
125
|
chris@52
|
126 int
|
chris@52
|
127 Yin::setThresholdDistr(float parameter)
|
chris@52
|
128 {
|
chris@52
|
129 m_threshDistr = static_cast<size_t>(parameter);
|
chris@52
|
130 return 0;
|
chris@52
|
131 }
|
chris@52
|
132
|
chris@52
|
133 int
|
chris@52
|
134 Yin::setFrameSize(size_t parameter)
|
chris@52
|
135 {
|
chris@52
|
136 m_frameSize = parameter;
|
chris@52
|
137 m_yinBufferSize = m_frameSize/2;
|
chris@52
|
138 return 0;
|
chris@52
|
139 }
|
chris@52
|
140
|
chris@52
|
141 // int
|
chris@52
|
142 // Yin::setRemoveUnvoiced(bool parameter)
|
chris@52
|
143 // {
|
chris@52
|
144 // m_removeUnvoiced = parameter;
|
chris@52
|
145 // return 0;
|
chris@52
|
146 // }
|
chris@54
|
147
|
chris@54
|
148 float
|
chris@54
|
149 Yin::constrainedMinPick(const double *in, const float minFreq, const int maxFreq) const {
|
chris@54
|
150
|
chris@54
|
151 double* yinBuffer = new double[m_yinBufferSize];
|
chris@54
|
152
|
chris@54
|
153 // calculate aperiodicity function for all periods
|
chris@54
|
154 YinUtil::fastDifference(in, yinBuffer, m_yinBufferSize);
|
chris@54
|
155 YinUtil::cumulativeDifference(yinBuffer, m_yinBufferSize);
|
chris@54
|
156
|
chris@54
|
157 int minPeriod = m_inputSampleRate / maxFreq;
|
chris@54
|
158 int maxPeriod = m_inputSampleRate / minFreq;
|
chris@54
|
159
|
chris@54
|
160 if (minPeriod < 0 || maxPeriod > m_yinBufferSize || minPeriod > maxPeriod) {
|
chris@54
|
161 delete [] yinBuffer;
|
chris@54
|
162 return 0.f;
|
chris@54
|
163 }
|
chris@54
|
164
|
chris@54
|
165 float bestVal = 1000;
|
chris@54
|
166 int bestTau = 0;
|
chris@54
|
167 for (int tau = minPeriod; tau <= maxPeriod; ++tau)
|
chris@54
|
168 {
|
chris@54
|
169 if (yinBuffer[tau] < bestVal)
|
chris@54
|
170 {
|
chris@54
|
171 bestVal = yinBuffer[tau];
|
chris@54
|
172 bestTau = tau;
|
chris@54
|
173 }
|
chris@54
|
174 }
|
chris@54
|
175
|
chris@54
|
176 float interpolatedTau =
|
chris@54
|
177 YinUtil::parabolicInterpolation(yinBuffer, bestTau, m_yinBufferSize);
|
chris@54
|
178
|
chris@54
|
179 delete [] yinBuffer;
|
chris@54
|
180 return m_inputSampleRate * (1.0 / interpolatedTau);
|
matthiasm@49
|
181 } |