1 // Licensed GNU LGPL v3 or later: http://www.gnu.org/licenses/lgpl.html
2
3 #include "smsinedecoder.hh"
4 #include "smmath.hh"
5 #include "smfft.hh"
6 #include "smifftsynth.hh"
7 #include "smaudio.hh"
8 #include "smutils.hh"
9 #include "smalignedarray.hh"
10 #include <assert.h>
11 #include <stdio.h>
12 #include <math.h>
13
14 using SpectMorph::SineDecoder;
15 using SpectMorph::AudioBlock;
16 using std::vector;
17
18 /**
19 * \brief Constructor setting up the various decoding parameters
20 *
21 * @param mix_freq sample rate to be used for reconstruction
22 * @param frame_size frame size (in samples)
23 * @param frame_step frame step (in samples)
24 * @param mode selects decoding algorithm to be used
25 */
SineDecoder(double fundamental_freq,double mix_freq,size_t frame_size,size_t frame_step,Mode mode)26 SineDecoder::SineDecoder (double fundamental_freq, double mix_freq, size_t frame_size, size_t frame_step, Mode mode) :
27 mix_freq (mix_freq),
28 fundamental_freq (fundamental_freq),
29 frame_size (frame_size),
30 frame_step (frame_step),
31 mode (mode)
32 {
33 ifft_synth = NULL;
34 }
35
~SineDecoder()36 SineDecoder::~SineDecoder()
37 {
38 if (ifft_synth)
39 {
40 delete ifft_synth;
41 ifft_synth = NULL;
42 }
43 }
44
45 /**
46 * \brief Function which decodes a part of the signal.
47 *
48 * This needs two adjecant frames as arguments.
49 *
50 * @param frame the current frame (the frame to be decoded)
51 * @param next_frame the frame after the current frame
52 * @param window the reconstruction window used for MODE_PHASE_SYNC_OVERLAP
53 */
54 void
process(const AudioBlock & block,const AudioBlock & next_block,const vector<double> & window,vector<float> & decoded_sines)55 SineDecoder::process (const AudioBlock& block,
56 const AudioBlock& next_block,
57 const vector<double>& window,
58 vector<float>& decoded_sines)
59 {
60 /* phase synchronous reconstruction (no loops) */
61 if (mode == MODE_PHASE_SYNC_OVERLAP)
62 {
63 AlignedArray<float, 16> aligned_decoded_sines (frame_size);
64 for (size_t i = 0; i < block.freqs.size(); i++)
65 {
66 const double SA = double (frame_step) / double (frame_size) * 2.0;
67 const double mag_epsilon = 1e-8;
68
69 VectorSinParams params;
70 params.mag = block.mags_f (i) * SA;
71 if (params.mag > mag_epsilon)
72 {
73 params.mix_freq = mix_freq;
74 params.freq = block.freqs_f (i) * fundamental_freq;
75 params.phase = block.phases_f (i);
76 params.mode = VectorSinParams::ADD;
77
78 fast_vector_sinf (params, &aligned_decoded_sines[0], &aligned_decoded_sines[frame_size]);
79 }
80 }
81 for (size_t t = 0; t < frame_size; t++)
82 decoded_sines[t] = aligned_decoded_sines[t] * window[t];
83 return;
84 }
85 else if (mode == MODE_PHASE_SYNC_OVERLAP_IFFT)
86 {
87 const size_t block_size = frame_size;
88
89 if (!ifft_synth)
90 ifft_synth = new IFFTSynth (block_size, mix_freq, IFFTSynth::WIN_HANNING);
91
92 ifft_synth->clear_partials();
93 for (size_t i = 0; i < block.freqs.size(); i++)
94 {
95 const double SA = double (frame_step) / double (frame_size) * 2.0;
96 const double mag_epsilon = 1e-8;
97
98 const double mag = block.mags_f (i) * SA;
99 const double freq = block.freqs_f (i) * fundamental_freq;
100 if (mag > mag_epsilon)
101 ifft_synth->render_partial (freq, mag, block.phases_f (i));
102 }
103 ifft_synth->get_samples (&decoded_sines[0]);
104 return;
105 }
106
107 zero_float_block (decoded_sines.size(), &decoded_sines[0]);
108
109 /* phase distorted reconstruction */
110 vector<float> freqs (block.freqs.size());
111 vector<float> nfreqs (next_block.freqs.size());
112
113 for (size_t i = 0; i < freqs.size(); i++)
114 freqs[i] = block.freqs_f (i) * fundamental_freq;
115
116 for (size_t i = 0; i < nfreqs.size(); i++)
117 nfreqs[i] = next_block.freqs_f (i) * fundamental_freq;
118
119 int todo = freqs.size() + nfreqs.size();
120
121 synth_fixed_phase = next_synth_fixed_phase;
122 synth_fixed_phase.resize (freqs.size());
123 next_synth_fixed_phase.resize (nfreqs.size());
124
125 const double SIN_AMP = 1.0;
126 const bool TRACKING_SYNTH = true;
127 while (todo)
128 {
129 double best_delta = 1e10;
130 int best_i = 0, best_j = 0; /* init to get rid of gcc warning */
131 for (size_t i = 0; i < freqs.size(); i++)
132 {
133 for (size_t j = 0; j < nfreqs.size(); j++)
134 {
135 double delta = fabs (freqs[i] - nfreqs[j]) / freqs[i];
136 if (freqs[i] >= 0 && nfreqs[j] >= 0 && delta < best_delta && delta < 0.1)
137 {
138 best_delta = delta;
139 best_i = i;
140 best_j = j;
141 }
142 }
143 }
144 if (best_delta < 0.1)
145 {
146 double freq = freqs[best_i];
147 freqs[best_i] = -1;
148 double nfreq = nfreqs[best_j];
149 nfreqs[best_j] = -1;
150 double mag = block.mags_f (best_i);
151 double nmag = block.mags_f (best_j);
152
153 // fprintf (stderr, "%f | %f ==> %f | %f\n", freq, mag, nfreq, nmag);
154 assert (fabs (nfreq - freq) / freq < 0.1);
155
156 double phase_delta = 2 * M_PI * freq / mix_freq;
157 double nphase_delta = 2 * M_PI * nfreq / mix_freq;
158 double phase = synth_fixed_phase[best_i];
159 if (TRACKING_SYNTH)
160 {
161 for (size_t i = 0; i < frame_step; i++)
162 {
163 double inter = i / double (frame_step);
164
165 decoded_sines [i] += sin (phase) * ((1 - inter) * mag + inter * nmag) * SIN_AMP;
166 phase += (1 - inter) * phase_delta + inter * nphase_delta;
167 while (phase > 2 * M_PI)
168 phase -= 2 * M_PI;
169 }
170 next_synth_fixed_phase[best_j] = phase;
171 }
172 else
173 {
174 for (size_t i = 0; i < frame_size; i++)
175 {
176 decoded_sines [i] += sin (phase) * window[i] * mag * SIN_AMP;
177 phase += phase_delta;
178 while (phase > 2 * M_PI)
179 phase -= 2 * M_PI;
180 // nfreq phase required -> ramp
181 if (i == frame_step - 1)
182 next_synth_fixed_phase[best_j] = phase;
183 }
184 }
185 todo -= 2;
186 }
187 else
188 {
189 for (size_t from = 0; from < freqs.size(); from++)
190 {
191 if (freqs[from] > -1)
192 {
193 double freq = freqs[from];
194 freqs[from] = -1;
195 double mag = block.mags_f (from);
196
197 // fprintf (stderr, "%f | %f >>> \n", freq, mag);
198
199 double phase_delta = 2 * M_PI * freq / mix_freq;
200 double phase = synth_fixed_phase[from];
201 if (TRACKING_SYNTH)
202 {
203 for (size_t i = 0; i < frame_step; i++)
204 {
205 double inter = i / double (frame_step);
206
207 decoded_sines [i] += sin (phase) * (1 - inter) * mag * SIN_AMP;
208 phase += phase_delta;
209 while (phase > 2 * M_PI)
210 phase -= 2 * M_PI;
211 }
212 }
213 else
214 {
215 for (size_t i = 0; i < frame_size; i++)
216 {
217 decoded_sines [i] += sin (phase) * window[i] * mag * SIN_AMP;
218 phase += phase_delta;
219 while (phase > 2 * M_PI)
220 phase -= 2 * M_PI;
221 }
222 }
223 todo--;
224 }
225 }
226 for (size_t to = 0; to < nfreqs.size(); to++)
227 {
228 if (nfreqs[to] > -1)
229 {
230 double freq = nfreqs[to];
231 nfreqs[to] = -1;
232 double mag = next_block.mags_f (to);
233
234 // fprintf (stderr, "%f | %f <<< \n", freq, mag);
235
236 double phase_delta = 2 * M_PI * freq / mix_freq;
237 double phase = 0;
238 if (TRACKING_SYNTH)
239 {
240 for (size_t i = 0; i < frame_step; i++)
241 {
242 double inter = i / double (frame_step);
243
244 decoded_sines[i] += sin (phase) * inter * mag * SIN_AMP; /* XXX */
245 phase += phase_delta;
246 while (phase > 2 * M_PI)
247 phase -= 2 * M_PI;
248 }
249 next_synth_fixed_phase[to] = phase;
250 }
251 todo--;
252 }
253 }
254 }
255 }
256 }
257