comparison dsp/rateconversion/Resampler.cpp @ 366:767947956fc1

More resampler fixes (particularly to latency calculation) and tests
author Chris Cannam <c.cannam@qmul.ac.uk>
date Mon, 14 Oct 2013 16:15:32 +0100
parents b73dad5e6201
children f8fc21365a8c
comparison
equal deleted inserted replaced
365:b73dad5e6201 366:767947956fc1
8 8
9 #include <iostream> 9 #include <iostream>
10 #include <vector> 10 #include <vector>
11 11
12 using std::vector; 12 using std::vector;
13
14 //#define DEBUG_RESAMPLER 1
13 15
14 Resampler::Resampler(int sourceRate, int targetRate) : 16 Resampler::Resampler(int sourceRate, int targetRate) :
15 m_sourceRate(sourceRate), 17 m_sourceRate(sourceRate),
16 m_targetRate(targetRate) 18 m_targetRate(targetRate)
17 { 19 {
50 kw.cut(filter); 52 kw.cut(filter);
51 53
52 int inputSpacing = m_targetRate / m_gcd; 54 int inputSpacing = m_targetRate / m_gcd;
53 int outputSpacing = m_sourceRate / m_gcd; 55 int outputSpacing = m_sourceRate / m_gcd;
54 56
55 m_latency = int(ceil((m_filterLength / 2.0) / outputSpacing)); 57 #ifdef DEBUG_RESAMPLER
56 58 std::cerr << "resample " << m_sourceRate << " -> " << m_targetRate
57 int bufferLength = 0; 59 << ": inputSpacing " << inputSpacing << ", outputSpacing "
60 << outputSpacing << ": filter length " << m_filterLength
61 << std::endl;
62 #endif
58 63
59 m_phaseData = new Phase[inputSpacing]; 64 m_phaseData = new Phase[inputSpacing];
60 65
61 for (int phase = 0; phase < inputSpacing; ++phase) { 66 for (int phase = 0; phase < inputSpacing; ++phase) {
62 67
64 69
65 p.nextPhase = phase - outputSpacing; 70 p.nextPhase = phase - outputSpacing;
66 while (p.nextPhase < 0) p.nextPhase += inputSpacing; 71 while (p.nextPhase < 0) p.nextPhase += inputSpacing;
67 p.nextPhase %= inputSpacing; 72 p.nextPhase %= inputSpacing;
68 73
69 p.drop = int(ceil(std::max(0, outputSpacing - phase) / inputSpacing)); 74 p.drop = int(ceil(std::max(0.0, double(outputSpacing - phase))
70 p.take = int((outputSpacing + 75 / inputSpacing));
71 ((m_filterLength - 1 - phase) % inputSpacing)) 76
72 / outputSpacing); 77 int filtZipLength = int(ceil(double(m_filterLength - phase)
73 78 / inputSpacing));
74 int filtZipLength = int(ceil((m_filterLength - phase) / inputSpacing));
75 if (filtZipLength > bufferLength) {
76 bufferLength = filtZipLength;
77 }
78
79 for (int i = 0; i < filtZipLength; ++i) { 79 for (int i = 0; i < filtZipLength; ++i) {
80 p.filter.push_back(filter[i * inputSpacing + phase]); 80 p.filter.push_back(filter[i * inputSpacing + phase]);
81 } 81 }
82 82
83 m_phaseData[phase] = p; 83 m_phaseData[phase] = p;
84 } 84 }
85
86 #ifdef DEBUG_RESAMPLER
87 for (int phase = 0; phase < inputSpacing; ++phase) {
88 std::cerr << "filter for phase " << phase << " of " << inputSpacing << " (with length " << m_phaseData[phase].filter.size() << "):";
89 for (int i = 0; i < m_phaseData[phase].filter.size(); ++i) {
90 if (i % 4 == 0) {
91 std::cerr << std::endl << i << ": ";
92 }
93 float v = m_phaseData[phase].filter[i];
94 if (v == 1) {
95 std::cerr << " *** " << v << " *** ";
96 } else {
97 std::cerr << v << " ";
98 }
99 }
100 std::cerr << std::endl;
101 }
102 #endif
85 103
86 delete[] filter; 104 delete[] filter;
87 105
88 // The May implementation of this uses a pull model -- we ask the 106 // The May implementation of this uses a pull model -- we ask the
89 // resampler for a certain number of output samples, and it asks 107 // resampler for a certain number of output samples, and it asks
100 // either return a very variable number of samples (none at all 118 // either return a very variable number of samples (none at all
101 // until the filter fills, then half the filter length at once) or 119 // until the filter fills, then half the filter length at once) or
102 // else have a lengthy declared latency on the output. We do the 120 // else have a lengthy declared latency on the output. We do the
103 // latter. (What do other implementations do?) 121 // latter. (What do other implementations do?)
104 122
105 m_phase = m_filterLength % inputSpacing; 123 m_phase = (m_filterLength/2) % inputSpacing;
106 m_buffer = vector<double>(bufferLength, 0); 124
125 m_buffer = vector<double>(m_phaseData[0].filter.size(), 0);
126
127 m_latency =
128 ((m_buffer.size() * inputSpacing) - (m_filterLength/2)) / outputSpacing
129 + m_phase;
130
131 #ifdef DEBUG_RESAMPLER
132 std::cerr << "initial phase " << m_phase << " (as " << (m_filterLength/2) << " % " << inputSpacing << ")"
133 << ", latency " << m_latency << std::endl;
134 #endif
107 } 135 }
108 136
109 double 137 double
110 Resampler::reconstructOne(const double *src) 138 Resampler::reconstructOne()
111 { 139 {
112 Phase &pd = m_phaseData[m_phase]; 140 Phase &pd = m_phaseData[m_phase];
113 double *filt = pd.filter.data(); 141 double *filt = pd.filter.data();
142 double v = 0.0;
114 int n = pd.filter.size(); 143 int n = pd.filter.size();
115 double v = 0.0;
116 for (int i = 0; i < n; ++i) { 144 for (int i = 0; i < n; ++i) {
117 v += m_buffer[i] * filt[i]; 145 v += m_buffer[i] * filt[i];
118 } 146 }
119 m_buffer = vector<double>(m_buffer.begin() + pd.drop, m_buffer.end()); 147 m_buffer = vector<double>(m_buffer.begin() + pd.drop, m_buffer.end());
120 for (int i = 0; i < pd.take; ++i) { 148 m_phase = pd.nextPhase;
149 return v;
150 }
151
152 int
153 Resampler::process(const double *src, double *dst, int n)
154 {
155 for (int i = 0; i < n; ++i) {
121 m_buffer.push_back(src[i]); 156 m_buffer.push_back(src[i]);
122 } 157 }
123 return v; 158
124 } 159 int maxout = int(ceil(double(n) * m_targetRate / m_sourceRate));
125 160 int outidx = 0;
126 int 161
127 Resampler::process(const double *src, double *dst, int remaining) 162 #ifdef DEBUG_RESAMPLER
128 { 163 std::cerr << "process: buf siz " << m_buffer.size() << " filt siz for phase " << m_phase << " " << m_phaseData[m_phase].filter.size() << std::endl;
129 int m = 0; 164 #endif
130 int offset = 0; 165
131 166 while (outidx < maxout &&
132 while (remaining >= m_phaseData[m_phase].take) { 167 m_buffer.size() >= m_phaseData[m_phase].filter.size()) {
133 // std::cerr << "remaining = " << remaining << ", m = " << m << ", take = " << m_phaseData[m_phase].take << std::endl; 168 dst[outidx] = reconstructOne();
134 int advance = m_phaseData[m_phase].take; 169 outidx++;
135 dst[m] = reconstructOne(src + offset); 170 }
136 offset += advance; 171
137 remaining -= advance; 172 return outidx;
138 m_phase = m_phaseData[m_phase].nextPhase; 173 }
139 // std::cerr << "remaining -> " << remaining << ", new phase has advance " << m_phaseData[m_phase].take << std::endl; 174
140 ++m;
141 }
142
143 // if (remaining > 0) {
144 // std::cerr << "have " << remaining << " spare, pushing to buffer" << std::endl;
145 // }
146
147 for (int i = 0; i < remaining; ++i) {
148 m_buffer.push_back(src[offset + i]);
149 }
150
151 return m;
152 }
153
154 std::vector<double> 175 std::vector<double>
155 Resampler::resample(int sourceRate, int targetRate, const double *data, int n) 176 Resampler::resample(int sourceRate, int targetRate, const double *data, int n)
156 { 177 {
157 Resampler r(sourceRate, targetRate); 178 Resampler r(sourceRate, targetRate);
158 179
159 int latency = r.getLatency(); 180 int latency = r.getLatency();
160 181
161 int m = int(ceil((n * targetRate) / sourceRate)); 182 int m = int(ceil(double(n * targetRate) / sourceRate));
162 int m1 = m + latency; 183 int m1 = m + latency;
163 int n1 = int((m1 * sourceRate) / targetRate); 184 int n1 = int(double(m1 * sourceRate) / targetRate);
164 185
165 vector<double> pad(n1 - n, 0.0); 186 vector<double> pad(n1 - n, 0.0);
166 vector<double> out(m1, 0.0); 187 vector<double> out(m1, 0.0);
167 188
168 int got = r.process(data, out.data(), n); 189 int got = r.process(data, out.data(), n);
169 got += r.process(pad.data(), out.data() + got, pad.size()); 190 got += r.process(pad.data(), out.data() + got, pad.size());
170 191
192 #ifdef DEBUG_RESAMPLER
193 std::cerr << "resample: " << n << " in, " << got << " out" << std::endl;
194 for (int i = 0; i < got; ++i) {
195 if (i % 5 == 0) std::cout << std::endl << i << "... ";
196 std::cout << (float) out[i] << " ";
197 }
198 std::cout << std::endl;
199 #endif
200
171 return vector<double>(out.begin() + latency, out.begin() + got); 201 return vector<double>(out.begin() + latency, out.begin() + got);
172 } 202 }
173 203