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_2M5_8K , T.W_2M5_16K, T.W_2M5_24K,
28            T.W_2M5_32K, T.W_2M5_48K, T.W_2M5_48K_HR, T.W_2M5_96K_HR ],
29
30          [ T.W_5M_8K  , T.W_5M_16K , T.W_5M_24K ,
31            T.W_5M_32K , T.W_5M_48K , T.W_5M_48K_HR , T.W_5M_96K_HR  ],
32
33          [ T.W_7M5_8K , T.W_7M5_16K, T.W_7M5_24K,
34            T.W_7M5_32K, T.W_7M5_48K, None, None ],
35
36          [ T.W_10M_8K , T.W_10M_16K, T.W_10M_24K,
37            T.W_10M_32K, T.W_10M_48K, T.W_10M_48K_HR, T.W_10M_96K_HR ] ]
38
39    def __init__(self, dt, sr):
40
41        self.ns = T.NS[dt][sr]
42        self.nd = T.ND[dt][sr]
43
44        self.t = np.zeros(2*self.ns)
45        self.w = Mdct.W[dt][sr]
46
47
48class MdctForward(Mdct):
49
50    def __init__(self, dt, sr):
51
52        super().__init__(dt, sr)
53
54    def run(self, x):
55
56        ns = self.ns
57        nd = self.nd
58
59        self.t[nd:nd+ns] = x
60        t = self.t * self.w
61        self.t[0:nd] = x[ns-nd:]
62
63        n  = len(t)
64        n2 = n // 2
65
66        z = t * np.exp(-2j * np.pi * np.arange(n) / (2*n))
67        z = scipy.fft.fft(z)[:n2]
68        z = z * np.exp(-2j * np.pi *
69                (n2/2 + 0.5) * (np.arange(n2) + 0.5) / (2 * n2))
70        return np.real(z) * np.sqrt(2/n2)
71
72
73class MdctInverse(Mdct):
74
75    def __init__(self, dt, sr):
76
77        super().__init__(dt, sr)
78
79    def run(self, x):
80
81        ns = self.ns
82        nd = self.nd
83
84        n = len(x)
85
86        x = np.append(x, -x[::-1])
87        z = x * np.exp(2j * np.pi * (n/2 + 0.5) * np.arange(2*n) / (2*n))
88        z = scipy.fft.ifft(z) * n
89        z = z * np.exp(2j * np.pi * (np.arange(2*n) + (n/2 + 0.5)) / (4*n))
90        t = np.real(z) * np.sqrt(2/n)
91
92        t = t * self.w[::-1]
93
94        y = np.empty(ns)
95        y[:nd] = t[ns-nd:ns] + self.t[2*ns-nd:]
96        y[nd:] = t[ns:2*ns-nd]
97        self.t = t
98
99        return y
100
101### ------------------------------------------------------------------------ ###
102
103def check_forward_unit(rng, dt, sr):
104
105    ns = T.NS[dt][sr]
106    nd = T.ND[dt][sr]
107    ok = True
108
109    x = (2 * rng.random(ns)) - 1
110
111    y   = [ None ] * 2
112    y_c = [ None ] * 2
113
114    mdct = MdctForward(dt, sr)
115    y[0] = mdct.run(x)
116    y[1] = mdct.run(x)
117
118    (y_c[0], d_c) = lc3.mdct_forward(dt, sr, x, np.zeros(nd))
119    y_c[1] = lc3.mdct_forward(dt, sr, x, d_c)[0]
120
121    ok = ok and np.amax(np.abs(y[0] - y_c[0])) < 1e-5
122    ok = ok and np.amax(np.abs(y[1] - y_c[1])) < 1e-5
123
124    return ok
125
126
127def check_forward_appendix_c(dt):
128
129    i0 = dt - T.DT_7M5
130    sr = T.SRATE_16K
131
132    ok = True
133
134    ns = T.NS[dt][sr]
135    nd = T.ND[dt][sr]
136
137    (y, d) = lc3.mdct_forward(dt, sr, C.X_PCM[i0][0], np.zeros(nd))
138    ok = ok and np.amax(np.abs(y - C.X[i0][0])) < 1e-1
139
140    (y, d) = lc3.mdct_forward(dt, sr, C.X_PCM[i0][1], d)
141    ok = ok and np.amax(np.abs(y - C.X[i0][1])) < 1e-1
142
143    return ok
144
145
146def check_inverse_unit(rng, dt, sr):
147
148    ns = T.NS[dt][sr]
149    nd = T.ND[dt][sr]
150    ok = True
151
152    x  = (2 * rng.random(ns)) - 1
153
154    y   = [ None ] * 2
155    y_c = [ None ] * 2
156
157    mdct = MdctInverse(dt, sr)
158    y[0] = mdct.run(x)
159    y[1] = mdct.run(x)
160
161    (y_c[0], d_c) = lc3.mdct_inverse(dt, sr, x, np.zeros(nd))
162    y_c[1] = lc3.mdct_inverse(dt, sr, x, d_c)[0]
163
164    ok = ok and np.amax(np.abs(y[0] - y_c[0])) < 1e-5
165    ok = ok and np.amax(np.abs(y[1] - y_c[1])) < 1e-5
166
167    return ok
168
169
170def check_inverse_appendix_c(dt):
171
172    i0 = dt - T.DT_7M5
173    sr = T.SRATE_16K
174
175    ok = True
176
177    ns = T.NS[dt][sr]
178    nd = T.ND[dt][sr]
179
180    (y, d0) = lc3.mdct_inverse(dt, sr, C.X_HAT_SNS[i0][0], np.zeros(nd))
181    yr = C.T_HAT_MDCT[i0][0][ns-nd:2*ns-nd]
182    dr = C.T_HAT_MDCT[i0][0][2*ns-nd:]
183
184    ok = ok and np.amax(np.abs(yr - y )) < 1e-1
185    ok = ok and np.amax(np.abs(dr - d0)) < 1e-1
186
187    (y, d1) = lc3.mdct_inverse(dt, sr, C.X_HAT_SNS[i0][1], d0)
188    yr[  :nd] = C.T_HAT_MDCT[i0][1][ns-nd:ns] + d0
189    yr[nd:ns] = C.T_HAT_MDCT[i0][1][ns:2*ns-nd]
190    dr        = C.T_HAT_MDCT[i0][1][2*ns-nd:]
191
192    ok = ok and np.amax(np.abs(yr - y )) < 1e-1
193    ok = ok and np.amax(np.abs(dr - d1)) < 1e-1
194
195    return ok
196
197
198def check():
199
200    rng = np.random.default_rng(1234)
201
202    ok  = True
203
204    for dt in range(T.NUM_DT):
205        for sr in range(T.SRATE_8K, T.SRATE_48K + 1):
206            ok = ok and check_forward_unit(rng, dt, sr)
207            ok = ok and check_inverse_unit(rng, dt, sr)
208
209    for dt in ( T.DT_2M5, T.DT_5M, T.DT_10M ):
210        for sr in ( T.SRATE_48K_HR, T.SRATE_96K_HR ):
211            ok = ok and check_forward_unit(rng, dt, sr)
212            ok = ok and check_inverse_unit(rng, dt, sr)
213
214    for dt in ( T.DT_7M5, T.DT_10M ):
215        ok = ok and check_forward_appendix_c(dt)
216        ok = ok and check_inverse_appendix_c(dt)
217
218    return ok
219