bd4e695568809345bfa23be63a6c9fe077a40fce
[ffmpeg.git] / libavcodec / aacps_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 Libav.
7  *
8  * Libav 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  * Libav 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 Libav; if not, write to the Free Software
20  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
21  */
22
23 #ifndef AACPS_TABLEGEN_H
24 #define AACPS_TABLEGEN_H
25
26 #include <math.h>
27 #include <stdint.h>
28
29 #if CONFIG_HARDCODED_TABLES
30 #define ps_tableinit()
31 #include "libavcodec/aacps_tables.h"
32 #else
33 #include "libavutil/common.h"
34 #include "libavutil/mathematics.h"
35 #include "libavutil/mem.h"
36 #define NR_ALLPASS_BANDS20 30
37 #define NR_ALLPASS_BANDS34 50
38 #define PS_AP_LINKS 3
39 static float pd_re_smooth[8*8*8];
40 static float pd_im_smooth[8*8*8];
41 static float HA[46][8][4];
42 static float HB[46][8][4];
43 static DECLARE_ALIGNED(16, float, f20_0_8) [ 8][8][2];
44 static DECLARE_ALIGNED(16, float, f34_0_12)[12][8][2];
45 static DECLARE_ALIGNED(16, float, f34_1_8) [ 8][8][2];
46 static DECLARE_ALIGNED(16, float, f34_2_4) [ 4][8][2];
47 static DECLARE_ALIGNED(16, float, Q_fract_allpass)[2][50][3][2];
48 static DECLARE_ALIGNED(16, float, phi_fract)[2][50][2];
49
50 static const float g0_Q8[] = {
51     0.00746082949812f, 0.02270420949825f, 0.04546865930473f, 0.07266113929591f,
52     0.09885108575264f, 0.11793710567217f, 0.125f
53 };
54
55 static const float g0_Q12[] = {
56     0.04081179924692f, 0.03812810994926f, 0.05144908135699f, 0.06399831151592f,
57     0.07428313801106f, 0.08100347892914f, 0.08333333333333f
58 };
59
60 static const float g1_Q8[] = {
61     0.01565675600122f, 0.03752716391991f, 0.05417891378782f, 0.08417044116767f,
62     0.10307344158036f, 0.12222452249753f, 0.125f
63 };
64
65 static const float g2_Q4[] = {
66     -0.05908211155639f, -0.04871498374946f, 0.0f,   0.07778723915851f,
67      0.16486303567403f,  0.23279856662996f, 0.25f
68 };
69
70 static void make_filters_from_proto(float (*filter)[8][2], const float *proto, int bands)
71 {
72     int q, n;
73     for (q = 0; q < bands; q++) {
74         for (n = 0; n < 7; n++) {
75             double theta = 2 * M_PI * (q + 0.5) * (n - 6) / bands;
76             filter[q][n][0] = proto[n] *  cos(theta);
77             filter[q][n][1] = proto[n] * -sin(theta);
78         }
79     }
80 }
81
82 static void ps_tableinit(void)
83 {
84     static const float ipdopd_sin[] = { 0, M_SQRT1_2, 1,  M_SQRT1_2,  0, -M_SQRT1_2, -1, -M_SQRT1_2 };
85     static const float ipdopd_cos[] = { 1, M_SQRT1_2, 0, -M_SQRT1_2, -1, -M_SQRT1_2,  0,  M_SQRT1_2 };
86     int pd0, pd1, pd2;
87
88     static const float iid_par_dequant[] = {
89         //iid_par_dequant_default
90         0.05623413251903, 0.12589254117942, 0.19952623149689, 0.31622776601684,
91         0.44668359215096, 0.63095734448019, 0.79432823472428, 1,
92         1.25892541179417, 1.58489319246111, 2.23872113856834, 3.16227766016838,
93         5.01187233627272, 7.94328234724282, 17.7827941003892,
94         //iid_par_dequant_fine
95         0.00316227766017, 0.00562341325190, 0.01,             0.01778279410039,
96         0.03162277660168, 0.05623413251903, 0.07943282347243, 0.11220184543020,
97         0.15848931924611, 0.22387211385683, 0.31622776601684, 0.39810717055350,
98         0.50118723362727, 0.63095734448019, 0.79432823472428, 1,
99         1.25892541179417, 1.58489319246111, 1.99526231496888, 2.51188643150958,
100         3.16227766016838, 4.46683592150963, 6.30957344480193, 8.91250938133745,
101         12.5892541179417, 17.7827941003892, 31.6227766016838, 56.2341325190349,
102         100,              177.827941003892, 316.227766016837,
103     };
104     static const float icc_invq[] = {
105         1, 0.937,      0.84118,    0.60092,    0.36764,   0,      -0.589,    -1
106     };
107     static const float acos_icc_invq[] = {
108         0, 0.35685527, 0.57133466, 0.92614472, 1.1943263, M_PI/2, 2.2006171, M_PI
109     };
110     int iid, icc;
111
112     int k, m;
113     static const int8_t f_center_20[] = {
114         -3, -1, 1, 3, 5, 7, 10, 14, 18, 22,
115     };
116     static const int8_t f_center_34[] = {
117          2,  6, 10, 14, 18, 22, 26, 30,
118         34,-10, -6, -2, 51, 57, 15, 21,
119         27, 33, 39, 45, 54, 66, 78, 42,
120        102, 66, 78, 90,102,114,126, 90,
121     };
122     static const float fractional_delay_links[] = { 0.43f, 0.75f, 0.347f };
123     const float fractional_delay_gain = 0.39f;
124
125     for (pd0 = 0; pd0 < 8; pd0++) {
126         float pd0_re = ipdopd_cos[pd0];
127         float pd0_im = ipdopd_sin[pd0];
128         for (pd1 = 0; pd1 < 8; pd1++) {
129             float pd1_re = ipdopd_cos[pd1];
130             float pd1_im = ipdopd_sin[pd1];
131             for (pd2 = 0; pd2 < 8; pd2++) {
132                 float pd2_re = ipdopd_cos[pd2];
133                 float pd2_im = ipdopd_sin[pd2];
134                 float re_smooth = 0.25f * pd0_re + 0.5f * pd1_re + pd2_re;
135                 float im_smooth = 0.25f * pd0_im + 0.5f * pd1_im + pd2_im;
136                 float pd_mag = 1 / sqrt(im_smooth * im_smooth + re_smooth * re_smooth);
137                 pd_re_smooth[pd0*64+pd1*8+pd2] = re_smooth * pd_mag;
138                 pd_im_smooth[pd0*64+pd1*8+pd2] = im_smooth * pd_mag;
139             }
140         }
141     }
142
143     for (iid = 0; iid < 46; iid++) {
144         float c = iid_par_dequant[iid]; ///< Linear Inter-channel Intensity Difference
145         float c1 = (float)M_SQRT2 / sqrtf(1.0f + c*c);
146         float c2 = c * c1;
147         for (icc = 0; icc < 8; icc++) {
148             /*if (PS_BASELINE || ps->icc_mode < 3)*/ {
149                 float alpha = 0.5f * acos_icc_invq[icc];
150                 float beta  = alpha * (c1 - c2) * (float)M_SQRT1_2;
151                 HA[iid][icc][0] = c2 * cosf(beta + alpha);
152                 HA[iid][icc][1] = c1 * cosf(beta - alpha);
153                 HA[iid][icc][2] = c2 * sinf(beta + alpha);
154                 HA[iid][icc][3] = c1 * sinf(beta - alpha);
155             } /* else */ {
156                 float alpha, gamma, mu, rho;
157                 float alpha_c, alpha_s, gamma_c, gamma_s;
158                 rho = FFMAX(icc_invq[icc], 0.05f);
159                 alpha = 0.5f * atan2f(2.0f * c * rho, c*c - 1.0f);
160                 mu = c + 1.0f / c;
161                 mu = sqrtf(1 + (4 * rho * rho - 4)/(mu * mu));
162                 gamma = atanf(sqrtf((1.0f - mu)/(1.0f + mu)));
163                 if (alpha < 0) alpha += M_PI/2;
164                 alpha_c = cosf(alpha);
165                 alpha_s = sinf(alpha);
166                 gamma_c = cosf(gamma);
167                 gamma_s = sinf(gamma);
168                 HB[iid][icc][0] =  M_SQRT2 * alpha_c * gamma_c;
169                 HB[iid][icc][1] =  M_SQRT2 * alpha_s * gamma_c;
170                 HB[iid][icc][2] = -M_SQRT2 * alpha_s * gamma_s;
171                 HB[iid][icc][3] =  M_SQRT2 * alpha_c * gamma_s;
172             }
173         }
174     }
175
176     for (k = 0; k < NR_ALLPASS_BANDS20; k++) {
177         double f_center, theta;
178         if (k < FF_ARRAY_ELEMS(f_center_20))
179             f_center = f_center_20[k] * 0.125;
180         else
181             f_center = k - 6.5f;
182         for (m = 0; m < PS_AP_LINKS; m++) {
183             theta = -M_PI * fractional_delay_links[m] * f_center;
184             Q_fract_allpass[0][k][m][0] = cos(theta);
185             Q_fract_allpass[0][k][m][1] = sin(theta);
186         }
187         theta = -M_PI*fractional_delay_gain*f_center;
188         phi_fract[0][k][0] = cos(theta);
189         phi_fract[0][k][1] = sin(theta);
190     }
191     for (k = 0; k < NR_ALLPASS_BANDS34; k++) {
192         double f_center, theta;
193         if (k < FF_ARRAY_ELEMS(f_center_34))
194             f_center = f_center_34[k] / 24.;
195         else
196             f_center = k - 26.5f;
197         for (m = 0; m < PS_AP_LINKS; m++) {
198             theta = -M_PI * fractional_delay_links[m] * f_center;
199             Q_fract_allpass[1][k][m][0] = cos(theta);
200             Q_fract_allpass[1][k][m][1] = sin(theta);
201         }
202         theta = -M_PI*fractional_delay_gain*f_center;
203         phi_fract[1][k][0] = cos(theta);
204         phi_fract[1][k][1] = sin(theta);
205     }
206
207     make_filters_from_proto(f20_0_8,  g0_Q8,   8);
208     make_filters_from_proto(f34_0_12, g0_Q12, 12);
209     make_filters_from_proto(f34_1_8,  g1_Q8,   8);
210     make_filters_from_proto(f34_2_4,  g2_Q4,   4);
211 }
212 #endif /* CONFIG_HARDCODED_TABLES */
213
214 #endif /* AACPS_TABLEGEN_H */