]> git.sesse.net Git - nageru/blob - resampler.cpp
Update loop filter to take into account ratio; based on email from Fons Adriaensen.
[nageru] / resampler.cpp
1 // Parts of the code is adapted from Adriaensen's project Zita-ajbridge,
2 // although it has been heavily reworked for this use case. Original copyright follows:
3 //
4 //  Copyright (C) 2012-2015 Fons Adriaensen <fons@linuxaudio.org>
5 //    
6 //  This program is free software; you can redistribute it and/or modify
7 //  it under the terms of the GNU General Public License as published by
8 //  the Free Software Foundation; either version 3 of the License, or
9 //  (at your option) any later version.
10 //
11 //  This program is distributed in the hope that it will be useful,
12 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
13 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14 //  GNU General Public License for more details.
15 //
16 //  You should have received a copy of the GNU General Public License
17 //  along with this program.  If not, see <http://www.gnu.org/licenses/>.
18
19 #include "resampler.h"
20
21 #include <stdio.h>
22 #include <math.h>
23 #include <zita-resampler/vresampler.h>
24
25 Resampler::Resampler(unsigned freq_in, unsigned freq_out, unsigned num_channels)
26         : freq_in(freq_in), freq_out(freq_out), num_channels(num_channels),
27           ratio(double(freq_out) / double(freq_in))
28 {
29         vresampler.setup(ratio, num_channels, /*hlen=*/32);
30
31         // Prime the resampler so there's no more delay.
32         vresampler.inp_count = vresampler.inpsize() / 2 - 1;
33         vresampler.out_count = 1048576;
34         vresampler.process ();
35 }
36
37 void Resampler::add_input_samples(double pts, const float *samples, ssize_t num_samples)
38 {
39         if (first_input) {
40                 // Synthesize a fake length.
41                 last_input_len = double(num_samples) / freq_in;
42                 first_input = false;
43         } else {
44                 last_input_len = pts - last_input_pts;
45         }
46
47         last_input_pts = pts;
48
49         k_a0 = k_a1;
50         k_a1 += num_samples;
51
52         for (ssize_t i = 0; i < num_samples * num_channels; ++i) {
53                 buffer.push_back(samples[i]);
54         }
55 }
56
57 void Resampler::get_output_samples(double pts, float *samples, ssize_t num_samples)
58 {
59         double last_output_len;
60         if (first_output) {
61                 // Synthesize a fake length.
62                 last_output_len = double(num_samples) / freq_out;
63         } else {
64                 last_output_len = pts - last_output_pts;
65         }
66         last_output_pts = pts;
67
68         // Using the time point since just before the last call to add_input_samples() as a base,
69         // estimate actual delay based on activity since then, measured in number of input samples:
70         double actual_delay = 0.0;
71         actual_delay += (k_a1 - k_a0) * last_output_len / last_input_len;    // Inserted samples since k_a0, rescaled for the different time periods.
72         actual_delay += k_a0 - total_consumed_samples;                       // Samples inserted before k_a0 but not consumed yet.
73         actual_delay += vresampler.inpdist();                                // Delay in the resampler itself.
74         double err = actual_delay - expected_delay;
75         if (first_output && err < 0.0) {
76                 // Before the very first block, insert artificial delay based on our initial estimate,
77                 // so that we don't need a long period to stabilize at the beginning.
78                 int delay_samples_to_add = lrintf(-err);
79                 for (ssize_t i = 0; i < delay_samples_to_add * num_channels; ++i) {
80                         buffer.push_front(0.0f);
81                 }
82                 total_consumed_samples -= delay_samples_to_add;  // Equivalent to increasing k_a0 and k_a1.
83                 err += delay_samples_to_add;
84                 first_output = false;
85         }
86
87         // Compute loop filter coefficients for the two filters. We need to compute them
88         // every time, since they depend on the number of samples the user asked for.
89         //
90         // The loop bandwidth starts at 1.0 Hz, then goes down to 0.05 Hz after four seconds.
91         double loop_bandwidth_hz = (k_a0 < 4 * freq_in) ? 1.0 : 0.05;
92
93         // Set filters. The first filter much wider than the first one (20x as wide).
94         double w = (2.0 * M_PI) * loop_bandwidth_hz * num_samples / freq_out;
95         double w0 = 1.0 - exp(-20.0 * w);
96         double w1 = w * 1.5 / num_samples / ratio;
97         double w2 = w / 1.5;
98
99         // Filter <err> through the loop filter to find the correction ratio.
100         z1 += w0 * (w1 * err - z1);
101         z2 += w0 * (z1 - z2);
102         z3 += w2 * z2;
103         double rcorr = 1.0 - z2 - z3;
104         if (rcorr > 1.05) rcorr = 1.05;
105         if (rcorr < 0.95) rcorr = 0.95;
106         vresampler.set_rratio(rcorr);
107
108         // Finally actually resample, consuming exactly <num_samples> output samples.
109         vresampler.out_data = samples;
110         vresampler.out_count = num_samples;
111         while (vresampler.out_count > 0) {
112                 if (buffer.empty()) {
113                         // This should never happen unless delay is set way too low,
114                         // or we're dropping a lot of data.
115                         fprintf(stderr, "PANIC: Out of input samples to resample, still need %d output samples!\n",
116                                 int(vresampler.out_count));
117                         break;
118                 }
119
120                 float inbuf[1024];
121                 size_t num_input_samples = sizeof(inbuf) / (sizeof(float) * num_channels);
122                 if (num_input_samples * num_channels > buffer.size()) {
123                         num_input_samples = buffer.size() / num_channels;
124                 }
125                 for (size_t i = 0; i < num_input_samples * num_channels; ++i) {
126                         inbuf[i] = buffer[i];
127                 }
128
129                 vresampler.inp_count = num_input_samples;
130                 vresampler.inp_data = inbuf;
131
132                 vresampler.process();
133
134                 size_t consumed_samples = num_input_samples - vresampler.inp_count;
135                 total_consumed_samples += consumed_samples;
136                 buffer.erase(buffer.begin(), buffer.begin() + consumed_samples * num_channels);
137         }
138 }