Line data Source code
1 : /*
2 : * Copyright (c) 2011 The WebRTC project authors. All Rights Reserved.
3 : *
4 : * Use of this source code is governed by a BSD-style license
5 : * that can be found in the LICENSE file in the root of the source
6 : * tree. An additional intellectual property rights grant can be found
7 : * in the file PATENTS. All contributing project authors may
8 : * be found in the AUTHORS file in the root of the source tree.
9 : */
10 :
11 :
12 : /*
13 : * This file contains the function WebRtcSpl_Sqrt().
14 : * The description header can be found in signal_processing_library.h
15 : *
16 : */
17 :
18 : #include "webrtc/base/checks.h"
19 : #include "webrtc/common_audio/signal_processing/include/signal_processing_library.h"
20 :
21 : int32_t WebRtcSpl_SqrtLocal(int32_t in);
22 :
23 0 : int32_t WebRtcSpl_SqrtLocal(int32_t in)
24 : {
25 :
26 : int16_t x_half, t16;
27 : int32_t A, B, x2;
28 :
29 : /* The following block performs:
30 : y=in/2
31 : x=y-2^30
32 : x_half=x/2^31
33 : t = 1 + (x_half) - 0.5*((x_half)^2) + 0.5*((x_half)^3) - 0.625*((x_half)^4)
34 : + 0.875*((x_half)^5)
35 : */
36 :
37 0 : B = in / 2;
38 :
39 0 : B = B - ((int32_t)0x40000000); // B = in/2 - 1/2
40 0 : x_half = (int16_t)(B >> 16); // x_half = x/2 = (in-1)/2
41 0 : B = B + ((int32_t)0x40000000); // B = 1 + x/2
42 0 : B = B + ((int32_t)0x40000000); // Add 0.5 twice (since 1.0 does not exist in Q31)
43 :
44 0 : x2 = ((int32_t)x_half) * ((int32_t)x_half) * 2; // A = (x/2)^2
45 0 : A = -x2; // A = -(x/2)^2
46 0 : B = B + (A >> 1); // B = 1 + x/2 - 0.5*(x/2)^2
47 :
48 0 : A >>= 16;
49 0 : A = A * A * 2; // A = (x/2)^4
50 0 : t16 = (int16_t)(A >> 16);
51 0 : B += -20480 * t16 * 2; // B = B - 0.625*A
52 : // After this, B = 1 + x/2 - 0.5*(x/2)^2 - 0.625*(x/2)^4
53 :
54 0 : A = x_half * t16 * 2; // A = (x/2)^5
55 0 : t16 = (int16_t)(A >> 16);
56 0 : B += 28672 * t16 * 2; // B = B + 0.875*A
57 : // After this, B = 1 + x/2 - 0.5*(x/2)^2 - 0.625*(x/2)^4 + 0.875*(x/2)^5
58 :
59 0 : t16 = (int16_t)(x2 >> 16);
60 0 : A = x_half * t16 * 2; // A = x/2^3
61 :
62 0 : B = B + (A >> 1); // B = B + 0.5*A
63 : // After this, B = 1 + x/2 - 0.5*(x/2)^2 + 0.5*(x/2)^3 - 0.625*(x/2)^4 + 0.875*(x/2)^5
64 :
65 0 : B = B + ((int32_t)32768); // Round off bit
66 :
67 0 : return B;
68 : }
69 :
70 0 : int32_t WebRtcSpl_Sqrt(int32_t value)
71 : {
72 : /*
73 : Algorithm:
74 :
75 : Six term Taylor Series is used here to compute the square root of a number
76 : y^0.5 = (1+x)^0.5 where x = y-1
77 : = 1+(x/2)-0.5*((x/2)^2+0.5*((x/2)^3-0.625*((x/2)^4+0.875*((x/2)^5)
78 : 0.5 <= x < 1
79 :
80 : Example of how the algorithm works, with ut=sqrt(in), and
81 : with in=73632 and ut=271 (even shift value case):
82 :
83 : in=73632
84 : y= in/131072
85 : x=y-1
86 : t = 1 + (x/2) - 0.5*((x/2)^2) + 0.5*((x/2)^3) - 0.625*((x/2)^4) + 0.875*((x/2)^5)
87 : ut=t*(1/sqrt(2))*512
88 :
89 : or:
90 :
91 : in=73632
92 : in2=73632*2^14
93 : y= in2/2^31
94 : x=y-1
95 : t = 1 + (x/2) - 0.5*((x/2)^2) + 0.5*((x/2)^3) - 0.625*((x/2)^4) + 0.875*((x/2)^5)
96 : ut=t*(1/sqrt(2))
97 : ut2=ut*2^9
98 :
99 : which gives:
100 :
101 : in = 73632
102 : in2 = 1206386688
103 : y = 0.56176757812500
104 : x = -0.43823242187500
105 : t = 0.74973506527313
106 : ut = 0.53014274874797
107 : ut2 = 2.714330873589594e+002
108 :
109 : or:
110 :
111 : in=73632
112 : in2=73632*2^14
113 : y=in2/2
114 : x=y-2^30
115 : x_half=x/2^31
116 : t = 1 + (x_half) - 0.5*((x_half)^2) + 0.5*((x_half)^3) - 0.625*((x_half)^4)
117 : + 0.875*((x_half)^5)
118 : ut=t*(1/sqrt(2))
119 : ut2=ut*2^9
120 :
121 : which gives:
122 :
123 : in = 73632
124 : in2 = 1206386688
125 : y = 603193344
126 : x = -470548480
127 : x_half = -0.21911621093750
128 : t = 0.74973506527313
129 : ut = 0.53014274874797
130 : ut2 = 2.714330873589594e+002
131 :
132 : */
133 :
134 : int16_t x_norm, nshift, t16, sh;
135 : int32_t A;
136 :
137 0 : int16_t k_sqrt_2 = 23170; // 1/sqrt2 (==5a82)
138 :
139 0 : A = value;
140 :
141 : // The convention in this function is to calculate sqrt(abs(A)). Negate the
142 : // input if it is negative.
143 0 : if (A < 0) {
144 0 : if (A == WEBRTC_SPL_WORD32_MIN) {
145 : // This number cannot be held in an int32_t after negating.
146 : // Map it to the maximum positive value.
147 0 : A = WEBRTC_SPL_WORD32_MAX;
148 : } else {
149 0 : A = -A;
150 : }
151 0 : } else if (A == 0) {
152 0 : return 0; // sqrt(0) = 0
153 : }
154 :
155 0 : sh = WebRtcSpl_NormW32(A); // # shifts to normalize A
156 0 : A = WEBRTC_SPL_LSHIFT_W32(A, sh); // Normalize A
157 0 : if (A < (WEBRTC_SPL_WORD32_MAX - 32767))
158 : {
159 0 : A = A + ((int32_t)32768); // Round off bit
160 : } else
161 : {
162 0 : A = WEBRTC_SPL_WORD32_MAX;
163 : }
164 :
165 0 : x_norm = (int16_t)(A >> 16); // x_norm = AH
166 :
167 0 : nshift = (sh / 2);
168 0 : RTC_DCHECK_GE(nshift, 0);
169 :
170 0 : A = (int32_t)WEBRTC_SPL_LSHIFT_W32((int32_t)x_norm, 16);
171 0 : A = WEBRTC_SPL_ABS_W32(A); // A = abs(x_norm<<16)
172 0 : A = WebRtcSpl_SqrtLocal(A); // A = sqrt(A)
173 :
174 0 : if (2 * nshift == sh) {
175 : // Even shift value case
176 :
177 0 : t16 = (int16_t)(A >> 16); // t16 = AH
178 :
179 0 : A = k_sqrt_2 * t16 * 2; // A = 1/sqrt(2)*t16
180 0 : A = A + ((int32_t)32768); // Round off
181 0 : A = A & ((int32_t)0x7fff0000); // Round off
182 :
183 0 : A >>= 15; // A = A>>16
184 :
185 : } else
186 : {
187 0 : A >>= 16; // A = A>>16
188 : }
189 :
190 0 : A = A & ((int32_t)0x0000ffff);
191 0 : A >>= nshift; // De-normalize the result.
192 :
193 0 : return A;
194 : }
|
