Line data Source code
1 : /*
2 : * SPDX-License-Identifier: GPL-3.0-only
3 : * MuseScore-CLA-applies
4 : *
5 : * MuseScore Studio
6 : * Music Composition & Notation
7 : *
8 : * Copyright (C) 2021 MuseScore Limited and others
9 : *
10 : * This program is free software: you can redistribute it and/or modify
11 : * it under the terms of the GNU General Public License version 3 as
12 : * published by the Free Software Foundation.
13 : *
14 : * This program is distributed in the hope that it will be useful,
15 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : * GNU General Public License for more details.
18 : *
19 : * You should have received a copy of the GNU General Public License
20 : * along with this program. If not, see <https://www.gnu.org/licenses/>.
21 : */
22 : #include "samplerateconvertor.h"
23 :
24 : #include <cmath>
25 :
26 : #include "log.h"
27 :
28 : using namespace muse::audio::engine;
29 :
30 0 : SampleRateConvertor::SampleRateConvertor(const std::vector<float>& data,
31 : unsigned int channelsCount,
32 : unsigned int sampleRateIn,
33 0 : unsigned int sampleRateOut)
34 0 : : m_data(data), m_fir(), m_channelsCount(channelsCount), m_sampleRateIn(sampleRateIn), m_sampleRateOut(sampleRateOut)
35 : {
36 0 : m_fir.resize(FIR_LENGTH, 0);
37 0 : m_y.resize(FIR_LENGTH, 0);
38 0 : initWindow();
39 0 : }
40 :
41 0 : std::vector<float> SampleRateConvertor::convert()
42 : {
43 0 : std::vector<float> out;
44 0 : auto resultSamples = m_data.size() * m_sampleRateOut / (m_channelsCount * m_sampleRateIn);
45 :
46 0 : out.resize(resultSamples * m_channelsCount);
47 0 : for (unsigned sample = 0; sample < resultSamples; ++sample) {
48 0 : for (unsigned int channel = 0; channel < m_channelsCount; ++channel) {
49 0 : out[sample * m_channelsCount + channel] = y(sample, channel);
50 : }
51 : }
52 :
53 0 : return out;
54 0 : }
55 :
56 0 : unsigned int SampleRateConvertor::convert(float* buffer, unsigned int from, unsigned int count)
57 : {
58 0 : unsigned int sample = from;
59 0 : unsigned int converted = 0;
60 :
61 0 : for (converted = 0; converted < count; ++converted, ++sample) {
62 0 : for (unsigned int channel = 0; channel < m_channelsCount; ++channel) {
63 0 : if (!availableSamples(sample)) {
64 0 : return converted;
65 : }
66 0 : buffer[converted * m_channelsCount + channel] = y(sample, channel);
67 : }
68 : }
69 0 : return converted;
70 : }
71 :
72 0 : float SampleRateConvertor::y(unsigned int sample, unsigned int channel) const
73 : {
74 0 : switch (m_method) {
75 0 : case SINC:
76 0 : return ySinc(sample, channel);
77 0 : case FIR:
78 0 : return yFIR(sample, channel);
79 : }
80 0 : return 0.f;
81 : }
82 :
83 0 : void SampleRateConvertor::setChannelCount(unsigned int count)
84 : {
85 0 : m_channelsCount = count;
86 0 : }
87 :
88 0 : void SampleRateConvertor::setSampleRateIn(unsigned int sampleRate)
89 : {
90 0 : if (m_sampleRateIn != sampleRate) {
91 0 : m_sampleRateIn = sampleRate;
92 0 : initWindow();
93 : }
94 0 : }
95 :
96 0 : void SampleRateConvertor::setSampleRateOut(unsigned int sampleRate)
97 : {
98 0 : if (m_sampleRateOut != sampleRate) {
99 0 : m_sampleRateOut = sampleRate;
100 0 : initWindow();
101 : }
102 0 : }
103 :
104 0 : float SampleRateConvertor::ySinc(unsigned int sample, unsigned int channel) const
105 : {
106 0 : double currentOutputSampleTime = sample / static_cast<double>(m_sampleRateOut);
107 0 : auto zeroInputSample = std::floor(currentOutputSampleTime * m_sampleRateIn);
108 :
109 0 : float out = 0.f;
110 0 : for (unsigned int i = 0; i < USE_SAMPLES; ++i) {
111 0 : int currentSample = zeroInputSample + i - USE_SAMPLES / 2;
112 0 : double currentInputSampleTime = currentSample / static_cast<double>(m_sampleRateIn);
113 :
114 0 : if (currentSample < 0 || currentSample * m_channelsCount + channel >= m_data.size()) {
115 0 : continue;
116 : }
117 0 : out += m_data.at(currentSample * m_channelsCount + channel) * sinc(currentOutputSampleTime - currentInputSampleTime);
118 : }
119 0 : return out;
120 : }
121 :
122 0 : float SampleRateConvertor::yFIR(unsigned int sample, unsigned int channel) const
123 : {
124 0 : auto x = [this, &channel](int m) -> float {
125 0 : int pos = m * m_M / m_L * m_channelsCount + channel;
126 0 : if (pos >= 0 && pos < static_cast<int>(m_data.size())) {
127 0 : return m_data.at(pos);
128 : }
129 0 : return 0.f;
130 0 : };
131 :
132 0 : float y = 0.f;
133 0 : for (unsigned int i = 0; i < FIR_LENGTH; ++i) {
134 0 : int t = sample - static_cast<int>(i);
135 0 : y += x(t) * m_fir[i];
136 : }
137 0 : return y;
138 : }
139 :
140 0 : bool SampleRateConvertor::availableSamples(unsigned int sample) const
141 : {
142 0 : float currentOutputSampleTime = sample / static_cast<float>(m_sampleRateOut);
143 0 : auto firstInputSample = std::floor(currentOutputSampleTime * m_sampleRateIn) - USE_SAMPLES / 2;
144 :
145 : // first sample for conversion points out of the input buffer
146 0 : if (firstInputSample * m_channelsCount >= m_data.size()) {
147 0 : return false;
148 : }
149 0 : return true;
150 : }
151 :
152 0 : inline double zeroBessel(double x)
153 : {
154 : //return std::cyl_bessel_i(0, x);
155 :
156 0 : double s = 1, y = 1;
157 0 : int m = 1;
158 :
159 0 : while (y > std::numeric_limits<float>::min()) {
160 0 : m += 2;
161 0 : y *= x * x / (4 * m * m);
162 0 : s += s * y;
163 : }
164 :
165 0 : return s;
166 : }
167 :
168 0 : void SampleRateConvertor::initWindow()
169 : {
170 0 : int max = std::max(m_sampleRateIn, m_sampleRateOut), min = std::min(m_sampleRateIn, m_sampleRateOut);
171 0 : int maxCommonDivider = 1;
172 0 : while (max % min != 0) {
173 0 : maxCommonDivider = max % min;
174 0 : max = min;
175 0 : min = maxCommonDivider;
176 : }
177 0 : m_M = m_sampleRateIn / maxCommonDivider;
178 0 : m_L = m_sampleRateOut / maxCommonDivider;
179 0 : double fStop = std::min(m_sampleRateIn, m_sampleRateOut) / 2,
180 0 : fIntermediateSampleRate = m_sampleRateIn * m_M,
181 0 : attenuation = 96 /*dB*/;
182 0 : int M = FIR_LENGTH;
183 0 : int Np = (M - 1) / 2;
184 :
185 0 : double alpha = 0.1102 * (attenuation - 8.7);
186 : double A[FIR_LENGTH];
187 :
188 0 : A[0] = 2 * fStop / fIntermediateSampleRate;
189 0 : for (int j = 1; j <= Np; j++) {
190 0 : A[j] = std::sin(2 * j * M_PI * fStop / fIntermediateSampleRate) / j * M_PI;
191 : }
192 :
193 0 : for (int j = 0; j <= Np; j++) {
194 0 : m_fir[Np + j] = A[j];
195 0 : m_fir[Np + j] *= zeroBessel(alpha * std::sqrt(1 - (j * j / (Np * Np))));
196 0 : m_fir[Np + j] /= zeroBessel(alpha);
197 : }
198 :
199 0 : for (int j = 0; j < Np; j++) {
200 0 : m_fir[j] = m_fir[M - 1 - j];
201 : }
202 0 : }
203 :
204 0 : double SampleRateConvertor::sinc(double value) const
205 : {
206 0 : if (value == 0) {
207 0 : return 1.f;
208 : }
209 0 : double sincArg = M_PI * value * m_sampleRateIn;
210 0 : double gain = std::min(1.0, m_sampleRateOut / static_cast<double>(m_sampleRateIn));
211 0 : return gain * std::sin(sincArg) / sincArg;
212 : }
|
