1#
2# Copyright 2022 Google LLC
3#
4# Licensed under the Apache License, Version 2.0 (the "License");
5# you may not use this file except in compliance with the License.
6# You may obtain a copy of the License at
7#
8#     http://www.apache.org/licenses/LICENSE-2.0
9#
10# Unless required by applicable law or agreed to in writing, software
11# distributed under the License is distributed on an "AS IS" BASIS,
12# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13# See the License for the specific language governing permissions and
14# limitations under the License.
15#
16
17import numpy as np
18import scipy.fft
19
20import lc3
21import tables as T, appendix_c as C
22
23### ------------------------------------------------------------------------ ###
24
25class Mdct:
26
27    W = [ [ T.W_7M5_60, T.W_7M5_120, T.W_7M5_180, T.W_7M5_240, T.W_7M5_360 ],
28          [ T.W_10M_80, T.W_10M_160, T.W_10M_240, T.W_10M_320, T.W_10M_480 ] ]
29
30    def __init__(self, dt, sr):
31
32        self.ns = T.NS[dt][sr]
33        self.nd = T.ND[dt][sr]
34
35        self.t = np.zeros(2*self.ns)
36        self.w = Mdct.W[dt][sr]
37
38
39class MdctForward(Mdct):
40
41    def __init__(self, dt, sr):
42
43        super().__init__(dt, sr)
44
45    def run(self, x):
46
47        ns = self.ns
48        nd = self.nd
49
50        self.t[nd:nd+ns] = x
51        t = self.t * self.w
52        self.t[0:nd] = x[ns-nd:]
53
54        n  = len(t)
55        n2 = n // 2
56
57        z = t * np.exp(-2j * np.pi * np.arange(n) / (2*n))
58        z = scipy.fft.fft(z)[:n2]
59        z = z * np.exp(-2j * np.pi *
60                (n2/2 + 0.5) * (np.arange(n2) + 0.5) / (2 * n2))
61        return np.real(z) * np.sqrt(2/n2)
62
63
64class MdctInverse(Mdct):
65
66    def __init__(self, dt, sr):
67
68        super().__init__(dt, sr)
69
70    def run(self, x):
71
72        ns = self.ns
73        nd = self.nd
74
75        n = len(x)
76
77        x = np.append(x, -x[::-1])
78        z = x * np.exp(2j * np.pi * (n/2 + 0.5) * np.arange(2*n) / (2*n))
79        z = scipy.fft.ifft(z) * n
80        z = z * np.exp(2j * np.pi * (np.arange(2*n) + (n/2 + 0.5)) / (4*n))
81        t = np.real(z) * np.sqrt(2/n)
82
83        t = t * self.w[::-1]
84
85        y = np.empty(ns)
86        y[:nd] = t[ns-nd:ns] + self.t[2*ns-nd:]
87        y[nd:] = t[ns:2*ns-nd]
88        self.t = t
89
90        return y
91
92### ------------------------------------------------------------------------ ###
93
94def check_forward_unit(rng, dt, sr):
95
96    ns = T.NS[dt][sr]
97    nd = T.ND[dt][sr]
98    ok = True
99
100    x = (2 * rng.random(ns)) - 1
101
102    y   = [ None ] * 2
103    y_c = [ None ] * 2
104
105    mdct = MdctForward(dt, sr)
106    y[0] = mdct.run(x)
107    y[1] = mdct.run(x)
108
109    (y_c[0], d_c) = lc3.mdct_forward(dt, sr, x, np.zeros(nd))
110    y_c[1] = lc3.mdct_forward(dt, sr, x, d_c)[0]
111
112    ok = ok and np.amax(np.abs(y[0] - y_c[0])) < 1e-5
113    ok = ok and np.amax(np.abs(y[1] - y_c[1])) < 1e-5
114
115    return ok
116
117
118def check_forward_appendix_c(dt):
119
120    sr = T.SRATE_16K
121    ns = T.NS[dt][sr]
122    nd = T.ND[dt][sr]
123    ok = True
124
125    (y, d) = lc3.mdct_forward(dt, sr, C.X_PCM[dt][0], np.zeros(nd))
126    ok = ok and np.amax(np.abs(y - C.X[dt][0])) < 1e-1
127
128    (y, d) = lc3.mdct_forward(dt, sr, C.X_PCM[dt][1], d)
129    ok = ok and np.amax(np.abs(y - C.X[dt][1])) < 1e-1
130
131    return ok
132
133
134def check_inverse_unit(rng, dt, sr):
135
136    ns = T.NS[dt][sr]
137    nd = [ (23 * ns) // 30, (5 * ns) // 8 ][dt]
138    ok = True
139
140    x  = (2 * rng.random(ns)) - 1
141
142    y   = [ None ] * 2
143    y_c = [ None ] * 2
144
145    mdct = MdctInverse(dt, sr)
146    y[0] = mdct.run(x)
147    y[1] = mdct.run(x)
148
149    (y_c[0], d_c) = lc3.mdct_inverse(dt, sr, x, np.zeros(nd))
150    y_c[1] = lc3.mdct_inverse(dt, sr, x, d_c)[0]
151
152    ok = ok and np.amax(np.abs(y[0] - y_c[0])) < 1e-5
153    ok = ok and np.amax(np.abs(y[1] - y_c[1])) < 1e-5
154
155    return ok
156
157
158def check_inverse_appendix_c(dt):
159
160    sr = T.SRATE_16K
161    ns = T.NS[dt][sr]
162    nd = [ (23 * ns) // 30, (5 * ns) // 8 ][dt]
163    ok = True
164
165    (y, d0) = lc3.mdct_inverse(dt, sr, C.X_HAT_SNS[dt][0], np.zeros(nd))
166    yr = C.T_HAT_MDCT[dt][0][ns-nd:2*ns-nd]
167    dr = C.T_HAT_MDCT[dt][0][2*ns-nd:]
168
169    ok = ok and np.amax(np.abs(yr - y )) < 1e-1
170    ok = ok and np.amax(np.abs(dr - d0)) < 1e-1
171
172    (y, d1) = lc3.mdct_inverse(dt, sr, C.X_HAT_SNS[dt][1], d0)
173    yr[  :nd] = C.T_HAT_MDCT[dt][1][ns-nd:ns] + d0
174    yr[nd:ns] = C.T_HAT_MDCT[dt][1][ns:2*ns-nd]
175    dr        = C.T_HAT_MDCT[dt][1][2*ns-nd:]
176
177    ok = ok and np.amax(np.abs(yr - y )) < 1e-1
178    ok = ok and np.amax(np.abs(dr - d1)) < 1e-1
179
180    return ok
181
182
183def check():
184
185    rng = np.random.default_rng(1234)
186
187    ok  = True
188
189    for dt in range(T.NUM_DT):
190        for sr in range(T.NUM_SRATE):
191            ok = ok and check_forward_unit(rng, dt, sr)
192            ok = ok and check_inverse_unit(rng, dt, sr)
193
194    for dt in range(T.NUM_DT):
195        ok = ok and check_forward_appendix_c(dt)
196        ok = ok and check_inverse_appendix_c(dt)
197
198    return ok
199