1 // Licensed GNU LGPL v3 or later: http://www.gnu.org/licenses/lgpl.html
2
3 #include <string>
4 #include <vector>
5 #include "smfft.hh"
6 #include "smmain.hh"
7 #include "smpolyphaseinter.hh"
8 #include "smminiresampler.hh"
9 #include "smmath.hh"
10 #include <assert.h>
11 #include <math.h>
12 #include <string.h>
13
14 using std::string;
15 using std::vector;
16 using std::max;
17
18 using namespace SpectMorph;
19
20 void
get_error_signal(int high_sr,int sr,double freq,bool use_ppi,vector<float> & xout)21 get_error_signal (int high_sr, int sr, double freq, bool use_ppi, vector<float>& xout)
22 {
23 vector<float> audio_in (high_sr), out (sr);
24 for (size_t i = 0; i < audio_in.size(); i++)
25 {
26 audio_in[i] = sin (2 * M_PI * i * freq / high_sr);
27 }
28
29 if (use_ppi)
30 {
31 PolyPhaseInter *ppi = PolyPhaseInter::the();
32 for (size_t i = 0; i < out.size(); i++)
33 {
34 out[i] = ppi->get_sample (audio_in, double (high_sr) / double (sr) * i);
35 }
36 }
37 else
38 {
39 WavData wav_data (audio_in, 1, high_sr, 32);
40 MiniResampler mini_resampler (wav_data, double (high_sr) / sr);
41 mini_resampler.read (0, out.size(), &out[0]);
42 }
43 assert (xout.size() == 0);
44 for (size_t i = 100; i < out.size() - 100; i++)
45 {
46 double expect = sin (2 * M_PI * i * freq / sr);
47 xout.push_back (out[i] - expect);
48 }
49 }
50
51 void
scan(int high_sr,int sr,bool use_ppi)52 scan (int high_sr, int sr, bool use_ppi)
53 {
54 printf ("# scan: high_sr=%d, sr=%d, use_ppi=%s\n", high_sr, sr, use_ppi ? "true" : "false");
55
56 for (double freq = 100; freq < (sr/2.0 - 1.0); freq += 100)
57 {
58 vector<float> out;
59 get_error_signal (high_sr, sr, freq, use_ppi, out);
60
61 double max_diff = 0;
62 for (size_t i = 0; i < out.size(); i++)
63 {
64 max_diff = max (max_diff, double (out[i]));
65 }
66 double max_diff_db = 20 * log (max_diff) / log (10);
67 printf ("%.17g %.17g\n", freq, max_diff_db);
68 }
69 }
70
71 double
complex_abs(double re,double im)72 complex_abs (double re, double im)
73 {
74 return sqrt (re * re + im * im);
75 }
76
77 void
error_spectrum(int high_sr,int sr,double freq,bool use_ppi)78 error_spectrum (int high_sr, int sr, double freq, bool use_ppi)
79 {
80 vector<float> out;
81 get_error_signal (high_sr, sr, freq, use_ppi, out);
82
83 size_t FFT_SIZE = 32768;
84 float *fft_in = FFT::new_array_float (FFT_SIZE);
85 float *fft_out = FFT::new_array_float (FFT_SIZE);
86
87 double normalize = 0;
88
89 for (guint i = 0; i < FFT_SIZE; i++)
90 {
91 const double w = window_blackman ((2.0 * i - FFT_SIZE) / FFT_SIZE);
92
93 normalize += w;
94 fft_in[i] = out[i] * w;
95 }
96 normalize *= 0.5;
97
98 FFT::fftar_float (FFT_SIZE, fft_in, fft_out);
99
100 for (guint i = 0; i < FFT_SIZE/2; i++)
101 {
102 const double normalized_error = complex_abs (fft_out[i * 2], fft_out[i * 2 + 1]) / normalize;
103 const double normalized_error_db = 20 * log (normalized_error) / log (10);
104
105 printf ("%f %f\n", i / double (FFT_SIZE) * 44100, normalized_error_db);
106 }
107
108 FFT::free_array_float (fft_in);
109 FFT::free_array_float (fft_out);
110 }
111
112 int
main(int argc,char ** argv)113 main (int argc, char **argv)
114 {
115 Main main (&argc, &argv);
116
117 if (argc == 4 && strcmp (argv[1], "scan") == 0)
118 {
119 scan (atoi (argv[2]), atoi (argv[3]), false);
120 }
121 else if (argc == 4 && strcmp (argv[1], "scan-pp") == 0)
122 {
123 scan (atoi (argv[2]), atoi (argv[3]), true);
124 }
125 else if (argc == 5 && strcmp (argv[1], "error-spectrum") == 0)
126 {
127 error_spectrum (atoi (argv[2]), atoi (argv[3]), sm_atof (argv[4]), false);
128 }
129 else if (argc == 5 && strcmp (argv[1], "error-spectrum-pp") == 0)
130 {
131 error_spectrum (atoi (argv[2]), atoi (argv[3]), sm_atof (argv[4]), true);
132 }
133 }
134