1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
|
/* GMR-1 AMBE vocoder - Math functions */
/* (C) 2011-2016 by Sylvain Munaut <tnt@246tNt.com>
* All Rights Reserved
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published by
* the Free Software Foundation; either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
/*! \addtogroup codec_private
* @{
*/
/*! \file codec/math.c
* \brief Osmocom GMR-1 AMBE vocoder math functions
*/
#include <math.h>
#include "private.h"
/*! \brief Table for \ref cosf_fast and \ref sinf_fast */
static float cos_tbl[1024];
/*! \brief Initializes \ref cos_tbl for \ref cosf_fast */
static void __attribute__ ((constructor))
cos_tbl_init(void)
{
int i;
for (i=0; i<1024; i++)
cos_tbl[i] = cosf((M_PIf * i) / 512.0f);
}
/*! \brief Fast Cosinus approximation using a simple table
* \param[in] angle The angle value
* \returns The cosinus of the angle
*/
float
cosf_fast(float angle)
{
const float f = 512.0f / M_PIf;
return cos_tbl[(int)(angle*f) & 1023];
}
/*! \brief Fast Sinus approximation using a simple table
* \param[in] angle The angle value
* \returns The sinus of the angle
*/
float
sinf_fast(float angle)
{
const float f = 512.0f / M_PIf;
return cos_tbl[((int)(angle*f) + 768) & 1023];
}
/*! \brief Forward Discrete Cosine Transform (fDCT)
* \param[out] out fDCT result buffer (freq domain, M elements)
* \param[in] in fDCT input buffer (time domain, N elements)
* \param[in] N Number of points of the DCT
* \param[in] M Limit to the number of frequency components (M <= N)
*/
void
ambe_fdct(float *out, float *in, int N, int M)
{
int i, j;
for (i=0; i<M; i++)
{
float v = 0.0f;
for (j=0; j<N; j++)
{
v += in[j] * cosf_fast( (M_PIf / N) * (j + .5f) * i );
}
out[i] = v / (float)N;
}
}
/*! \brief Inverse Discrete Cosine Transform (iDCT)
* \param[out] out iDCT result buffer (time domain, N elements)
* \param[in] in iDCT input buffer (freq domain, M elements)
* \param[in] N Number of points of the DCT
* \param[in] M Limit to the number of frequency components (M <= N)
*/
void
ambe_idct(float *out, float *in, int N, int M)
{
int i, j;
for (i=0; i<N; i++)
{
float v = in[0];
for (j=1; j<M; j++)
{
v += 2.0f * in[j] * cosf_fast( (M_PIf / N) * j * (i + .5f) );
}
out[i] = v;
}
}
/*! \brief Forward Discrete Fourrier Transform (float->complex)
* \param[out] out_i Real component result buffer (freq domain, N/2+1 elements)
* \param[out] out_q Imag component result buffer (freq domain, N/2+1 elements)
* \param[in] in Input buffer (time domain, M elements)
* \param[in] N Number of points of the DFT
* \param[in] M Limit to to the number of available time domain elements
*
* Since the input is float, the result is symmetric and so only one side
* is computed. The output index 0 is DC.
*/
void
ambe_fdft_fc(float *out_i, float *out_q, float *in, int N, int M)
{
int fb, ts;
for (fb=0; fb<=(N/2); fb++)
{
float i=0.0f, q=0.0f;
for (ts=0; ts<M; ts++)
{
float angle = (- 2.0f * M_PIf / N) * fb * ts;
i += in[ts] * cosf_fast(angle);
q += in[ts] * sinf_fast(angle);
}
out_i[fb] = i;
out_q[fb] = q;
}
}
/*! \brief Inverse Discret Fourrier Transform (complex->float)
* \param[out] out Result buffer (time domain, M
* \param[in] in_i Real component input buffer (freq domain, N/2+1 elements)
* \param[in] in_q Imag component input buffer (freq domain, N/2+1 elements)
* \param[in] N Number of points of the DFT
* \param[in] M Limit to the number of time domain elements to generate
*
* The input is assumed to be symmetric and so only N/2+1 inputs are
* needed. DC component must be input index 0.
*/
void
ambe_idft_cf(float *out, float *in_i, float *in_q, int N, int M)
{
int fb, ts;
for (ts=0; ts<M; ts++)
{
float r=0.0f;
for (fb=0; fb<=(N/2); fb++)
{
float angle = (- 2.0f * M_PIf / N) * fb * ts;
float m = (fb == 0 || fb == (N/2)) ? 1.0f : 2.0f;
r += m * (in_i[fb] * cosf_fast(angle) + in_q[fb] * sinf_fast(angle));
}
out[ts] = r / N;
}
}
/*! @} */
|
