1 // Licensed GNU LGPL v3 or later: http://www.gnu.org/licenses/lgpl.html
2
3 #include "smmain.hh"
4 #include "smrandom.hh"
5 #include "smencoder.hh"
6 #include "smlivedecoder.hh"
7 #include "smminiresampler.hh"
8 #include "smfft.hh"
9 #include "smdebug.hh"
10
11 #include <assert.h>
12 #include <unistd.h>
13
14 using namespace SpectMorph;
15
16 using std::vector;
17 using std::string;
18 using std::max;
19
20 /* this generates noise signals for any sample rate >= 48000 with noise in frequency range [0,24000] */
21 class NoiseGenerator
22 {
23 vector<float> noise_96000;
24 public:
25 vector<float> gen_noise (int sr);
26 };
27
28 vector<float>
gen_noise(int sr)29 NoiseGenerator::gen_noise (int sr)
30 {
31 assert (sr >= 48000);
32
33 Random random;
34 if (!noise_96000.size())
35 {
36 const size_t FFT_SIZE = 512 * 1024; // should be power of two, greater than 5 seconds at 96000 (480000)
37
38 float *in = FFT::new_array_float (FFT_SIZE);
39 float *out = FFT::new_array_float (FFT_SIZE);
40
41 const double MAG = sqrt (0.5);
42 for (size_t i = 0; i < FFT_SIZE; i++)
43 in[i] = random.random_double_range (-MAG, MAG);
44
45 vector<float> old_in (in, in+FFT_SIZE);
46
47 FFT::fftar_float (FFT_SIZE, in, out, FFT::PLAN_ESTIMATE);
48
49 // band limit noise spectrum, highest permitted frequency: 24000 Hz
50 for (size_t i = 0; i < FFT_SIZE / 2; i++)
51 {
52 if (i >= FFT_SIZE / 4)
53 out[i * 2] = out[i * 2 + 1] = 0;
54 }
55
56 FFT::fftsr_float (FFT_SIZE, out, in, FFT::PLAN_ESTIMATE);
57
58 // normalization, output
59 noise_96000.resize (FFT_SIZE);
60 for (size_t i = 0; i < FFT_SIZE; i++)
61 noise_96000[i] = in[i] / FFT_SIZE;
62
63 FFT::free_array_float (in);
64 FFT::free_array_float (out);
65 }
66
67 vector<float> out (5 * sr);
68 WavData wav_data (noise_96000, 1, 96000, 32);
69 MiniResampler mini_resampler (wav_data, 96000 / double (sr));
70 mini_resampler.read (0, out.size(), &out[0]);
71 return out;
72 }
73
74 void
avg_spectrum(const char * label,vector<float> & signal,int sr)75 avg_spectrum (const char *label, vector<float>& signal, int sr)
76 {
77 const size_t block_size = 2048;
78 size_t offset = 0;
79
80 vector<float> avg (block_size / 2 + 1);
81 vector<float> window (block_size);
82
83 for (guint i = 0; i < window.size(); i++)
84 window[i] = window_cos (2.0 * i / block_size - 1.0);
85
86 float *in = FFT::new_array_float (block_size);
87 float *out = FFT::new_array_float (block_size + 2);
88
89 while (offset + block_size < signal.size())
90 {
91 std::copy (signal.begin() + offset, signal.begin() + offset + block_size, in);
92 for (size_t i = 0; i < block_size; i++)
93 in[i] *= window[i];
94
95 FFT::fftar_float (block_size, in, out);
96
97 out[block_size] = out[1];
98 out[block_size+1] = 0;
99 out[1] = 0;
100
101 for (size_t d = 0; d < block_size + 2; d += 2)
102 avg[d/2] += out[d] * out[d] + out[d+1] * out[d+1];
103 offset += block_size / 4;
104 }
105
106 for (size_t i = 0; i < avg.size(); i++)
107 printf ("%s %g %g\n", label, i * double (sr) / block_size, avg[i]);
108 }
109
110
111 size_t
make_odd(size_t n)112 make_odd (size_t n)
113 {
114 if (n & 1)
115 return n;
116 return n - 1;
117 }
118
119 void
encode(vector<float> & audio_in,int sr,const string & win,float fundamental_freq)120 encode (vector<float>& audio_in, int sr, const string& win, float fundamental_freq)
121 {
122 EncoderParams enc_params;
123 enc_params.mix_freq = sr;
124 enc_params.zeropad = 4;
125 enc_params.fundamental_freq = fundamental_freq;
126 enc_params.frame_size_ms = 40;
127 enc_params.frame_size_ms = max<double> (enc_params.frame_size_ms, 1000 / enc_params.fundamental_freq * 4);
128 enc_params.frame_step_ms = enc_params.frame_size_ms / 4.0;
129 enc_params.frame_size = make_odd (enc_params.mix_freq * 0.001 * enc_params.frame_size_ms);
130 enc_params.frame_step = enc_params.mix_freq * 0.001 * enc_params.frame_step_ms;
131
132 /* compute block size from frame size (smallest 2^k value >= frame_size) */
133 uint64 block_size = 1;
134 while (block_size < enc_params.frame_size)
135 block_size *= 2;
136 enc_params.block_size = block_size;
137 vector<float> window (block_size);
138
139 for (guint i = 0; i < window.size(); i++)
140 {
141 if (i < enc_params.frame_size)
142 {
143 if (win == "rect")
144 {
145 window[i] = 1;
146 }
147 else if (win == "cos")
148 {
149 window[i] = window_cos (2.0 * i / enc_params.frame_size - 1.0);
150 }
151 else
152 {
153 assert (false); // not reached
154 }
155 }
156 else
157 window[i] = 0;
158 }
159 enc_params.window = window;
160
161 Encoder encoder (enc_params);
162
163 WavData wav_data (audio_in, 1, enc_params.mix_freq, 32);
164
165 const char *sm_file = "testnoisesr.tmp.sm";
166 encoder.encode (wav_data, 0, 1, /*attack*/ false, /*sines*/ false);
167 encoder.save (sm_file);
168
169 WavSet wav_set;
170
171 WavSetWave new_wave;
172 new_wave.midi_note = 60; // doesn't matter
173 new_wave.channel = 0;
174 new_wave.path = sm_file;
175 wav_set.waves.push_back (new_wave);
176
177 wav_set.save ("testnoisesr.tmp.smset", true);
178 }
179
180 void
decode(vector<float> & audio_out,int sr)181 decode (vector<float>& audio_out, int sr)
182 {
183 WavSet wav_set;
184 Error error = wav_set.load ("testnoisesr.tmp.smset");
185 assert (!error);
186
187 LiveDecoder decoder (&wav_set);
188
189 float freq = 440;
190 decoder.retrigger (0, freq, 127, sr);
191 decoder.process (audio_out.size(), nullptr, &audio_out[0]);
192 }
193
194 void
decode_resample(vector<float> & audio_out,int sr)195 decode_resample (vector<float>& audio_out, int sr)
196 {
197 vector<float> audio_high (sr * 5);
198
199 decode (audio_high, sr);
200
201 WavData wav_data (audio_high, 1, sr, 32);
202 MiniResampler mini_resampler (wav_data, double (sr) / 48000.);
203 mini_resampler.read (0, audio_out.size(), &audio_out[0]);
204 }
205
206 void
delete_tmp_files()207 delete_tmp_files()
208 {
209 unlink ("testnoisesr.tmp.smset");
210 unlink ("testnoisesr.tmp.sm");
211 }
212
213 double
energy(vector<float> & audio)214 energy (vector<float>& audio)
215 {
216 double e = 0;
217 for (size_t i = 0; i < audio.size(); i++)
218 e += audio[i] * audio[i];
219
220 return e;
221 }
222
223 vector<float>
get_noise_envelope()224 get_noise_envelope()
225 {
226 WavSet wav_set;
227 Error error = wav_set.load ("testnoisesr.tmp.smset");
228 assert (!error);
229
230 assert (wav_set.waves.size() == 1);
231 Audio *audio = wav_set.waves[0].audio;
232
233 size_t count = 0;
234 vector<float> noise_env (32);
235 for (size_t pos = 4; pos < audio->contents.size() - 4; pos++)
236 {
237 AudioBlock *block = &audio->contents[pos];
238 assert (noise_env.size() == block->noise.size());
239 for (size_t i = 0; i < block->noise.size(); i++)
240 noise_env[i] += block->noise_f (i);
241 count++;
242 }
243 for (size_t i = 0; i < noise_env.size(); i++)
244 noise_env[i] /= count;
245 return noise_env;
246 }
247
248 void
dump_noise_envelope()249 dump_noise_envelope()
250 {
251 vector<float> noise_env = get_noise_envelope();
252
253 for (size_t i = 0; i < noise_env.size(); i++)
254 printf ("noise-envelope %zd %.17g\n", i, noise_env[i]);
255 }
256
257 int
reencode(int argc,char ** argv)258 reencode (int argc, char **argv)
259 {
260 assert (argc >= 4);
261 int from_sr = atoi (argv[2]);
262 int to_sr = atoi (argv[3]);
263 string win = (argc >= 5) ? argv[4] : "cos";
264 float fundamental_freq = (argc >= 6) ? sm_atof (argv[5]) : 440;
265
266 printf ("# reencode from_sr = %d to_sr = %d, win = %s, fundamental = %f\n", from_sr, to_sr, win.c_str(), fundamental_freq);
267
268 Random random;
269
270 vector<float> noise (from_sr * 5);
271 for (size_t i = 0; i < noise.size(); i++)
272 noise[i] = random.random_double_range (-0.5, 0.5);
273
274 printf ("noise-from-level %.17g\n", sqrt (energy (noise) / noise.size()));
275
276 encode (noise, from_sr, win, fundamental_freq);
277 vector<float> from_env = get_noise_envelope();
278
279 vector<float> audio_out (to_sr * 5);
280 decode (audio_out, to_sr);
281
282 printf ("noise-to-level %.17g\n", sqrt (energy (audio_out) / audio_out.size()));
283
284 encode (audio_out, to_sr, win, fundamental_freq);
285 vector<float> to_env = get_noise_envelope();
286
287 assert (from_env.size() == to_env.size());
288 for (size_t i = 0; i < from_env.size(); i++)
289 {
290 printf ("noise-envelope %2zd %6.2f %.17g %.17g\n", i, to_env[i] / from_env[i] * 100., from_env[i], to_env[i]);
291 }
292
293 delete_tmp_files();
294 return 0;
295 }
296
297 int
main(int argc,char ** argv)298 main (int argc, char **argv)
299 {
300 Random random;
301
302 Main main (&argc, &argv);
303
304 // Debug::debug_enable ("encoder");
305
306 if (argc >= 2 && strcmp (argv[1], "reencode") == 0) // noise reencode test
307 {
308 return reencode (argc, argv);
309 }
310
311 // default noise sr test
312 int sr = (argc >= 2) ? atoi (argv[1]) : 48000;
313 string win = (argc >= 3) ? argv[2] : "cos";
314 float fundamental_freq = (argc >= 4) ? sm_atof (argv[3]) : 440;
315
316 printf ("# sr = %d, win = %s, fundamental_freq = %f\n", sr, win.c_str(), fundamental_freq);
317
318 NoiseGenerator noise_generator;
319 vector<float> noise;
320 vector<float> noise_48000;
321 if (sr == 48000)
322 {
323 noise.resize (sr * 5);
324 for (size_t i = 0; i < noise.size(); i++)
325 noise[i] = random.random_double_range (-0.5, 0.5);
326 noise_48000 = noise;
327 }
328 else
329 {
330 noise = noise_generator.gen_noise (sr);
331 noise_48000 = noise_generator.gen_noise (48000);
332 }
333
334 encode (noise, sr, win, fundamental_freq);
335 dump_noise_envelope();
336
337 printf ("noise-in-energy %.17g\n", energy (noise_48000));
338 printf ("noise-in-level %.17g\n", sqrt (energy (noise_48000) / noise_48000.size()));
339 printf ("noise-sr-level %.17g\n", sqrt (energy (noise) / noise.size()));
340
341 vector<float> audio_out (noise_48000.size());
342 for (int i = 48000; i < 197000; i += 12000)
343 {
344 decode_resample (audio_out, i);
345 printf ("noise-out-energy %d %.17g\n", i, energy (audio_out));
346 }
347 delete_tmp_files();
348 }
349