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
|
matthiasm@70
|
28 Yin::Yin(size_t frameSize, size_t inputSampleRate, double thresh, bool fast) :
|
chris@52
|
29 m_frameSize(frameSize),
|
chris@52
|
30 m_inputSampleRate(inputSampleRate),
|
chris@52
|
31 m_thresh(thresh),
|
chris@52
|
32 m_threshDistr(2),
|
matthiasm@70
|
33 m_yinBufferSize(frameSize/2),
|
matthiasm@70
|
34 m_fast(fast)
|
chris@52
|
35 {
|
chris@52
|
36 if (frameSize & (frameSize-1)) {
|
chris@52
|
37 // throw "N must be a power of two";
|
chris@52
|
38 }
|
chris@52
|
39 }
|
chris@52
|
40
|
chris@52
|
41 Yin::~Yin()
|
chris@52
|
42 {
|
chris@52
|
43 }
|
chris@52
|
44
|
chris@52
|
45 Yin::YinOutput
|
chris@52
|
46 Yin::process(const double *in) const {
|
chris@52
|
47
|
chris@52
|
48 double* yinBuffer = new double[m_yinBufferSize];
|
chris@52
|
49 // calculate aperiodicity function for all periods
|
matthiasm@70
|
50 if (m_fast) YinUtil::fastDifference(in, yinBuffer, m_yinBufferSize);
|
matthiasm@70
|
51 else YinUtil::slowDifference(in, yinBuffer, m_yinBufferSize);
|
chris@52
|
52 YinUtil::cumulativeDifference(yinBuffer, m_yinBufferSize);
|
chris@52
|
53
|
chris@52
|
54 int tau = 0;
|
chris@52
|
55 tau = YinUtil::absoluteThreshold(yinBuffer, m_yinBufferSize, m_thresh);
|
chris@52
|
56
|
chris@52
|
57 double interpolatedTau;
|
chris@52
|
58 double aperiodicity;
|
chris@52
|
59 double f0;
|
chris@52
|
60
|
Chris@61
|
61 if (tau!=0 && tau!=int(m_yinBufferSize)-1)
|
chris@52
|
62 {
|
chris@52
|
63 interpolatedTau = YinUtil::parabolicInterpolation(yinBuffer, abs(tau), m_yinBufferSize);
|
chris@52
|
64 f0 = m_inputSampleRate * (1.0 / interpolatedTau);
|
chris@52
|
65 } else {
|
chris@52
|
66 interpolatedTau = 0;
|
chris@52
|
67 f0 = 0;
|
chris@52
|
68 }
|
chris@52
|
69 double rms = std::sqrt(YinUtil::sumSquare(in, 0, m_yinBufferSize)/m_yinBufferSize);
|
chris@52
|
70 aperiodicity = yinBuffer[abs(tau)];
|
chris@52
|
71 // std::cerr << aperiodicity << std::endl;
|
chris@52
|
72 if (tau < 0) f0 = -f0;
|
chris@52
|
73
|
chris@52
|
74 Yin::YinOutput yo(f0, 1-aperiodicity, rms);
|
chris@52
|
75 for (size_t iBuf = 0; iBuf < m_yinBufferSize; ++iBuf)
|
chris@52
|
76 {
|
chris@52
|
77 yo.salience.push_back(yinBuffer[iBuf] < 1 ? 1-yinBuffer[iBuf] : 0); // why are the values sometimes < 0 if I don't check?
|
chris@52
|
78 }
|
chris@52
|
79
|
chris@52
|
80 delete [] yinBuffer;
|
chris@52
|
81 return yo;
|
chris@52
|
82 }
|
chris@52
|
83
|
chris@52
|
84 Yin::YinOutput
|
chris@52
|
85 Yin::processProbabilisticYin(const double *in) const {
|
chris@54
|
86
|
chris@52
|
87 double* yinBuffer = new double[m_yinBufferSize];
|
chris@52
|
88
|
chris@52
|
89 // calculate aperiodicity function for all periods
|
matthiasm@70
|
90 if (m_fast) YinUtil::fastDifference(in, yinBuffer, m_yinBufferSize);
|
matthiasm@70
|
91 else YinUtil::slowDifference(in, yinBuffer, m_yinBufferSize);
|
matthiasm@70
|
92
|
chris@52
|
93 YinUtil::cumulativeDifference(yinBuffer, m_yinBufferSize);
|
chris@52
|
94
|
chris@52
|
95 vector<double> peakProbability = YinUtil::yinProb(yinBuffer, m_threshDistr, m_yinBufferSize);
|
chris@54
|
96
|
chris@54
|
97 // basic yin output
|
chris@54
|
98 Yin::YinOutput yo(0,0,0);
|
Chris@61
|
99 for (int iBuf = 1; iBuf < int(m_yinBufferSize)-1; ++iBuf)
|
chris@52
|
100 {
|
chris@52
|
101 if (peakProbability[iBuf] > 0)
|
chris@52
|
102 {
|
chris@52
|
103 double currentF0 =
|
chris@52
|
104 m_inputSampleRate * (1.0 /
|
chris@52
|
105 YinUtil::parabolicInterpolation(yinBuffer, iBuf, m_yinBufferSize));
|
chris@52
|
106 yo.freqProb.push_back(pair<double, double>(currentF0, peakProbability[iBuf]));
|
chris@52
|
107 }
|
chris@52
|
108 }
|
chris@52
|
109
|
chris@54
|
110 // add salience
|
chris@54
|
111 for (size_t iBuf = 0; iBuf < m_yinBufferSize; ++iBuf) {
|
chris@54
|
112 yo.salience.push_back(peakProbability[iBuf]);
|
chris@54
|
113 }
|
chris@54
|
114
|
chris@52
|
115 // std::cerr << yo.freqProb.size() << std::endl;
|
chris@52
|
116
|
chris@52
|
117 delete [] yinBuffer;
|
chris@52
|
118 return yo;
|
chris@52
|
119 }
|
chris@52
|
120
|
chris@52
|
121
|
chris@52
|
122 int
|
chris@52
|
123 Yin::setThreshold(double parameter)
|
chris@52
|
124 {
|
chris@52
|
125 m_thresh = static_cast<float>(parameter);
|
chris@52
|
126 return 0;
|
chris@52
|
127 }
|
chris@52
|
128
|
chris@52
|
129 int
|
chris@52
|
130 Yin::setThresholdDistr(float parameter)
|
chris@52
|
131 {
|
chris@52
|
132 m_threshDistr = static_cast<size_t>(parameter);
|
chris@52
|
133 return 0;
|
chris@52
|
134 }
|
chris@52
|
135
|
chris@52
|
136 int
|
chris@52
|
137 Yin::setFrameSize(size_t parameter)
|
chris@52
|
138 {
|
chris@52
|
139 m_frameSize = parameter;
|
chris@52
|
140 m_yinBufferSize = m_frameSize/2;
|
chris@52
|
141 return 0;
|
chris@52
|
142 }
|
chris@52
|
143
|
matthiasm@70
|
144 int
|
matthiasm@70
|
145 Yin::setFast(bool fast)
|
matthiasm@70
|
146 {
|
matthiasm@70
|
147 m_fast = fast;
|
matthiasm@70
|
148 return 0;
|
matthiasm@70
|
149 }
|
matthiasm@70
|
150
|
chris@52
|
151 // int
|
chris@52
|
152 // Yin::setRemoveUnvoiced(bool parameter)
|
chris@52
|
153 // {
|
chris@52
|
154 // m_removeUnvoiced = parameter;
|
chris@52
|
155 // return 0;
|
chris@52
|
156 // }
|
chris@54
|
157
|
chris@54
|
158 float
|
chris@54
|
159 Yin::constrainedMinPick(const double *in, const float minFreq, const int maxFreq) const {
|
chris@54
|
160
|
chris@54
|
161 double* yinBuffer = new double[m_yinBufferSize];
|
chris@54
|
162
|
chris@54
|
163 // calculate aperiodicity function for all periods
|
matthiasm@60
|
164 YinUtil::slowDifference(in, yinBuffer, m_yinBufferSize);
|
chris@54
|
165 YinUtil::cumulativeDifference(yinBuffer, m_yinBufferSize);
|
chris@54
|
166
|
chris@54
|
167 int minPeriod = m_inputSampleRate / maxFreq;
|
chris@54
|
168 int maxPeriod = m_inputSampleRate / minFreq;
|
chris@54
|
169
|
Chris@61
|
170 if (minPeriod < 0 || maxPeriod > int(m_yinBufferSize) || minPeriod > maxPeriod) {
|
chris@54
|
171 delete [] yinBuffer;
|
chris@54
|
172 return 0.f;
|
chris@54
|
173 }
|
chris@54
|
174
|
chris@54
|
175 float bestVal = 1000;
|
chris@54
|
176 int bestTau = 0;
|
chris@54
|
177 for (int tau = minPeriod; tau <= maxPeriod; ++tau)
|
chris@54
|
178 {
|
chris@54
|
179 if (yinBuffer[tau] < bestVal)
|
chris@54
|
180 {
|
chris@54
|
181 bestVal = yinBuffer[tau];
|
chris@54
|
182 bestTau = tau;
|
chris@54
|
183 }
|
chris@54
|
184 }
|
chris@54
|
185
|
chris@54
|
186 float interpolatedTau =
|
chris@54
|
187 YinUtil::parabolicInterpolation(yinBuffer, bestTau, m_yinBufferSize);
|
chris@54
|
188
|
chris@54
|
189 delete [] yinBuffer;
|
chris@54
|
190 return m_inputSampleRate * (1.0 / interpolatedTau);
|
Chris@61
|
191 }
|