a1773c987caa82cc98b607df7953cafd11f2ec52
[ffmpeg.git] / libavcodec / ps_tablegen.h
1 /*
2  * Header file for hardcoded Parametric Stereo tables
3  *
4  * Copyright (c) 2010 Alex Converse <alex.converse@gmail.com>
5  *
6  * This file is part of FFmpeg.
7  *
8  * FFmpeg is free software; you can redistribute it and/or
9  * modify it under the terms of the GNU Lesser General Public
10  * License as published by the Free Software Foundation; either
11  * version 2.1 of the License, or (at your option) any later version.
12  *
13  * FFmpeg is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
16  * Lesser General Public License for more details.
17  *
18  * You should have received a copy of the GNU Lesser General Public
19  * License along with FFmpeg; if not, write to the Free Software
20  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
21  */
22
23 #ifndef PS_TABLEGEN_H
24 #define PS_TABLEGEN_H
25
26 #include <stdint.h>
27 #include <math.h>
28
29 #if CONFIG_HARDCODED_TABLES
30 #define ps_tableinit()
31 #include "libavcodec/ps_tables.h"
32 #else
33 #include "../libavutil/common.h"
34 #ifndef M_SQRT1_2
35 #define M_SQRT1_2      0.70710678118654752440  /* 1/sqrt(2) */
36 #endif
37 #ifndef M_PI
38 #define M_PI           3.14159265358979323846  /* pi */
39 #endif
40 #ifndef M_SQRT2
41 #define M_SQRT2        1.41421356237309504880  /* sqrt(2) */
42 #endif
43 #define NR_ALLPASS_BANDS20 30
44 #define NR_ALLPASS_BANDS34 50
45 #define PS_AP_LINKS 3
46 static float pd_re_smooth[8*8*8];
47 static float pd_im_smooth[8*8*8];
48 static float HA[46][8][4];
49 static float HB[46][8][4];
50 static float f20_0_8 [ 8][7][2];
51 static float f34_0_12[12][7][2];
52 static float f34_1_8 [ 8][7][2];
53 static float f34_2_4 [ 4][7][2];
54 static float Q_fract_allpass[2][50][3][2];
55 static float phi_fract[2][50][2];
56
57 static const float g0_Q8[] = {
58     0.00746082949812f, 0.02270420949825f, 0.04546865930473f, 0.07266113929591f,
59     0.09885108575264f, 0.11793710567217f, 0.125f
60 };
61
62 static const float g0_Q12[] = {
63     0.04081179924692f, 0.03812810994926f, 0.05144908135699f, 0.06399831151592f,
64     0.07428313801106f, 0.08100347892914f, 0.08333333333333f
65 };
66
67 static const float g1_Q8[] = {
68     0.01565675600122f, 0.03752716391991f, 0.05417891378782f, 0.08417044116767f,
69     0.10307344158036f, 0.12222452249753f, 0.125f
70 };
71
72 static const float g2_Q4[] = {
73     -0.05908211155639f, -0.04871498374946f, 0.0f,   0.07778723915851f,
74      0.16486303567403f,  0.23279856662996f, 0.25f
75 };
76
77 static void make_filters_from_proto(float (*filter)[7][2], const float *proto, int bands)
78 {
79     int q, n;
80     for (q = 0; q < bands; q++) {
81         for (n = 0; n < 7; n++) {
82             double theta = 2 * M_PI * (q + 0.5) * (n - 6) / bands;
83             filter[q][n][0] = proto[n] *  cos(theta);
84             filter[q][n][1] = proto[n] * -sin(theta);
85         }
86     }
87 }
88
89 static void ps_tableinit(void)
90 {
91     static const float ipdopd_sin[] = { 0, M_SQRT1_2, 1,  M_SQRT1_2,  0, -M_SQRT1_2, -1, -M_SQRT1_2 };
92     static const float ipdopd_cos[] = { 1, M_SQRT1_2, 0, -M_SQRT1_2, -1, -M_SQRT1_2,  0,  M_SQRT1_2 };
93     int pd0, pd1, pd2;
94
95     static const float iid_par_dequant[] = {
96         //iid_par_dequant_default
97         0.05623413251903, 0.12589254117942, 0.19952623149689, 0.31622776601684,
98         0.44668359215096, 0.63095734448019, 0.79432823472428, 1,
99         1.25892541179417, 1.58489319246111, 2.23872113856834, 3.16227766016838,
100         5.01187233627272, 7.94328234724282, 17.7827941003892,
101         //iid_par_dequant_fine
102         0.00316227766017, 0.00562341325190, 0.01,             0.01778279410039,
103         0.03162277660168, 0.05623413251903, 0.07943282347243, 0.11220184543020,
104         0.15848931924611, 0.22387211385683, 0.31622776601684, 0.39810717055350,
105         0.50118723362727, 0.63095734448019, 0.79432823472428, 1,
106         1.25892541179417, 1.58489319246111, 1.99526231496888, 2.51188643150958,
107         3.16227766016838, 4.46683592150963, 6.30957344480193, 8.91250938133745,
108         12.5892541179417, 17.7827941003892, 31.6227766016838, 56.2341325190349,
109         100,              177.827941003892, 316.227766016837,
110     };
111     static const float icc_invq[] = {
112         1, 0.937,      0.84118,    0.60092,    0.36764,   0,      -0.589,    -1
113     };
114     static const float acos_icc_invq[] = {
115         0, 0.35685527, 0.57133466, 0.92614472, 1.1943263, M_PI/2, 2.2006171, M_PI
116     };
117     int iid, icc;
118
119     int k, m;
120     static const int8_t f_center_20[] = {
121         -3, -1, 1, 3, 5, 7, 10, 14, 18, 22,
122     };
123     static const int8_t f_center_34[] = {
124          2,  6, 10, 14, 18, 22, 26, 30,
125         34,-10, -6, -2, 51, 57, 15, 21,
126         27, 33, 39, 45, 54, 66, 78, 42,
127        102, 66, 78, 90,102,114,126, 90,
128     };
129     static const float fractional_delay_links[] = { 0.43f, 0.75f, 0.347f };
130     const float fractional_delay_gain = 0.39f;
131
132     for (pd0 = 0; pd0 < 8; pd0++) {
133         float pd0_re = ipdopd_cos[pd0];
134         float pd0_im = ipdopd_sin[pd0];
135         for (pd1 = 0; pd1 < 8; pd1++) {
136             float pd1_re = ipdopd_cos[pd1];
137             float pd1_im = ipdopd_sin[pd1];
138             for (pd2 = 0; pd2 < 8; pd2++) {
139                 float pd2_re = ipdopd_cos[pd2];
140                 float pd2_im = ipdopd_sin[pd2];
141                 float re_smooth = 0.25f * pd0_re + 0.5f * pd1_re + pd2_re;
142                 float im_smooth = 0.25f * pd0_im + 0.5f * pd1_im + pd2_im;
143                 float pd_mag = 1 / sqrt(im_smooth * im_smooth + re_smooth * re_smooth);
144                 pd_re_smooth[pd0*64+pd1*8+pd2] = re_smooth * pd_mag;
145                 pd_im_smooth[pd0*64+pd1*8+pd2] = im_smooth * pd_mag;
146             }
147         }
148     }
149
150     for (iid = 0; iid < 46; iid++) {
151         float c = iid_par_dequant[iid]; //<Linear Inter-channel Intensity Difference
152         float c1 = (float)M_SQRT2 / sqrtf(1.0f + c*c);
153         float c2 = c * c1;
154         for (icc = 0; icc < 8; icc++) {
155             /*if (PS_BASELINE || ps->icc_mode < 3)*/ {
156                 float alpha = 0.5f * acos_icc_invq[icc];
157                 float beta  = alpha * (c1 - c2) * (float)M_SQRT1_2;
158                 HA[iid][icc][0] = c2 * cosf(beta + alpha);
159                 HA[iid][icc][1] = c1 * cosf(beta - alpha);
160                 HA[iid][icc][2] = c2 * sinf(beta + alpha);
161                 HA[iid][icc][3] = c1 * sinf(beta - alpha);
162             } /* else */ {
163                 float alpha, gamma, mu, rho;
164                 float alpha_c, alpha_s, gamma_c, gamma_s;
165                 rho = FFMAX(icc_invq[icc], 0.05f);
166                 alpha = 0.5f * atan2f(2.0f * c * rho, c*c - 1.0f);
167                 mu = c + 1.0f / c;
168                 mu = sqrtf(1 + (4 * rho * rho - 4)/(mu * mu));
169                 gamma = atanf(sqrtf((1.0f - mu)/(1.0f + mu)));
170                 if (alpha < 0) alpha += M_PI/2;
171                 alpha_c = cosf(alpha);
172                 alpha_s = sinf(alpha);
173                 gamma_c = cosf(gamma);
174                 gamma_s = sinf(gamma);
175                 HB[iid][icc][0] =  M_SQRT2 * alpha_c * gamma_c;
176                 HB[iid][icc][1] =  M_SQRT2 * alpha_s * gamma_c;
177                 HB[iid][icc][2] = -M_SQRT2 * alpha_s * gamma_s;
178                 HB[iid][icc][3] =  M_SQRT2 * alpha_c * gamma_s;
179             }
180         }
181     }
182
183     for (k = 0; k < NR_ALLPASS_BANDS20; k++) {
184         double f_center, theta;
185         if (k < FF_ARRAY_ELEMS(f_center_20))
186             f_center = f_center_20[k] * 0.125;
187         else
188             f_center = k - 6.5f;
189         for (m = 0; m < PS_AP_LINKS; m++) {
190             theta = -M_PI * fractional_delay_links[m] * f_center;
191             Q_fract_allpass[0][k][m][0] = cos(theta);
192             Q_fract_allpass[0][k][m][1] = sin(theta);
193         }
194         theta = -M_PI*fractional_delay_gain*f_center;
195         phi_fract[0][k][0] = cos(theta);
196         phi_fract[0][k][1] = sin(theta);
197     }
198     for (k = 0; k < NR_ALLPASS_BANDS34; k++) {
199         double f_center, theta;
200         if (k < FF_ARRAY_ELEMS(f_center_34))
201             f_center = f_center_34[k] / 24.;
202         else
203             f_center = k - 26.5f;
204         for (m = 0; m < PS_AP_LINKS; m++) {
205             theta = -M_PI * fractional_delay_links[m] * f_center;
206             Q_fract_allpass[1][k][m][0] = cos(theta);
207             Q_fract_allpass[1][k][m][1] = sin(theta);
208         }
209         theta = -M_PI*fractional_delay_gain*f_center;
210         phi_fract[1][k][0] = cos(theta);
211         phi_fract[1][k][1] = sin(theta);
212     }
213
214     make_filters_from_proto(f20_0_8,  g0_Q8,   8);
215     make_filters_from_proto(f34_0_12, g0_Q12, 12);
216     make_filters_from_proto(f34_1_8,  g1_Q8,   8);
217     make_filters_from_proto(f34_2_4,  g2_Q4,   4);
218 }
219 #endif /* CONFIG_HARDCODED_TABLES */
220
221 #endif /* PS_TABLEGEN_H */