Line data Source code
1 : /******************************************************************************************************
2 :
3 : (C) 2022-2025 IVAS codec Public Collaboration with portions copyright Dolby International AB, Ericsson AB,
4 : Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD.,
5 : Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange,
6 : Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other
7 : contributors to this repository. All Rights Reserved.
8 :
9 : This software is protected by copyright law and by international treaties.
10 : The IVAS codec Public Collaboration consisting of Dolby International AB, Ericsson AB,
11 : Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD.,
12 : Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange,
13 : Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other
14 : contributors to this repository retain full ownership rights in their respective contributions in
15 : the software. This notice grants no license of any kind, including but not limited to patent
16 : license, nor is any license granted by implication, estoppel or otherwise.
17 :
18 : Contributors are required to enter into the IVAS codec Public Collaboration agreement before making
19 : contributions.
20 :
21 : This software is provided "AS IS", without any express or implied warranties. The software is in the
22 : development stage. It is intended exclusively for experts who have experience with such software and
23 : solely for the purpose of inspection. All implied warranties of non-infringement, merchantability
24 : and fitness for a particular purpose are hereby disclaimed and excluded.
25 :
26 : Any dispute, controversy or claim arising under or in relation to providing this software shall be
27 : submitted to and settled by the final, binding jurisdiction of the courts of Munich, Germany in
28 : accordance with the laws of the Federal Republic of Germany excluding its conflict of law rules and
29 : the United Nations Convention on Contracts on the International Sales of Goods.
30 :
31 : *******************************************************************************************************/
32 :
33 : /*====================================================================================
34 : EVS Codec 3GPP TS26.443 Nov 04, 2021. Version 12.14.0 / 13.10.0 / 14.6.0 / 15.4.0 / 16.3.0
35 : ====================================================================================*/
36 :
37 : #include "cnst.h"
38 : #include <stdint.h>
39 : #include "options.h"
40 : #ifdef DEBUGGING
41 : #include "debug.h"
42 : #endif
43 : #include <math.h>
44 : #include "prot.h"
45 : #include "wmc_auto.h"
46 :
47 : /*-------------------------------------------------------------------*
48 : * four1()
49 : *
50 : * From "numerical recipes in C".
51 : * Replace data by its DFT, if isign is input as 1; or replace data
52 : * by nn times its inverse-DFT, if isign is input as -1.
53 : * data is a complex array of length nn, input as a real
54 : * array data[1...2nn]. nn must be an integer power of 2
55 : *-------------------------------------------------------------------*/
56 :
57 0 : static void four1(
58 : float *data, /* i/o: data array .......... */
59 : int16_t nn, /* i : length of data array */
60 : int16_t isign /* i : sign +1 or -1 */
61 : )
62 : {
63 : int16_t n, mmax, m, j, istep, i;
64 : float wtemp, wr, wpr, wpi, wi, theta;
65 : float tempr, tempi;
66 :
67 0 : n = nn << 1;
68 0 : j = 1;
69 :
70 : /* this is the bit-reversal section of the routine */
71 0 : for ( i = 1; i < n; i += 2 )
72 : {
73 0 : if ( j > i )
74 : {
75 : /* exchange the two complex numbers */
76 0 : SWAP( data[j], data[i] );
77 0 : SWAP( data[j + 1], data[i + 1] );
78 : }
79 0 : m = n >> 1;
80 0 : while ( m >= 2 && j > m )
81 : {
82 0 : j -= m;
83 0 : m >>= 1;
84 : }
85 0 : j += m;
86 : }
87 0 : mmax = 2;
88 : /* here begins the Danielson-Lanczos section of the routine */
89 : /* Outer loop executed log2(nn) times */
90 0 : while ( n > mmax )
91 : {
92 0 : istep = 2 * mmax;
93 : /* initialization for the trigonometric recurrence */
94 0 : theta = (float) ( 6.28318530717959 / ( isign * mmax ) );
95 0 : wtemp = (float) ( sin( 0.5f * theta ) );
96 0 : wpr = -2.0f * wtemp * wtemp;
97 0 : wpi = (float) sin( theta );
98 0 : wr = 1.0f;
99 0 : wi = 0.0f;
100 : /* here are the two nested loops */
101 0 : for ( m = 1; m < mmax; m += 2 )
102 : {
103 0 : for ( i = m; i <= n; i += istep )
104 : {
105 : /* this is Danielson-Lanczos formula */
106 0 : j = i + mmax;
107 0 : tempr = wr * data[j] - wi * data[j + 1];
108 0 : tempi = wr * data[j + 1] + wi * data[j];
109 0 : data[j] = data[i] - tempr;
110 0 : data[j + 1] = data[i + 1] - tempi;
111 0 : data[i] += tempr;
112 0 : data[i + 1] += tempi;
113 : }
114 : /* trigonometric recurrence */
115 0 : wr = ( wtemp = wr ) * wpr - wi * wpi + wr;
116 0 : wi = wi * wpr + wtemp * wpi + wi;
117 : }
118 0 : mmax = istep;
119 : }
120 :
121 0 : return;
122 : }
123 :
124 : /*-------------------------------------------------------------------------*
125 : * realft()
126 : *
127 : * from "numerical recipes in C".
128 : * Calculates the Fourier Transform of a set of 2*n real-valued data points.
129 : * Replaces this data (which is stored in the array data[1..2n]) by the
130 : * positive frequancy half of its complex Fourier Transform. The real-valued
131 : * first and last components of the complex transform are returned as elements
132 : * data[1] and data[2] respectively. n must be a power of 2. This routine
133 : * also calculates the inverse transform of a complex data array if it is the
134 : * tranform of real data. (Results in this case must be multiplied by 1/n.)
135 : *--------------------------------------------------------------------------*/
136 :
137 0 : void realft(
138 : float *data, /* i/o: data array */
139 : int16_t n, /* i : length of data array */
140 : int16_t isign /* i : sign +1 or -1 */
141 : )
142 : {
143 : int16_t i, i1, i2, i3, i4, n2p3;
144 0 : float c1 = 0.5, c2, h1r, h1i, h2r, h2i;
145 : float wr, wi, wpr, wpi, wtemp, theta;
146 :
147 0 : theta = (float) ( EVS_PI / (float) n );
148 0 : if ( isign == 1 )
149 : {
150 : /* the forward transorm here */
151 0 : c2 = -0.5;
152 0 : four1( data, n, 1 );
153 : }
154 : else
155 : {
156 : /* otherwise set up for the inverse transform */
157 0 : c2 = 0.5;
158 0 : theta = -theta;
159 : }
160 0 : wtemp = (float) sin( 0.5f * theta );
161 0 : wpr = -2.0f * wtemp * wtemp;
162 0 : wpi = (float) sin( theta );
163 0 : wr = 1.0f + wpr;
164 0 : wi = wpi;
165 0 : n2p3 = 2 * n + 3;
166 : /* case i=1 done separately below */
167 0 : for ( i = 2; i <= n / 2; i++ )
168 : {
169 0 : i4 = 1 + ( i3 = n2p3 - ( i2 = 1 + ( i1 = i + i - 1 ) ) );
170 : /* the two separate transforms are separated out of data */
171 0 : h1r = c1 * ( data[i1] + data[i3] );
172 0 : h1i = c1 * ( data[i2] - data[i4] );
173 0 : h2r = -c2 * ( data[i2] + data[i4] );
174 0 : h2i = c2 * ( data[i1] - data[i3] );
175 : /* here they are recombined to form the true transform
176 : of the original real data */
177 0 : data[i1] = h1r + wr * h2r - wi * h2i;
178 0 : data[i2] = h1i + wr * h2i + wi * h2r;
179 0 : data[i3] = h1r - wr * h2r + wi * h2i;
180 0 : data[i4] = -h1i + wr * h2i + wi * h2r;
181 : /* the recurrence */
182 0 : wr = ( wtemp = wr ) * wpr - wi * wpi + wr;
183 0 : wi = wi * wpr + wtemp * wpi + wi;
184 : }
185 0 : if ( isign == 1 )
186 : {
187 : /* squeeze the first and the last data together to get them
188 : all within the original data */
189 0 : data[1] = ( h1r = data[1] ) + data[2];
190 0 : data[2] = h1r - data[2];
191 : }
192 : else
193 : {
194 : /* this is the inverse transform for the case isign=-1 */
195 0 : data[1] = c1 * ( ( h1r = data[1] ) + data[2] );
196 0 : data[2] = c1 * ( h1r - data[2] );
197 0 : four1( data, n, -1 );
198 : }
199 :
200 0 : return;
201 : }
|
