From 16e60160f8b733000603eba0eb090887c1fc1782 Mon Sep 17 00:00:00 2001 From: Hui Zhang Date: Fri, 17 Sep 2021 02:36:00 -0500 Subject: [PATCH] Kaldi (#839) * can do frames, real stft * format * stft complex, powspec, magspec * add common utils * add window process func * using frames and matmul as stft * read with 2d; window process * test with dither, remove dc offset, preermphs * add doc string * more frontend utils * add logspec * fix typing * add delpoy mergify label --- deepspeech/io/dataset.py | 12 +- third_party/__init__.py | 0 third_party/paddle_audio/__init__.py | 0 third_party/paddle_audio/frontend.py | 146 ----- third_party/paddle_audio/frontend/common.py | 201 +++++++ third_party/paddle_audio/frontend/english.wav | Bin 0 -> 35824 bytes third_party/paddle_audio/frontend/kaldi.py | 266 +++++++++ .../paddle_audio/frontend/kaldi_test.py | 533 ++++++++++++++++++ 8 files changed, 1006 insertions(+), 152 deletions(-) create mode 100644 third_party/__init__.py create mode 100644 third_party/paddle_audio/__init__.py delete mode 100644 third_party/paddle_audio/frontend.py create mode 100644 third_party/paddle_audio/frontend/common.py create mode 100644 third_party/paddle_audio/frontend/english.wav create mode 100644 third_party/paddle_audio/frontend/kaldi.py create mode 100644 third_party/paddle_audio/frontend/kaldi_test.py diff --git a/deepspeech/io/dataset.py b/deepspeech/io/dataset.py index d1fe0470..e58e03b4 100644 --- a/deepspeech/io/dataset.py +++ b/deepspeech/io/dataset.py @@ -76,19 +76,19 @@ class ManifestDataset(Dataset): Args: manifest_path (str): manifest josn file path - max_input_len ([type], optional): maximum output seq length, + max_input_len ([type], optional): maximum output seq length, in seconds for raw wav, in frame numbers for feature data. Defaults to float('inf'). - min_input_len (float, optional): minimum input seq length, + min_input_len (float, optional): minimum input seq length, in seconds for raw wav, in frame numbers for feature data. Defaults to 0.0. - max_output_len (float, optional): maximum input seq length, + max_output_len (float, optional): maximum input seq length, in modeling units. Defaults to 500.0. - min_output_len (float, optional): minimum input seq length, + min_output_len (float, optional): minimum input seq length, in modeling units. Defaults to 0.0. - max_output_input_ratio (float, optional): maximum output seq length/output seq length ratio. + max_output_input_ratio (float, optional): maximum output seq length/output seq length ratio. Defaults to 10.0. min_output_input_ratio (float, optional): minimum output seq length/output seq length ratio. Defaults to 0.05. - + """ super().__init__() diff --git a/third_party/__init__.py b/third_party/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/third_party/paddle_audio/__init__.py b/third_party/paddle_audio/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/third_party/paddle_audio/frontend.py b/third_party/paddle_audio/frontend.py deleted file mode 100644 index 1b337732..00000000 --- a/third_party/paddle_audio/frontend.py +++ /dev/null @@ -1,146 +0,0 @@ -from typing import Tuple -import numpy as np -import paddle -from paddle import Tensor -from paddle import nn -from paddle.nn import functional as F - - -def frame(x: Tensor, - num_samples: Tensor, - win_length: int, - hop_length: int, - clip: bool = True) -> Tuple[Tensor, Tensor]: - """Extract frames from audio. - - Parameters - ---------- - x : Tensor - Shape (N, T), batched waveform. - num_samples : Tensor - Shape (N, ), number of samples of each waveform. - win_length : int - Window length. - hop_length : int - Number of samples shifted between ajancent frames. - clip : bool, optional - Whether to clip audio that does not fit into the last frame, by - default True - - Returns - ------- - frames : Tensor - Shape (N, T', win_length). - num_frames : Tensor - Shape (N, ) number of valid frames - """ - assert hop_length <= win_length - num_frames = (num_samples - win_length) // hop_length - padding = (0, 0) - if not clip: - num_frames += 1 - # NOTE: pad hop_length - 1 to the right to ensure that there is at most - # one frame dangling to the righe edge - padding = (0, hop_length - 1) - - weight = paddle.eye(win_length).unsqueeze(1) - - frames = F.conv1d(x.unsqueeze(1), - weight, - padding=padding, - stride=(hop_length, )) - return frames, num_frames - - -class STFT(nn.Layer): - """A module for computing stft transformation in a differentiable way. - - Parameters - ------------ - n_fft : int - Number of samples in a frame. - - hop_length : int - Number of samples shifted between adjacent frames. - - win_length : int - Length of the window. - - clip: bool - Whether to clip audio is necesaary. - """ - def __init__(self, - n_fft: int, - hop_length: int, - win_length: int, - window_type: str = None, - clip: bool = True): - super().__init__() - - self.hop_length = hop_length - self.n_bin = 1 + n_fft // 2 - self.n_fft = n_fft - self.clip = clip - - # calculate window - if window_type is None: - window = np.ones(win_length) - elif window_type == "hann": - window = np.hanning(win_length) - elif window_type == "hamming": - window = np.hamming(win_length) - else: - raise ValueError("Not supported yet!") - - if win_length < n_fft: - window = F.pad(window, (0, n_fft - win_length)) - elif win_length > n_fft: - window = window[:n_fft] - - # (n_bins, n_fft) complex - kernel_size = min(n_fft, win_length) - weight = np.fft.fft(np.eye(n_fft))[:self.n_bin, :kernel_size] - w_real = weight.real - w_imag = weight.imag - - # (2 * n_bins, kernel_size) - w = np.concatenate([w_real, w_imag], axis=0) - w = w * window - - # (2 * n_bins, 1, kernel_size) # (C_out, C_in, kernel_size) - w = np.expand_dims(w, 1) - weight = paddle.cast(paddle.to_tensor(w), paddle.get_default_dtype()) - self.register_buffer("weight", weight) - - def forward(self, x: Tensor, num_samples: Tensor) -> Tuple[Tensor, Tensor]: - """Compute the stft transform. - Parameters - ------------ - x : Tensor [shape=(B, T)] - The input waveform. - num_samples : Tensor - Number of samples of each waveform. - Returns - ------------ - D : Tensor - Shape(N, T', n_bins, 2) Spectrogram. - - num_frames: Tensor - Shape (N,) number of samples of each spectrogram - """ - num_frames = (num_samples - self.win_length) // self.hop_length - padding = (0, 0) - if not self.clip: - num_frames += 1 - padding = (0, self.hop_length - 1) - - batch_size, _, _ = paddle.shape(x) - x = x.unsqueeze(-1) - D = F.conv1d(self.weight, - x, - stride=(self.hop_length, ), - padding=padding, - data_format="NLC") - D = paddle.reshape(D, [batch_size, -1, self.n_bin, 2]) - return D, num_frames - diff --git a/third_party/paddle_audio/frontend/common.py b/third_party/paddle_audio/frontend/common.py new file mode 100644 index 00000000..7638dae5 --- /dev/null +++ b/third_party/paddle_audio/frontend/common.py @@ -0,0 +1,201 @@ +import paddle +import numpy as np +from typing import Tuple, Optional, Union + + +# https://github.com/kaldi-asr/kaldi/blob/cbed4ff688/src/feat/feature-window.cc#L109 +def povey_window(frame_len:int) -> np.ndarray: + win = np.empty(frame_len) + a = 2 * np.pi / (frame_len -1) + for i in range(frame_len): + win[i] = (0.5 - 0.5 * np.cos(a * i) )**0.85 + return win + +def hann_window(frame_len:int) -> np.ndarray: + win = np.empty(frame_len) + a = 2 * np.pi / (frame_len -1) + for i in range(frame_len): + win[i] = 0.5 - 0.5 * np.cos(a * i) + return win + +def sine_window(frame_len:int) -> np.ndarray: + win = np.empty(frame_len) + a = 2 * np.pi / (frame_len -1) + for i in range(frame_len): + win[i] = np.sin(0.5 * a * i) + return win + +def hamm_window(frame_len:int) -> np.ndarray: + win = np.empty(frame_len) + a = 2 * np.pi / (frame_len -1) + for i in range(frame_len): + win[i] = 0.54 - 0.46 * np.cos(a * i) + return win + +def get_window(wintype:Optional[str], winlen:int) -> np.ndarray: + """get window function + + Args: + wintype (Optional[str]): window type. + winlen (int): window length in samples. + + Raises: + ValueError: not support window. + + Returns: + np.ndarray: window coeffs. + """ + # calculate window + if not wintype or wintype == 'rectangular': + window = np.ones(winlen) + elif wintype == "hann": + window = hann_window(winlen) + elif wintype == "hamm": + window = hamm_window(winlen) + elif wintype == "povey": + window = povey_window(winlen) + else: + msg = f"{wintype} Not supported yet!" + raise ValueError(msg) + return window + + +def dft_matrix(n_fft:int, winlen:int=None, n_bin:int=None) -> Tuple[np.ndarray, np.ndarray, int]: + # https://en.wikipedia.org/wiki/Discrete_Fourier_transform + # (n_bins, n_fft) complex + if n_bin is None: + n_bin = 1 + n_fft // 2 + if winlen is None: + winlen = n_bin + # https://github.com/numpy/numpy/blob/v1.20.0/numpy/fft/_pocketfft.py#L49 + kernel_size = min(n_fft, winlen) + + n = np.arange(0, n_fft, 1.) + wsin = np.empty((n_bin, kernel_size)) #[Cout, kernel_size] + wcos = np.empty((n_bin, kernel_size)) #[Cout, kernel_size] + for k in range(n_bin): # Only half of the bins contain useful info + wsin[k,:] = -np.sin(2*np.pi*k*n/n_fft)[:kernel_size] + wcos[k,:] = np.cos(2*np.pi*k*n/n_fft)[:kernel_size] + w_real = wcos + w_imag = wsin + return w_real, w_imag, kernel_size + + +def dft_matrix_fast(n_fft:int, winlen:int=None, n_bin:int=None) -> Tuple[np.ndarray, np.ndarray, int]: + # (n_bins, n_fft) complex + if n_bin is None: + n_bin = 1 + n_fft // 2 + if winlen is None: + winlen = n_bin + # https://github.com/numpy/numpy/blob/v1.20.0/numpy/fft/_pocketfft.py#L49 + kernel_size = min(n_fft, winlen) + + # https://en.wikipedia.org/wiki/DFT_matrix + # https://ccrma.stanford.edu/~jos/st/Matrix_Formulation_DFT.html + weight = np.fft.fft(np.eye(n_fft))[:self.n_bin, :kernel_size] + w_real = weight.real + w_imag = weight.imag + return w_real, w_imag, kernel_size + + +def bin2hz(bin:Union[List[int], np.ndarray], N:int, sr:int)->List[float]: + """FFT bins to Hz. + + http://practicalcryptography.com/miscellaneous/machine-learning/intuitive-guide-discrete-fourier-transform/ + + Args: + bins (List[int] or np.ndarray): bin index. + N (int): the number of samples, or FFT points. + sr (int): sampling rate. + + Returns: + List[float]: Hz's. + """ + hz = bin * float(sr) / N + + +def hz2mel(hz): + """Convert a value in Hertz to Mels + + :param hz: a value in Hz. This can also be a numpy array, conversion proceeds element-wise. + :returns: a value in Mels. If an array was passed in, an identical sized array is returned. + """ + return 1127 * np.log(1+hz/700.0) + + +def mel2hz(mel): + """Convert a value in Mels to Hertz + + :param mel: a value in Mels. This can also be a numpy array, conversion proceeds element-wise. + :returns: a value in Hertz. If an array was passed in, an identical sized array is returned. + """ + return 700 * (np.exp(mel/1127.0)-1) + + + +def rms_to_db(rms: float): + """Root Mean Square to dB. + + Args: + rms ([float]): root mean square + + Returns: + float: dB + """ + return 20.0 * math.log10(max(1e-16, rms)) + + +def rms_to_dbfs(rms: float): + """Root Mean Square to dBFS. + https://fireattack.wordpress.com/2017/02/06/replaygain-loudness-normalization-and-applications/ + Audio is mix of sine wave, so 1 amp sine wave's Full scale is 0.7071, equal to -3.0103dB. + + dB = dBFS + 3.0103 + dBFS = db - 3.0103 + e.g. 0 dB = -3.0103 dBFS + + Args: + rms ([float]): root mean square + + Returns: + float: dBFS + """ + return rms_to_db(rms) - 3.0103 + + +def max_dbfs(sample_data: np.ndarray): + """Peak dBFS based on the maximum energy sample. + + Args: + sample_data ([np.ndarray]): float array, [-1, 1]. + + Returns: + float: dBFS + """ + # Peak dBFS based on the maximum energy sample. Will prevent overdrive if used for normalization. + return rms_to_dbfs(max(abs(np.min(sample_data)), abs(np.max(sample_data)))) + + +def mean_dbfs(sample_data): + """Peak dBFS based on the RMS energy. + + Args: + sample_data ([np.ndarray]): float array, [-1, 1]. + + Returns: + float: dBFS + """ + return rms_to_dbfs( + math.sqrt(np.mean(np.square(sample_data, dtype=np.float64)))) + + +def gain_db_to_ratio(gain_db: float): + """dB to ratio + + Args: + gain_db (float): gain in dB + + Returns: + float: scale in amp + """ + return math.pow(10.0, gain_db / 20.0) \ No newline at end of file diff --git a/third_party/paddle_audio/frontend/english.wav b/third_party/paddle_audio/frontend/english.wav new file mode 100644 index 0000000000000000000000000000000000000000..bb28291f69123209e6b7cc46b584d0a1f2c7bb16 GIT binary patch literal 35824 zcmW(-19%-<7d__;F5;$cYTIsV+qP}n#;aa!o2kuL+cwiSxVd*`4*u2u%hxtZlY3|O z*;sq+wWo2d>eau_A*4m6=GD6Q8=5Jc5JK@Wu{=H*M+hU~q-&?4oi5>L{(fI6ljx)@ zanxK@n;cUq$Q@E!`BhD_L@mL0SE}OVt4d8ukk{(B(nv6gz(>4FQenhH!tlHCN+(|8 z!XK`La*+^X;Dg~iLHLPC(w9^u|ESsOf!akDkYDNoergywNH&vURQ@hC#(vHNC8hF*kq%hvICE2Ams!?Q= zDyvG8HTcO2B!fy&8ObZzQ%Us=SN%xcQ4>f7QjR=VSJY8*OO;S7NqsUFcie>JBfVHj zZ7aE?dMQFW&~TDjja8rV8ih!IvRjQIQ%P~so=hYC$T@Y1d?sI1N7$q+`Hx&z$JBdu zO)XKu1JzHKP!r_AjI;tT_Nz#xH>YBWxmXeyPJFc#)>OsTEH93|Rr*&ycBfHj; zOrtJxLvs^aM0wSltitszz-gS9E!bx*MjaAA!?9|nMk8f%X*}n z+9}GbIb;pFsUFKysy`h>o0I8k5~)Snk{39Ik?IXRr67Jcqk4&ZyrqiL!D@rdODd?R zsw(-0Ynn}Z(m!e|{IWE;Kr54+RKjXUXa_o*bz&uyBXX!^sMy&H z*|@(p#6x>iliH*PeM_HeK{P+kcmq4h4$)8=LfQ}mSN}lPRpm(*_}M*~9p`;it%fJp zAaN?3%zzgy!#Q^)!|8t{kyNMm$!W5JoTpbwb$LP_kr8T*tfj`tukxMhjBD?ys*}8G zy-Y)P;Ws1Hag|j?lj$UjDy*KtkE)W~q^G(=O4FXC529#4sZA%7mujE-iOBh*zN?w4 z0KB*uu~ZSdQ)O0R_>FtjLJnw6&XT(%GfhV=Qis%`iSX@_qyRkV9`5Rjs)W4ol{6>Wkpaq+ zz3|Risv5o*N4CTA8X160a1F6@h2$Vli4WhsjjtWmj^sjCsEd=@qi*3Gmf*@K!{4*O z+g7NBD!-blE+K-_AoE;Q0Yo^(DYsC$V3%?@$7#6NLZl|?NZe!_@=Qxo2Dzj>?xQ6# z(pf}%E!go5IgjjAg>--q6ed4a5h8Jy8P!s?84*|x|1MQ2;i32CeYrrbLMD0)TlXMC z$t(3qWkw#FLu!+5Bs20=ZZZ=mm>1FVPQ4{3Nq5Bf3)xByB>&*Vi>c=7iJFOg`9lqZ z9}FQ2NhUHHmOibvA;UvgJ#19OG!a`l1x-(kj-8qVkhADX23?{)IH*& zb>L4oWH~aI-X=lHtNN0$G=xr4C**CFmM*4l>Vt1ABs=JCnuhKrzez1RpKhcKr`3qG zhKKzj-|*QywM8D5YmqrMvOxWd_cGw+$;o>}{a{22p~J~DrN1hqTB<6zhbHQ(3WZ<(hSwxRmD;OPk!++D8H9XTT$Lgh$yMC*LAe25a#d|n zHPm^zPUcb{RAn+r6_k6#MVWzgMkHmFRbhw0G#9JP#<31;Dx1yrvszjweVg7uOGEeL zYQt$SdYGKSwa&$TtVLG2i7JAM09*OVYP=$Yx85a1=yLdJHh95K)T5SUGCb~pDvCQ! zf&6e4`D(55!D`b~Fg|+Vw4b34)KZbC3HenkbxPI2dDkMbY7_2tF+9vb+=Rp9%c*xV zQAVlZi1S^zqMV4rp=ukV_l(Mb6Bw>;$z;-y$<|?$QKzAkN4C$Acl?PCSFw|n^4bks(y&*vWUSJ>Vdo| zlc{F-sICUUGBfb|v9L%bc#?&PTddl^GkPG7ub`SQ#rrImd*yGL2e0DBPrjDD(9i=a-%iZJn*4-rFUN_b5X{>eV&ER;j&2Ub^K@5=5?6|~b+)!o}Q z$}DAuo8z>7G>_WE2icbb7XyW@L)LG5k28j!7G8K>B|41V(x@@pXlm5gdueyrK=?v- z)E*9hY$qy+I6hH?%YWo^`CP50f^A^)5nKD%ac!}7pLu9g)QwfDq56f`T_*?2`S{sL zoa`A~Z4{}8>|UFsBe!wis}XTWP$5R+gk37VidKVFN5pw2#GZ{RxEyhn7Bw>^^2K!GCAML*NP}PUYv9?InA9j{DbI$-LX_fl zIOjXWcTq@|L1k~H7O98owTgueYmgh{HOWby%FQCRILP1d3Zjj8#Q!)Iof*zZ?iU;6 z0MX64Vh`qNRXrNd?rJsNM?#h+yB*Obq_g{=@ki^gRnzv9hhnE~`7b7(j_2_$d^E7x z&M!8qU9^^VU(amxGp3l8-PJrt+*izD`dhk5T^C#V8K<;U)VX6|IV-kx>QH zhm7jO^qKkqy}n*v|EW#Xk`O-y=^M3Lt`^xugg7g!)2}R=r2^LMEsLP4D{=yn^G799 zZRKOJSIn17fMN^K^E8HjMh9>MQJ|6a$R|EESEj{zmf}l!3sF_Bk?B=UctvG(S+zj+ znTIZ9J3OVlEGD}0=T1{6r8C+tY@f1oJFT66L1E_%?U6 z9?f#I7j&NNX%F@{_HB*t_^0yk)N%gA&sG7IT)SvIFzdSNy3)8xxC?m81by(Xa%VQC z)3maYQ^We{AL4)MzZ@uLb+8;u*%x>YwS$({7n*I{_uUsQQ^jlKF8wu8)<5h5)= z1*<)=z9ZVM@_{NVYpbU-D(ggFrDrt18He?(EWPrJZsHSH&RDygUC#N+SIDudDfvi_ zQ^~q%Q?zI;2a6>w@pCqB#=kproz2d7XAok4xa^JY<1Hyk$5Ic?L_8{5?i6ouUS;g` z*1CWdSYk!mv+NR1L9tGqV@c+*unwt@r>~i2MEE#Yc5MhPM;a@ab1tb-;*|J^xSzk$ z{Cf1Kg3s+pvdkFn>E&JKt?TXRF}yW{E(h)L_H=dE7N}Rwg1~kkNj#izJdyc-1v=Oh z9ZhJe9X0i>rsi4Zec(;>G;rrNm*_tBk(89N&KkR={msg3)d`dj^tT2%UF9}XpJmev zxfXbbggC*=JwuHmbiIlZMZ{F6p0y>g)9(tr3GA{?Ii+PKnp-br&T^%2?K0BqhiQ8< zL=BZS#Xcvgsd2d;Qq+w08j@k;AiDp!%MFL;Lv+&%`b!%MUm|yvK2QJ$W zogHF>N}!h6AiR9q2ALM6TO1MYT0zT@V{*Q@WvxwI6aO%-TwJC<|NhGIXJld}`#h;- z+TQ*l2SXZ%912SBW!@a#Ri3J@oZ41-(2n*uN~{*YI4)QGw#1@=Ce9(TTs5X`wNb`8 zS3i%>yDBI$XrHHwi)(@uk#SC6d!7|(tqw#7dIl<4%j``2n#fD0vIF`K^Q3FPySz8n z+aRc;dy~F^6co3e4VFLf%U|4|F|apK+WzEp6~Jojs=m|=fuH6vJLq-UWtB(9a^bXa z=GmXE`&MH+uhWiS5N%X|G-k=QfvAAGnl4`R*1Vxp(Kf9?0hbkG4|2NkouZ_wL0_1I z!n>p$oas}#LdgcZYLE_MxwFrH?f;W-GQMN{r1(Hw>p#uoUL;a`1es)J4;~gaIs9u_ z!;n_qy6#}tNb|Bjm7W!u?A-pAi525Z#X8Tk1yeCWA`Vs%+|rB!aF3(lx%R=%%IIK z)7Yg|Wf_QD+_CJWk%^b%_rxuZ%bUdF)kMFEQJ{H}OZ} z3f~%k%RoA-vt7^$6_e#^(nPCda`!NAHt%wGeq%JfDbCxqtoMOlR$hCX9dG;WfldiN zNJOfiWCeSs71G;kXX#v(PgHTwoUL)Db$ zgVM7$d?Bo{A*n>-vG@k@v*UlnABvAoDCJucXe+j}vhE5Yx5G~*b0R*69`|lCUuuk% zBT{5{S_H29J|@KDH_|2CNR08XaKI*L5w0bkltGPy?t1Tf{sCUeYR=YMgReNonf=5+ zC@IDl>-(HEEpXj#D!!;fELdM*gt+#(yLf*D9S;r;KIj=|v?e{A82>NdKHs*a6aI~X zZPsEtH_s=BlOL?S(ZsdZ{llHX9WbIa56L8`bKB}+Ewa+u%j^WZmh-?l%2UhNY7I@J zebe@7OW9-+EDt*2_DpN4b-{XMEwDaWbL|ZLmN-E=nSVkXrp}+CP5R;~3kAPp%XvNP zus>(gfW+wu6%x89v`NU8upyzFZ+KuV&q%YHB&c~*86IUhOv<8aG^pBpyJ<>ZqsAEtq?UeM9tfx$n#cU&vANY&9U;Lqdxl(@n7CaFiD zx%J(i!57Nf^sUy%tmp3HsqG=|SR+XLrD_Uc_pwd|)&@QYCR!!!2ljp^miLr}$!8j= z)zhA`DRj14Chj@|QOOoq^Q_%g1v{fNkFQomjg}!_Q;bZPCtbr7n}g45TSP6pGIIBS ze#5`dH!*Q&LQH(k1Uqq~Kh8cRzp_5A*q~iusUw2IMet{LSK~eFKzow{s+)LX=L~G~ z#Uw-|Tu4~$%Nz)Gvd9{=qn^=~!IQ=t>q&G!b^S5ZnGKBM`T|yutQDW^c>&_@o^&SZ ztADap!bu^}Ez$nkX``_Fz2~*}ruUX-iF>NqN-s&Di(R(k|Cuy6DW5-eppNy)ZpCZJ zbR>+W(a#yeEasYL<~15={pnJ5SCru|?Kf6q#K11AD{|yv9w8e5bG$-7zloeg|2cw@$mdoXSQot zSlX18(^OA&H9W*UO1-oC`|tV}TZJ6kId2yVFyFlRa`B%N^7tz|S=0`e%xvk&8oVuJ zamc%%AMX1`ZfzP326y)X?86RG$zG6jJK=2n+Jp>A-K^al%pfhQ_0`)OSDJ5gdbg_ zHwUz(q(&PGIe29hFn6@Lr= z^FTpoFyH62u6JGfCS&KzFI#+jG_q>~fn}oCoj`DVPE!S7FO>`1HMu*Zk z^^K>uR{4r09E~rK_%*4t-CKY;q_OOP)=JN3C}Xt=)anXy9WbNJRK`H9BpogHwa98ow$?a&)PUjK#fL|*N zF8d!coKyyXFqZ5j=g2ZrhzwKvB+ zbY}BkqBLo!M|yjQ=Sndk6h@O*CvRlk`bi0z(+Nx|naOR5M;LY08g|sIqI{0K;9T-rVvblLzJhyLhQ9k4_|^XSmQuqWZ-oS#Rn<~^ ze;}3h&Hllcin6?dbKz)z6uE+%8Wk&lUF@*GjXE(OD1CbLpn;(teVQB8#1r6zn^Z_{Ucy@Y*gej>`^U zpf`hWt4@>HX}!AH%yrLo$(7p`W$x79vM(e6Mvd@BP95i(bBR9^A?iPs7wkC&gOd+k zd2Q;WJ-~yMWJBqAvKWlacyJJJWDB)jJqN2Z9?aiGxm+%mM`eP{2F8M^gR-)GBo=_J zdM=uygV(@k9mV&Ki8?|<_m)Cr25a-3f8|-lLoryrXIb3Ap*xe6NNz=356)udBZc`= zJIFc0`+?Q3A)@R)Ny`&<#NSQ$C+VA&Ph=!pSd8A#9PE-V$2Hb<+`Mn3HLmHObxUuh zuV$B25q>W)G3k)6Yf{cY7rPhuom#Rsu3!UZJHzP+t(LLH?CM(ND(V_wRyS^IX<1kD zL9Q1M`3!!L7ZESSBycTTWG!T~Qs7;el7?)R#`HN_W#rcR%7e+qI`IBWkQrWrpGzyZ zi5QVb-U17^Qtbz;T~mdlCJdEPh@zv2iYTyGZ@_vzkWJ(|K|~DR=MT>>%EKbp`FWmR zRFTJttv~cm3il-c6PYIbm#4f|NKSHUI2-vOISO-@3(jGGe8RoB$oT1rnFCLqbZS36 zp_MVJm_N*H=n9{kE6mPjO|zX@$6RO})(#Sfx3``p1$;Y`+6UU%jKAj&KftGRm*_25 zkSMJOs{R$%6IWZ;TyvyRK)=A6)8r&XRh7F%VR4@S3wwOw14RnVcc!W=V8#2B5oh!>@0mqilU%!9n_T5wk>)*Z8TrgN zSeyKv{3FpZ4R@q-#NqY;dnLMHL*66Hv@^y>GoR~`S<8$uG8qH)<=Q69$gX1&`c@th zllU&Dt@F%z!}Ea4uSoV{2K0e^CXuuwy+k8fKK6<>qKnB3FufGhlpC1CypWDef$3E$ zWQB+F9K7j?l+vZ{$UgF?7!5ABpLh%|Y@3`WzaevW;*q>OKf+6jx#FJ422Y)?_R!a%ni1T3FN-`A*L&nHm8qNSAN;QC1F_Hm$;3IvjP#kjvR>DTSGdZz0iF6 zHDkQl+*RHE)BVRC?OqHk^atkOZEvx{?TSt=PWUUkomD$f0ciJ@Gec&g`LvOGoF1%R zy@m@V}n^WrR`wwj71U$uSi%%SHZPDktx+2TseQYa<=A@Q^@J5A1uNmf~3 zr+9aK(ZrAbS6t^sS zB=7KcP8p{+? zKT-9cfwLVVeu`?C;{1j^LPSXcSrT3{6WO7oD1vHpnS7x0wHUKTkP~_$YS7$xTC&#Q^$%l0xI{e_qn&)#D}SEAL2DFv-6}Yl@gh!qm&M5@^n(5LHt>c?Mh@ek z9;Nrx4uBJnr?uHZR*&0*sqWuw1?DB~@~utU5*T8)1{?ar z=C&6dYeV?t9V)d!W{|t0J2R{hte?^f=nM7c`V`H>KB@kqtJBZ^$8P3C@}hjX^U=Oz zk98*Uf5inEPTFGL-dlU3z1D7Pf3#wHC%vEER1aw7wDar`JH|NeK&FDd?#mDHsiFd| zznz%MH6Fz|pM#9nm~N-TX-OJI>(ghnF6+-q16>rOwJ^`#fz!_*ZE;of7Ei!_6CTW) z^Gm!eaMyD1+g|a2U*o3e56_s4dBQAkzc1BbvWa~(GX|Fp%M{wt(^s#?3ThjS+^+g& zU9G2T&okQ-16%y@{%lqbyQ7`czHHsF>exM;*LLBAqqBw&pZE(%0X@}UYZI|=Sv3Yix19KZ6CNw%t$*OAg^)I?%9M`Mp zx3tb$er=&PORr~K)!VaZW%C&RPJBf!>@5C;ZJXK^oaTJ2m?itFY@`!e1kSb{O~a0` z5||r&X2I+MW`%jdivOLpjsy;`4+Oo*iFfkye8`Wfc_yBc>pYCN<4^b$QBr=BE!9MI zNTsApbltlw^mXVGRp&NIMQy2~h*Z zDrc}AY0n3HumG9(IG>Iz7=V&yHfct~*a22p%b;Cm$-xTf^aW~e4Jb_J%4Nb2L&yyhfxoLSnu zg?kAx-stD_XuY;((P89|+$t=7gYSi>?sPucUBDeavgbOjcs}t^n;3t`xO-Imu&_LXQ)}n}-CxsY*ocSLg!|UTUPJ!p_gI9{;6ZlL%j$h_O z#WMMbDDAc@t2fqN-<+zaHQ%`Ncwc%Bx>g!OJH!^!wW_8}iHUl8vEQj|UkI!Q<6gvC zWiNLQIRjx~iJCG98TbU<#`gFA>J02njrcBji)1faLUzLE9BKqMm6({Il`P{bkViFWMT>r=U1HPj&gLK4Ib}8#fpJiY#U~;U!X4-kJ(5D=$`gLxmAvQRNXK`Y=#Ny zZy>B*=+VFNJR&c+|M9#SGV5_3DN=~F;;5>qoicm4F6+zKTNYvVbZ>HZcAYlT8>{sg ztpeK&_GK}qv7OZ-@zNP%7qH*jJDrsLvg37Pz?!;|hyD`_fCdLb2a}eTN9TBfz1NoO zhxIndhdZ_LT3gI&m^Oz!gH|aGX@UNrw=5+8h;<^n*vOl}W456a9D&TAPL9Tudx82x z3TV^x)cPm3mOavv8I6sadIGD8DfBLS30jq4=mRdw_HvU*xq7eu7o zjTvZPsK@F-SJp=L#=k9~UYP^ELQeHa?m)jX3W%bu%n0S(Y)p`sh|A&#=8IINr}wp& z#w+07V762v=2o+^nZp>T*VNBwYguhNn+W2jqewKSDAV{mCkx*XJbas%;YZ4H_9*agS-lrz&}8d>*Z;D9i%d#hg~K+0DZ29Mb_f8Y>1wp zBAF~N8{v1y$(-ai>t|d7|96GdmqLrnjq1MSRIOAFJ0n{4n(YXc7JaPnjkH6N~0aagR_}X##3^;tB+z&f1Ku0!Q zPQrA3nrtAm!aJhCT(p!8(Wk7H{mDyK$XID~(o*5-efoB@5vH1HbdMgY7uRmn8RR&r zOT($H3d>KtCLe|PEC4oWD6DZ9QGbXR2ENRSeA);)p$TA3Mv$d+6uZOru!(E|Iw-G} zN6VrGX%|^zc99l_T4^lKHaXc31;vDSqQZ$jFP+}5#o1UdJp#L-2 zEtXG<1U@gvPS76U3gYzo~tC36!#=m-w-dOWp&tU zKD^{N?6C>GS#Nm2LVUJJHj~+98g$F`fYH8*g813>vb0Jf?X(NVdZQ(|2b2TKcd7fQ64^-7CM?w(BzgQ z7Kvn^SY~*`G&YKDVZT`@vc@@9g#8OG))J`4pmBv}HW&DW0`f9650``+I*g3)TR#x; zzli#dWU97VZ^vw?Ny;$8=xCJH9zvJ%hAf~-P&{ozmwrz5f#xSjJ`}^yfsF-riiKxv z1InDJazVKhz*lKukJa#%kJ5`wGa25H3_6Q$a*sTQiF;YOMNCATH-|-MLwmYZ#>qZp zJ_|DTn6->y+IE&xuZ*rUx4B-Q0Usa4rqX=S%p|B2s8YYg9)1yZIuWzks!%MUqjvs) zY0d}5#Sm!b!k|%nMUvwLqu58*LTir9F$Zd{Vyqr(#!7&b{06<=CN&86SW4N5_3ko< z{0fHYHBM(Hddqvz2ZcaA`I6)WcXgjVry1y0+Fq-tpVq?JE;58%B?<71X+Ty#;6?Xj zPnjEcy#?pG1$}xb^qT+T^KkgzSEWH?+YvtV41CU4baUn8ckvclj48P9dWh2B=+A;= zEPOL76fR?AUqsnz=swHRky?yV+x5o?Xzy5*-rOt&hU&C_hRvcEf%c9->7%LE@{(Yp zGT(tI?MvqhKac)quSnwe_#u%8s8W&CP-nVm0Xl^;c*c7~{5tdv&*(^+gLG zjh0a{P}n?{4EOXKy1vQ4Efdivlz}?yq&g2~ijdu86MmZZg?uZ@iiPMCf`MyCLC;b_mPO?XLoanpzEnZdH%O==3s>#)!LGBnchTQs~ffLE}3E$~_&-Y#Y?YwdyQ1cl}TcBf!m-hHkAG z>!MZ0l*fa7(T1G|W_4&QG6gEULx_$ussS`d$H2TrL1*H?I^9$b==wb1h5E}CvY~pV zZbJuINiG&Wz;HGN^LG?>_?tf)qWV0NpTT@YAtN6~AJ?lXh3T zV{~-2Hpl8e*i3Dgp3zvUm(#-OP*nhI*&k66ezslI6O~|@a;P!;#00q){aYLI1a)P; zV$d!&K@LoYp6dkEPv4+H-$}N>&)bpeq$Jd5?V!g@BAIA5>c{^yfwrtTG`9V~Rvks< zJ%+k?0O$Q#{zCpsM895#9A^J&gTN-IhL*`@!?gKq7;@TaXqvO2cKX0yHAQ}mM7L1@ zn%gh(w3>@Kz+kch`pA^3HZ)3qnb1D+0+@w`LX!oNZ?mGFOh?sSiL9{$aaRPq;X5!s zhw$Df|IZO+#P=FN4fupCWkqyp)Nu7SH|jrGCC$*c=zaAQE202@rh>s=R% zP-Xj~7kmxX&=EC=%pw`cY9O8g$T#OmZs=hLLdDq#s=CzB2YrCPIRV<|m6%}cBhT^w z?V$;8K!3nSZ&0UqBR`(U&q`pVeBk7ogQ?02D-S?+zRcdSyU+l(rybY=ZJsucT_;_D zJ5s>UVt}(=;LaAwPcjuWvz_4;jqo1hp~S4Mge-!ecny`?ZRy}@2LkIAM|}3g9h3x1 zxEq!I8R8@(w5BES)BmAv4@3m-K=j;2%{&B4FHwo)9~PuvGeo_es%1-20(@8}bk6}1FZP3JuK-_Z1RkskrW1cP+e=gh;KJ4Q#vA>rKkJx-WR}XwUq=h z84bn!DzM;raIzUtb%T%_i=i&Xz+28p4|Waqp>h>cu^_`MG z(XD=#TY(Jvp@s}lGiY_Ky8b{L&U(pgBpdb z8-i#l3{AF9=aRMRKlofz%n#n9IzN^}P+6t{6Q9I6DmfSaIR;fD2Jv@A76Ou*hYYq8 z9nJ}ngbIEaYPg(GWv9T5<1@*}UTSU4(XK|Wjm8jtCy;j~t(hiS932eB;4@(D>0%yu zjV`DbMddy)t)s*fu@n4SZ7BL@gYB&jwl*ECpAD*eFr3I3MWN=+iw-LoQ?mE;6J1a1 zpeug_tNC&Lcc88R3DxE~)WF_|&lu$WcgWlS$=1kMJLv;V6|b{z^grO&2lQWBknV!c z`w}v3Ms)_t@GsCY&xTsx#+JZ#Xg;q%iCYWr$*?Dq8Cvk!sIRqAjoYCLQ!tNj(CthF z!p*NL!n=xrO^iZTCx{>muXYU8x;}7ADd5`kvH&vwOLdqm$LVx3j+x0^?~S>7H_eai zQUo}>AZq~S{A*x|rudz~zyrPTK6y~#uM6`3N|qTIKA1d4?MMdJ_AT~A;=p@e#O}me zG7!58ZSk)cD&ZsaC;dhb(q`DMI6#)dS`$bQWZVa+rw(Fy3A|<}t~ySpL9|n-k9xAt zEFJ3!X0S4CNm~QEo=3+10Hl@yXD|#?q_W5{8uoPV0M}dsdR_vb`2r2-aYXcE;P@G^ z!zRp@&ZEoOEUzP~BG4tK!?sBsRUcaPDzNZ3eE$@F_9z(Pq!h91Lal zTX^huPfZ9`o*Qf=}Jp;~JhKltH8x);MZ$xJ{OyhPVBjiC<&k6;5b>!=__$m{6r(@{* zM*&m+!W3pKG{argdgSwKxQ-Fv+dHBvt;R%ZJS;L0{9D+JqL1TIymB^_{06;T`0;%w3ZRTaW}4ui`1Hu8^#%mb1g9c?P$ znR-BpGszBgRy&ZP7eUoO4qG~7;Z0M>e0*mYcDY(%_JK_t%q9n7vuY+7>_|+Ymf$Wa z5YJKgPcSO#QXq{ZuviCd3DrQpJ%!V$ij3YGHrR#jlX=)r=z}QS0^d4^EsjZ0ZJ)(w zTaj%#pq4hpYt(@?t02aTLC@V9pH~H%I*IwyX;P0KqC06HS^?Lx0J(AyEkiSb%Pj(S zYAG}dln#QI{jg#kLZDx;NHFf%PvG@G{Wgz#CsgazwPib9PfAnle+%!m7%Cs>kwz7a6VJu zM@MkK+u-#RkfRpix5uK&H-Ufrz3OT>#R~BIhOpyg+}nL*^aUg(T}F4%?yyV+)R8gx z(~PEqZ9K3|cf?dUoes;?1CL(|SS*~}Mc&?`{_ZRdzz&HI6H*EvDThk;_kCL9YMbJ- zR=^h>@l&VJIWm$KaaDxmCwUM%6=C6)Kw%59ZMg=vSptt3Lb?FKH$$AYgExh4EI9NA?DEb$n5@(B2%VW_SLaVKwJuMM!;EL__@SYSD9HC6qC ze7X|Wm4*Mg1rr=fQlTl<4@KA5$pl?h{dAd>yC#NU|?lT}ckT)6VTnDTGE&H;a(3+)8OZ{7y>FO2@%!B1p> z@9bAiNGN?o8?(D?wpLBwpu6;xT0HLP0xHN~KA|e?atfT)EHGTzWnEyf1HgtYz(3{! zZb*fEYO5feW-;Uqj(GVFzevE25k=J0MpoE`GfluAUMzW!^Nm4IxdOSP26&JY@c+E( z2VQ*>I_twgF%h`ZdB_jDaj!X1_Xkm(RYZ;SVsl^*Z2`|+j7;N4wb_hQt%D37hFz9) z@P_Hw`pJQs&=!$c89S3f<_lHae<<@QbpDseOpbcjqB zeTT35z>eS07f(he%8I%3YvA*P;DxSAzbuEm_z+uAEfC#4G60==K9(OdimP-stw9rz zo&WA;v_f9bh>ya+K2_mk-4L&j(1mBid6j{6Y9n7~f)5>r725pYKH;cgLFg(AAiGXN z<#>*H8^_ALVR5s%K*fu6)Kec*obz*dc5_tdCcAHV}8(GJ=keU-7@_@%$s zwxO48K#Hj1@)cBS23U*dvgqIaR`n8=8iw=Eh|K3h&A(16fQPC<-y=s{M0Wp)`;5gb zX9n2q%BV6I;5V`4B~ac@MA~I+yZzNcTtn~O1f0nnFqfCXBh~_^+#YB94s6UE}? zPlw!e7;)DCwKE9)>MwXs7^={i|96@cu~-y&uP9K!b;MUW{KN(L@d9X*j^NaSpgm}V z+%+AZ(h`53s{`aTEw0Vemm3$2YQ}KAqLv8U6{Mbv-@FB;Y$50{Zpk5FL_Q$5YyxgL zh|?H{T%!Zwh9ZYgL|%zO#{C!gdkuK7R$$RQs3^lodUOl-=tJCNYsz7rmWaz(*zgNp z>lfJhE*y z)RhXbbXIuNG-SYJsDY#K9>>tbj3d#A|DEug5On8HV9{WBUprXzHc&+bI2;{W;Wv8q zOSq%az<;6Y6_DCac=TTM-5-zz^8p>s#ubjoJKo0TcYUDsQ`%>JoiWBpX;jj4X^&_N zvQ-Wi(fp3{!g zF*}AjRs?LE9jrS+MhU_;eL?)2$IH@!&}Z;m-c<$uEYb z3LvwOQdtQ{_jw9Bi(oUOdCUklu4{GK60%OF7rPxF^z~hxhJ1!Nji()q2E#cKSh^7= zTjkLU4h53<%QJQZL-7&ObqSa@oQwoZ+X;+XJ}_$4&>zN-!89FP0JsG&3`-r#>VA$KRP&S)BK&ESqezloQ4i+~8 zzkg3QLAGuHe0>E764N$d?wZJo5_{yw;LSTxujirWSHza0gTDU_7``_^_F0hG76U0g z17dg%J>)dxABC!a1s`|8q$}`3*g-)b60hcg6$oX;w2pdbFKG77aOaq;5Z*+@6u+nO99Y27TGQ$(vp#CM2 zKeV9MK@T_f>tT9s?13!BCSn07TPEA-?b-HP=ObSwvyg6>@Q=~B_Fb!}<-uv_?K0pQPPrf5hA3@yzz@#fBc#)rK z1zoQ#GTg52u3D~RW?kc&c9GUrX~ZAGtPCwxWk93f(Vqsu=+CYo>4?4d( z@XZKd@B%EoHj!9=##HJn#~jv=Gu7XHVMvXiXi zuuEz2NOmT_XdA7H{;ysN8|Sf@&$myWkAqMIIsioIEE6+wvfOr#PUi>X_BcV~AQ_axVSqnv(}{*aaUTWek* zM_@!?k(JkZ&QHo-K$69^tJ+!2T7y|A4aT-=dEn*xK;B{Cu!_-vv_E}?8Vd#u8q`jr z7MMa?yq5|yi$Scub_Kegq1tWsk{-wGY^975xuJW`3~sI{u4IGigKWQt-oiX~AMFez znFU@M4rZByPyG-5SY71Qf$BZbYF6wj=SDQ%z@%S5HChNr;0BnV1(?W8#^kq{+Jy}D z9dql)A`EQkWtkQ|!YpKf^@ybrh^V|kM_&*>Pw~#Nz!*vJq8H#U&Y@;k2bQ~tUG;G+ zmA)CeHtJsLdSQ+>7HKKibG1dxat>RU1KR>wtRD7PXM{+hhLUUaI4h$Kz*cT!wi%ss zHz30lU^QA`k~|XVrY-vRgRoX%Iax%D>S7!2`Lleeu9NTd2b-ro*7j&wwH~ZLIJ7Bh z6PVWN;1hdcf?gh1RSwtv7;3kMq!FstVN9H|>mRj&+EvW1B%Xe71=`7%&U>c}6t2_3 zyG>Ag$qhP(T||Y8V2u!G`+;TJgV{QZbH5KgK!4zXkGSWnsHFoiRoH;&E{RHY8C*eG z^o@*kRVSdJTZ9eCC@{+n<#9Ea&eX!q*6u=H?(OVN;dyT+Y1L>pnaNocxRjJVXH1n(OGqcV^U}`eWaPtz zkef*fwnMwDmxn5Hi9TG*$!_5JCYxk}*bEih9Z?MNUH}zy+5fxc)X1fKup3P=#~a4( zv2kLHO~tT)m-fGTJ|13R}MnGTc=&p8yaPd9pJ;ev7Gd;*5f1k1_vkB4Yn8v z%=J`u#=d$JOmSMsI%GMUp*JmOH@9tOsi}8g&5}NtPEq z`EP6ow}eWps|saaW4wDzP|uLlA>)Fpc}u${YhP6v{wz=DeZNESJlcmIMeg-P+G9ooJa7nPoDJ@Jtq3K!) zR-LxSQ=Z-dH}Apu&O_B`huks_m~A8wNIlHtv&lio42{)S)g7qnHJ-DQ8qYBrFP~s@ zx;QorW~19&i!5;s2yZ`XLuY8Tm>AAOxzDKrwf#F0s#cTf?4n-Q8xYTD*zvz(D% zyH28Ejgqzq)Cm+0@W4R34sS2>;aM#m*cnW2&q24+ne_mgy~Q#?AJ|VT3ysq;JjZC7 ze1_ezRydkr$XbC3Y_=!+yWOpp^jTT@0jm3w1vNsI~$AgyZR1<$+zg|G)2j zkJt$YOVJD&>@BkB2Vk0Gu!AOwi37+r;qbt1z_EkT=M0sHpeC(`-QYre1<#AEgBZD2 z6{A_SHO3tG$)G^U&Cs?XKfUW+$@H~km~h)o{55?lakH@cbAbGD4`DgL%s{D6VgaFVIJB#%!|(xT&)0 z4y+@@E>yqt(Av8&SzQSq>;~)~h52a%I!!xmp7Ev%{Sa0xY;!R4NWCKcB!=4y{Zo7% z-+SK-|8lF7(+3+cU!hFBjb}~mlM~^;uk|3;J$EkeDQ{cvO7|6G0L!5I@YA+w$6J~0 zC3bme;AV0^4;9a#=Kaohi9~su*esQC*}UlTxoW$@&B=NUCRZiluZ1xgtc|;Qik{`N zOaUF+cm5Byt8$CHatgN5YOo^OQsku=c8`7n(wPAaR0}A4F}R@r!1@+YD=;|H4-tPs;e}e{L;=B6t{&b-FkOp}ZU}DA~hm8)sanJqg|xL5scn z-Rq1pEQh+tyF1~~A0Kw!II+$ZXB?g!5wLICx9wrhEAGb=N;2tP%onbP?l{*&v$1g! z^|?5-;$wk_IA#Fr;pn z^Bc_PA9#9Nq`Q5AVDbhjrv@0_*H1G|Jxpbh-R zZ>p2GvUQL?}g=bA@sdH*S@pOuNq@bEAemQ6C@78wXM8M+KL=&Zk21CXzPdnmj7np zh_%&D=6ti;+Ec7xYk(E%Xfg|p)JH%``P(zwo8P<9UEDmZ&7f}LLX=j<>#v6*bh=Z) zu5En@Bn0YPm2jHH<#kd_a~pq*!X`I{7>D#YptLm5s%+7umPW6lg|MQesI0~h+aIl? zR&#r%^B8;Ro7D+C5pFH;!YF*~0VCI+ErS-nH#>}J8lLzk^NH@zNXI%wuun4-dq?xo z*I$8lr!5p3q2L>$ai>*q-^ryH_He**ycIgC^ii1`8RQHpqFph225k<#88#<$L{LHZ zbiEKYu0VO3h9`r%cC!_2d96OyYI_8*M?qRiyRGjwJmv$VJ~mAA>RYt*KuSe` zni^>bproR7KVo2|lfhYGudqvFV?P~sZfBD6m_(gMRZW9gY!jN94M101i0uJ3uSx2v z@8SY4&C~EU$T)?dt1J(d!cr(CtH1&&F;9PqP0iY12Im4_6_?@AkrKp30`!I#kgtoA zC-k6R!1GVYvalv$TY_JCwwlGY8`vq^YURXE_k@7g%57^-E<71xj#b=xZuN5h72A}b zPSrgx!Fk9SrI%@XpBj4(*z%KGsUJ+W0 z1lbJJB?2YWMP#STsJj#RRbCDCvMhS`w3xMQ25Q)Zn)sEjG_x%B*OfGzGG(oYsHV zU1u)(<=UP}zN``5B1ZTRc*kp%l|y0;SFzWbz0HecAkE+pYs~6&<~XAYqgN5_yl!WJ z?ebDtsP*;C@D%WLVjND&O{Hp7^ykFA%!gyr4%7gh+!Q;;YH5toeR@F8gPP-=Gn1-j z0og%L^hRCqDR@sNe2!k+8#2T)!V*x5z4lf#Y|&_1@}k$?Z~bWx0>z1PL&X1d)H+kF zY*u9FwQrBy$U#w)qH;v_i;VYo@m!GKyKAk9di!9d)CQ@o zgKDU4C@FMA_wc?&tv=Q>tG442EU|<#NSp8d#~0~8?Thg}@F?1N`M$8-Y0K_(-fUr3 zHtU#Qp+N{3Nko$vvo?JBch)Xvi_lN9;VTSdXUOcCuWbWM+l|I`rhH2JO^OBWj~7lk zN6w$MWa#^T>Uw&su`pON`g0UWOCHT=vvX| zqgq7H_I;&AO9Py;=F(8P;I7mz!EB)>p>3fOM1{||E|o@G5yQ?i!QCeoQr0pGn|yiw zIs6H}2HrB-YPpIy3)E`}YT{O?T0Ai8`>>9V7(=B{5nBlTDItJ6>j zjMYqayc(xwz<73*QbCz0pJSbO5=^Hs%9gT5Q$4rdSTASPVl_2#4zqf*v5Tg_-a0M~ zg-Lz^tU9x}oSZVj1>?0pp-9+d{$h?Wmzf{UA1yzbvUL8h*!uZ zy*QX>VY{8Hi-(n#9=|WiSKXh-Kiuc_-cx77;A%j%xRAA+jKKvZGit4Va5f5Sg+0YQXK=bHEG>4^+b{IS!Az;^ z0&h}#hb|i}>;~>NVUM^(Dj^q_pQ6xwicWJT%GJ8kLAo24+@8)=YnHjy_{#7brHrSB zZMvwy-;xJBrJCAmU14^9LC&_B-2bZGkD9up^htT(%@p-bOu5*_F-M}yMGf=U_e__I z3$^U=#@^6nJT*pz3WVka>j$R?!wA*g&`AA)ansV={$eTRv*vi8`CQ*sUt8ZGZ<6*- z$qkSEHr*wYT7Q8(&uV1)jUxKK&~KqZp+7@g^cH4qJFojd7!G4{yHZeXgf8<}wYMs( zJMqtYpsn|O=UIs_iXzEEPL#sa^=_dj!9JlZ`dOo^b--Tgl!D2WNxUakr;eXbUUghP z0rN9K8bbwk#A#|5w|d~=@EKj!2lEX1^7J%LRAc*p46^@V_ zan)@uyV|&js?pP8md3czBcqo18!-wm-SO5a{YbDwumXNtk3+peYA7aDKGZ0*Ce)GZ zxoXj~kh-evQQyq;XYt4Tx_NVZepde>@9n`VO^?gy1IE6sHPpBtdKp|Byc%p7YOAB( zvvZ<1{a%_Vmr&{}6_ql|R3)giRePz!)SK!YG~h)&ca?Hd83FHGYEdh+DRe$`Kp$tw zXc;qtD6MlwyK99@;%TyqpgdC9gVN-MG!1!Y0^Vww@JO0&-ZisZ<5>y4c)hpAS~bxC zmc})!D(m48nRJBR*$&y`sTDd<&y^HAD1Ck9qesNXqzlE=j_w$7%#%YIE&OU{F^A%! z(MKO-Ofs73%Yw~PF9co%wx?zfl`}kcZo!rgtNFa!dh(UfGVq z^oH=tNwIDk2SeF|qf>pszM&{%l6l>VaBR9I8NoQZinXbcyC}O=)syV`;3?#JrEXIt zd|MvKy~W4wGE~(Yt@36k{I$NA z(qYTX7-xKIwnkSy1#aj_ zGW}xiQ{je`Nx7)J2ebY|&MB=D;u(X6)SKyM{2#rER4@+w%t?Pd)p&x<>cZxJa2bX zsc!znm*$<(z&M19Q!XQqSo7#7* zHOvBj>5YO92mVtV-6!5W0OWP7U&t7EK{3qImIClLqTS%P2ctq$^k9~}`B9ks+~ zYGc#_YBX5sW4V(w5!CFa`AE+YTA#Wha4Rq~cvfF;ma{gRZwmoJ<>pch+V|{u*g2@zmys(>HM5xqs6vXN4o$U+g9W`6ijggL zhb24;Cjdq3jW@*ao(>G;+ZTDc)S zm|EaWbcD=fx0#blyED(_7M=|?sORE@n<(6()S+s6btWUJsf*Po>MPYz66Euu4>$9C zP7-^6Z~HzG|CTuo7tNREHtUFe(CH1jyn)`|A$a6V@Eson;{%UH`cA9q3qOFZ@s+zR zP4#ewbyX1;)<1B3>0%W_tyqstvN&qyRd#-|$wHD(tL*FMTjeRQRz)xTJ1%asy?=Xd zsDtHJ;sQE!9jQ!jv%=<}CEa1H)Gvolh1Tjh&4tz|6naq8KrDR}E zTt$ankNL4i`BgFGLg=(ie51$G&lTK*RJoJv!(@9W$o4Amq{@@s90w0N>C6&j*fod9 z5bxm`^sBN`>81RG;&!|Am9!ey-yG;&SI!RS3aFl9?nE6t*g%1dcQkirJ8tpt^;SO3ohRnhS3?65buaL6G7$6p!xZU&wgePE*ASqkO(tEK_ z#ljfK53j!j+@SA3K|YF8=z^VB(a+1^)z2Qso7IjfKF%ce0ZJZgZIx*DzGR^`A z&_7m&8F4wS3Yd!j`(04pzTzojK|Z({Q^g^$xL@IJ(Va@F8NI;Q?0U=S!&Ju0=i2k+8s?&@u9gD{=46=6tpqSQ1K^pDe@b zO!X?->zV7@t=4uO_$=?yEX;@P7Vao^KrMVPI{WLFAjrF&u)}jc*m0wYXno;aEyHo(DcE!zx{5^5whr{nCZ_QT za==&KNY}e9m~-Y2dKs0^wqZz+w;Tm+6)ux2bcvJ@qW|cna+Vv zmINoMDPw*o>J!z`ScmY%{$%V&a_aU_uF~?M~C( z>(8vr1U^!O=&Hj6i(~bEbPl+?#a4=@-Bv}Vs8o)qY)aFVMD>b#MF~nFVGKF*aq7Kf zkPk(OcVekvcj6uKyZr(@>kbhk3;wS@;!icYDFW;C1vTv=6s$kPhbv6C!&J;H}m$g(-f3B0~*bHaJ9m`y^`>W zJg~GUz|tuV^CrwB4Z`K9&Hrmg?fnKm;Bb0j-RaZygFV}t_2VPYJqj~jCyKWjoUG89oI~DzC8`M+rsEhB|SwVht(5rk!hp03wuQPK~g9khcu3IT)_C;pN zsWe`SDR!2o(zE#lo?@Vkc_e7mjv=_S1Jc-4VdiZ{`Yju{A5QUr7q>tZuhs^++pU0Z(Os}vcyub4F zDys5(MtGi&(HtaGhdZD)kLZ=}WA)tPX@)f)<>)R}6lRF&2=pl^VYS}%>6`Yj5f@Az{`difb(>)wPp(u6Z0PLMmCgr8T5E_ij+0^1px z2JDRu=;n9lJ$v(OAS~fk^yCim+5_YG7PBS}wnYbKN_ILo_uw3VViz&sh@Ay}nhyKw z7gkaN>v0&PGZfD8aCnjZc&`pzVL5mL1zrB|NPalKNl>m*Xx6!`L{FznVu(kDp0#@oQMkec6ZdF)x=>88m}s9mn%7!M>4^ zlPwguL?6>KhLw!N6rRW-=0_?VoC}QgMCQ<3#^ebR%i?7+M+6vTQOp-V)$uj{?I0t3 zmY)v0?C+u1vz2GMf+w~Vw(UID;8eQqBZ)n|S?P7)A%(}~OU5A!qZ7l3c!(D%@J>ST zw?44$65X3ZVJTjA<-BM=(#Z|LS!*!QOY^iZil2xDBgAZQZoBhry0YIDgN3{c1@M=Q z#vLa&YhNHjc67=U$(GPz8O7C?CAR#{dI_jo3@bo9xTmX|@TwES^G`S=W|qCmcp}Vx`3bzC4%|~nh)d&U-4YjbPl@nEZ@@%f zM^-TfUQH|1iC!XQI%a8S*3v$@4z*Eyt;0CEE5PyROKHzhZp$!a92c zS1q>~!MQTM;N^{j{riC3^EqSt6^zhT@N9kXK%2wYYR25k$DULZX4pkg&f@d{w!q}N z%Z!-Ld^*f%U8n0BK3!lLzy9RDSHg^T7QBsfQuz6{%W8Bd@$Ion{?{%U^e0Ia?R zthv$bbKks|kU8|m$wh3xD8xyp<-1B((77?n1i8C3 zNqolrB*U8AAZ{dveg=QMNkm%)>!%jo|9ND|lD!xuaVpHL)$AkL*;{@CqsRcGaEMTt zoU4q~L+TFCIzSdM03O?0VILg+DX2LQfa`8zPk#yP>nJPsBOI^v?4keRk2C{TXCY?H zV&Mp<+T37Qe#lOphnTeqM#B-*7yqz+f8vQJ(Gv}`&_^%=@yz8fxX(QB`#KO|CZl#Z z#=iaq^E4BuS?y=WbcJ=%jMy`sD?G+Ln#}jvs4QRcx!2gcPLMZkBHNqnj3&BHg1NJt zUH7z83=YF8aRKbMD_U*O9x|CrGS0V()CFnZD5aVU?%|1kzW{k~sL1f2kkH+yN{gMt;SSGkwLHQzGhIFDyMp)Ll1UoAu zqq>J!c#YbxJK0G+GMQuW6c!MnS`zs$GRB)=XegXe(h^2;A-*o*dY>~!9bj?vbTq1m zt-N-jU%1V4`e;Yfb!$MDJez$i3+s5K_)%V?rFz3&>chPGwUzQ1F&~Wa%xD;{qC;*X zK65iW9jtWb*T!@GGM&?X#$0nL$lf_;3rd9*WHfb!5D{~_uv07xM)()}(~at3rLXM3 z^^)P+{ta8iK^@Q;4oMLw!G34;Mrode%g$7*F_8@x!OpR@1u$6Jgtb}W36F8&)y6c47 zVzd&_Y+s?s^^ubz621M^ic)}nj_dePe!nCKj&pyt&grLvvBCDi=AkKi&}eVXwiBFg zHU@$M z|AGJ6N&CL@62G6pVi4BkZ!lPX=IaTnt4H#B**$Suer1ff(3xPxgZ3Xchr#HX zkE7on>lj$oY#bZA!xufo=`$VW74j50x7X7zFa2wc47?2zO-v-Rm?YpXx@wF-$xtRZy3B0}%zHJ<<0fByW;DTbK_U zZkn@_T&}z@-dT;0Y+P#FK=#zr!M|W^+_4IiRhI^7N)!i5+2ll&x|@`Wpl=_v@i0y! zw5xJ0QF0ekTl!FzU&n1Eqy64oYb?~8>L2xH=0|I|gOdhoKrc0OPBN=~@+bLArIJz` zH>h;VI^1m%#o@4FKar*85oh3*)51P!nnqQ)b5+r2f1kGZ7a{wKbIUPL-bqtIq5q^q zvypr%5nkW#?tSoro_rW{VG9vpA9 zv=y?}667j7gfFQsno0TOn)uw6)J|!U_%&=&7s~;mpK}O2qqN--eh3Ua`@UJuD8@Mp zmGMCD;&ee-Yf*vR;Oj7)zqZJVQdQ}OGhBKU>P@8-_>W)1s&wtC?j-Sulq7C)XIVec zZWqo!iTHj)vvmQ7jgwgQl{$)&>Fd{GU{W=NS zK8kZu?&8|B2~Na#r@4I)7t+a2Z$ZGRB&#seZg2YNTn%zl#Y*)5FN;<1Vr|7f+XTFa zAchCqa_awgX-u5}SzNLsvR@fgj8Dlq?pfEZCU!P@w@=}NjdOd$r^*GJbQSySpUjg{ zu)^}TNNFM-5x!`H%yLn@gZ?%i`nIM(l#N=z7)zASo(QR zO5tE%!?yCXoAiWt@Jxu7PD|zFyK*wjV-rC;pwL5Z>&=-$6g}J$fn(a@Vdz4!blkDKdZ?Gyl8IU;s~!BtPbBgZRCl+#MGau>GXA95>cm@vaBX0PUKhI}|jTxA6g zaPBf|vr1*eFNHYrk&{AQ{^vU?)rBy>R)JDaa2C3wg{vyfy%|IbG3%v zMc{n1JcpD|;#~JB$nhY!Xp`-MoXf5Z37@&e?isq)8SzuAN&U7_iNaxS z5a%u(r@MJsOm*@0qdz*I&is2i(m9;j^rVX7z!PxCi}T?U=8~6!R3(tHOYG-kz$Pz| zH^#c(kf|rq#T!Vz--Y_17k_u1?7?x`kt4g*^h3!G$~yb(&9DK&stFwrgNn`zcIl(^ z2&&;lSr^PZmYuo=8P*!-n(L8TsRjLWqL0TuiazFRuXGVIpx>In4jnJ-v~L*CgUbSa zQ*s4frIy!sm`CjFuo?S`yWpy|lCR4p)%2dH-hi*M?}fHUelNV`UPPRtM0>8a!W?6? z)>nmIhYA>b%o_GBXDUyjsa$|H(_N{hOjq*aLA6MYQD4esq$a{Drv-lFTj^_UC8s!J z{{^3|E!kuzXBa1Uic+k+R9+~@%j>1jczhHk{~VjzW?;0d>eXnNUAvg#X0g9ii4QX&G&Q_qHSncZ!{ZB3T)Nt!EP z#QQE;z3myoi8`OWy|s8bQ9!K)=YEM*$y#adHje6_LyJQ(dNbppIo7U(3(g+tGd$dT zbZssvCDkLU2p@O3l0$AH4ssvTOP8z*rb_=l#p;67lt#X`hz`w8SliR&hMdgwA5~sU zX(cu4OVFyN=o{*}#c?Y=WZks7!#n?*U2-y!Z!NXbVk(FtV0^R4D34P+cEY3QFiu>p z$Z;pYAzEY?VwMzVcf-?DSHFAgO!cHty4J1w#RR+s{EUrQ(3KE^sM&P^0)PcwA*qXv7KAO zsc!#nm9mza|KQ;1hGvD%h3@H-&E|G)_``ooZ{Y) zrRP})PnE79q7%U@R5ZeEVIxmNA5sbYtFQZzY~=&1I~=hR?d2fl%|M&lfjMP^ncoPn zK>>xsaVHPzpdn&)CEi;gs&ed!n4A$))otPgDz_YP*3P*LtZGJ6Tw6mads6NNB12P+ zZ|(j(xtZu3a-t>Zp>)$Wcp(sQ29=~`kRQW`uI@Nw2bn;%E_;SVLe?LSzEo5;Q*u$-ZW=pdny$%m7(f&k%gU%@TS9BGx zB?mX63rbJw-P@8P^`T;)M}D2xtpRuB5!m}xxWn6-2bMFAoF*rz)<~+oXXpoxQkh<- z4j4lIHlO)$j+K&3cKi!->M>Pnp+!Vgj2RqTB&wBX z0bX72gjejKErm|@VIzzFbMS27Y|7Zc55b4}WGl((3M0K7{(`M}a{s7VyvKc(?~u2b z_PgAik^J51V0Wj3`!hbx&-GmT793AM>3^DwZCGvM4e2KRl0mT4v!e^m&!`Mh%PV`O z3}SEh9yz^()9Mmj;IqRRYhhOepPU32u9#Ss+~~BN1$9qB%Y82Tltt!q)iPlfwFpY_ull=u|7%SYrd&&ceSF{5_b3t2PGdG%p;_>n&k zAmh6ND{~dx-;B6pCdhG~(vjz5uEmUuc&YlOmt?-3rHkTb_cyDdu{3l$wOZhDN@~Cp z8en{B7jq}!_BL8_QSz{>Ykps8|8-x4_j~oU6cje2`SEj>$uanPt&O4jk6aX%Ib;opA7`FfRt$Q~C2gDxx;%IEoNcK2x0)!1W~x=e2L$Q$iUS+G8VH#%Qqv zBiTu$VpQ~AX&0vxeq}#1YU<5{JpyA=1_$N^%W#rH4Qj&nIH+x;Zt1B^(@uLun3;9G zuhcej1)wpX4kde zm@SzjWx%+Ex+(w$6{E)abVQ>kVE!S=zHYzKpTGto-5!|wCR)>3wyX><)1 z+J3sPHRxQ`MhQ|C1n^g(1KjJLzV%TT;Us#BUF6SkE-ub^;-(_;qc?}If`AEcxY z>`l$4?=*jK@QM&i!w{p=R&r>Sz3Y94eT%)DwT4Q7Zs;JK91k*2bk2%9qigymPD6RF zZ!#0@EA9%hh+IQCN1oDJ-3TAa#|eiK>UFt|v_*K~RIvkQ99*{7hHGRrhnn3j%U> z!%hHq?~X^t*Hpe^QI?Hjr7fmzKW&x3$0*#xo5R|vLA0*Jx~hhaxs6>`adGUg+3?H--{OASXZ1p+q;+FZyaEw#t zK*r=}kTsj=T%KyEBeD0eQ-KPhKZr{LtEVEq0sX`NkiJgQuVdek-By*;39IQ|X(FdG z+Ov(#p?txifki1@11*9F^yStO_X~1vuhLqntISj4vYhRxqyU z9%F^^4XlAJ;0^WRx&0-TQ);X4)tp*Tovfz7c#M!wh~EmIaVB14T{Gk1qs=qxv4*~& z0@_T4R)CzMle7*+@IW{mWnl8);6n9L&_0Vx^h)cP^`4GceRiPZb|i|^MowNZ$SW{i zLuC9HoQ~jy_|=PhF4 zs&rThse|!Ry{cAIFH!xZpgN0|jtHgQv3AfrXh_C8y*K=nWb=@<(!PvqdS{$ht9vSY zn|XSu>*O=w4{M}qoN74Neaudm2bNPUxW~U*1+xRirIpoIs9DtO_V8d`w_L2ws>4|Hcu z1dZk9JZp!&9Sooms^i>r^IwR48IyfdX2~Pobx*=hXl`f24gNWk$O$Dk-nylC3|&ff z1EW%R2Y-USdEPeML7d8zN`3c@@>=Pw-Ql#l(cXWx+)4+rk$V{~E-YD2DBXdQBne(q zN$#?<_1fOb$yU|mL5dH)>rG`dPU_F7l{+hK~lZ9JO694y2_3AeF&#T6p?$V z$Gp=bt3-bnxz|%s-UIgO2uI1NL8r{IdYhmXD4g0Um|IUW{)C@95JlZq9A36bUE~Ng zsIByt^TiQYqTxp`q#8XA)*7Ny5O13Bdb7iT%WICZF4&pj71e>AG#+Mef4PuiD9hBN zYE8AMvR%p~UZiV2hN|rznfF*ch|^p3L5Ux_)j;vzFbf(`W0aEnz>`V9P1Hd9+Z5-K znkbhCxG$(`df>^>22RHbry>pl52=$^kmI+dgE<=9|5v&LSE&u3(-F)ApW!U^QXjgg z=fPQ~knLUq4;q0Ie>UjD}_=WOw90)x+8)Ed`$Ic>1^_rKND(CQCnZ0$*Euj(L%& ze~`1oZjl`%!vDAfYkeRXWi)loE2#rcm0>SWvJaeeP3eI&81~m_p^f_^eeY9l9ensI zp*yXI9JZ%j(*2FxtUB7}OVqXaE0JqflUhl8a3LCl-na>FO!4%cO3~*^PcHcnDD^`- z6Npe@#{er`Km|Sv{ckO}47cfj&%jym2^e=idLmsw9Inxgh++hLh$obMz7Y}M_!nvW zrMq}X&zFy>Uuol%kWj`ZOAj5_Uz^*Ro%^kTaTBaRH}1IeoG;P!Z^L=efy0}~j5#Z3 z0xy{ZUv!*Y3=Q%VFphsfUJts8`w)a zblwU&Mc_0xu`fB3SZ8HWl5eGNJC~EsR!MR4PsGDrc+70Xt1gRp2hMC4`d(Y;2+txL zaBPVRz5zX*Z^7EL(&hdUe#KkPeVK-)s;3yky6p#FW{Q-I8nv@DMf?+dx-E)|C)}+c zmEsz_5`SQ1%7B-Yb+6I88wgwDE1ttOdR23Am#IOJ_e zz1Rbni3ePdQ=f=E&v^xp)Gr)%%G!NVE=W!j*b7!3-YL3EIjv41q6@l;(|NQbCP{OVSV& zA{8d|T%K2uc^H8QSaq0RM?mx2;1#l&&-ok0ZEJGg+A!A^)4Lo-*X0eUok4$iAWtZc zuY2jAo4ohoG#e4X^OOyyTt2KYNpzc9s7BdB*G! zpYEp*nvLttj^|AtbRX4-^?j!?%@9}r zk4zmkjU!{?|znFr_#PVOAq~DM(Hk}|CH-^$MqT9F}X4C5RIEw4)nRXdF7&8 zUl2ZHNk*+SpB4V3@aL6C`%^)FA}d#$9(8vNTuple[int, np.ndarray]: + """load wav file. + + Args: + wavpath (str): wav path. + sr (int, optional): expect sample rate. Defaults to None. + dtype (str, optional): wav data bits. Defaults to 'int16'. + + Returns: + Tuple[int, np.ndarray]: sr (int), wav (int16) [T, C]. + """ + wav, r_sr = sf.read(wavpath, start=start, stop=stop, dtype=dtype, always_2d=always_2d) + if sr: + assert sr == r_sr + return r_sr, wav + + +def write(wavpath:str, wav:np.ndarray, sr:int, dtype='PCM_16'): + """write wav file. + + Args: + wavpath (str): file path to save. + wav (np.ndarray): wav data. + sr (int): data samplerate. + dtype (str, optional): wav bit format. Defaults to 'PCM_16'. + """ + sf.write(wavpath, wav, sr, subtype=dtype) + + +def frames(x: Tensor, + num_samples: Tensor, + sr: int, + win_length: float, + stride_length: float, + clip: bool = False) -> Tuple[Tensor, Tensor]: + """Extract frames from audio. + + Parameters + ---------- + x : Tensor + Shape (B, T), batched waveform. + num_samples : Tensor + Shape (B, ), number of samples of each waveform. + sr: int + Sampling Rate. + win_length : float + Window length in ms. + stride_length : float + Stride length in ms. + clip : bool, optional + Whether to clip audio that does not fit into the last frame, by + default True + + Returns + ------- + frames : Tensor + Shape (B, T', win_length). + num_frames : Tensor + Shape (B, ) number of valid frames + """ + assert stride_length <= win_length + stride_length = int(stride_length * sr) + win_length = int(win_length * sr) + + num_frames = (num_samples - win_length) // stride_length + padding = (0, 0) + if not clip: + num_frames += 1 + need_samples = num_frames * stride_length + win_length + padding = (0, need_samples - num_samples - 1) + + weight = paddle.eye(win_length).unsqueeze(1) #[win_length, 1, win_length] + + frames = F.conv1d(x.unsqueeze(-1), + weight, + padding=padding, + stride=(stride_length, ), + data_format='NLC') + return frames, num_frames + + +def dither(signal:Tensor, dither_value=1.0)->Tensor: + """dither frames for log compute. + + Args: + signal (Tensor): [B, T, D] + dither_value (float, optional): [scalar]. Defaults to 1.0. + + Returns: + Tensor: [B, T, D] + """ + D = paddle.shape(signal)[-1] + signal += paddle.normal(shape=[1, 1, D]) * dither_value + return signal + + +def remove_dc_offset(signal:Tensor)->Tensor: + """remove dc. + + Args: + signal (Tensor): [B, T, D] + + Returns: + Tensor: [B, T, D] + """ + signal -= paddle.mean(signal, axis=-1, keepdim=True) + return signal + +def preemphasis(signal:Tensor, coeff=0.97)->Tensor: + """perform preemphasis on the input signal. + + Args: + signal (Tensor): [B, T, D], The signal to filter. + coeff (float, optional): [scalar].The preemphasis coefficient. 0 is no filter, Defaults to 0.97. + + Returns: + Tensor: [B, T, D] + """ + return paddle.concat([ + (1-coeff)*signal[:, :, 0:1], + signal[:, :, 1:] - coeff * signal[:, :, :-1] + ], axis=-1) + + +class STFT(nn.Layer): + """A module for computing stft transformation in a differentiable way. + + http://practicalcryptography.com/miscellaneous/machine-learning/intuitive-guide-discrete-fourier-transform/ + + Parameters + ------------ + n_fft : int + Number of samples in a frame. + + sr: int + Number of Samplilng rate. + + stride_length : float + Number of samples shifted between adjacent frames. + + win_length : float + Length of the window. + + clip: bool + Whether to clip audio is necesaary. + """ + def __init__(self, + n_fft: int, + sr: int, + win_length: float, + stride_length: float, + dither:float=0.0, + preemph_coeff:float=0.97, + remove_dc_offset:bool=True, + window_type: str = 'povey', + clip: bool = False): + super().__init__() + self.sr = sr + self.win_length = win_length + self.stride_length = stride_length + self.dither = dither + self.preemph_coeff = preemph_coeff + self.remove_dc_offset = remove_dc_offset + self.window_type = window_type + self.clip = clip + + self.n_fft = n_fft + self.n_bin = 1 + n_fft // 2 + + w_real, w_imag, kernel_size = dft_matrix( + self.n_fft, int(self.win_length * self.sr), self.n_bin + ) + + # calculate window + window = get_window(window_type, kernel_size) + + # (2 * n_bins, kernel_size) + w = np.concatenate([w_real, w_imag], axis=0) + w = w * window + # (kernel_size, 2 * n_bins) + w = np.transpose(w) + weight = paddle.cast(paddle.to_tensor(w), paddle.get_default_dtype()) + self.register_buffer("weight", weight) + + def forward(self, x: Tensor, num_samples: Tensor) -> Tuple[Tensor, Tensor]: + """Compute the stft transform. + Parameters + ------------ + x : Tensor [shape=(B, T)] + The input waveform. + num_samples : Tensor [shape=(B,)] + Number of samples of each waveform. + Returns + ------------ + C : Tensor + Shape(B, T', n_bins, 2) Spectrogram. + + num_frames: Tensor + Shape (B,) number of samples of each spectrogram + """ + batch_size = paddle.shape(num_samples) + F, nframe = frames(x, num_samples, self.sr, self.win_length, self.stride_length, clip=self.clip) + if self.dither: + F = dither(F, self.dither) + if self.remove_dc_offset: + F = remove_dc_offset(F) + if self.preemph_coeff: + F = preemphasis(F) + C = paddle.matmul(F, self.weight) # [B, T, K] [K, 2 * n_bins] + C = paddle.reshape(C, [batch_size, -1, 2, self.n_bin]) + C = C.transpose([0, 1, 3, 2]) + return C, nframe + + +def powspec(C:Tensor) -> Tensor: + """Compute the power spectrum |X_k|^2. + + Args: + C (Tensor): [B, T, C, 2] + + Returns: + Tensor: [B, T, C] + """ + real, imag = paddle.chunk(C, 2, axis=-1) + return paddle.square(real.squeeze(-1)) + paddle.square(imag.squeeze(-1)) + + +def magspec(C: Tensor, eps=1e-10) -> Tensor: + """Compute the magnitude spectrum |X_k|. + + Args: + C (Tensor): [B, T, C, 2] + eps (float): epsilon. + + Returns: + Tensor: [B, T, C] + """ + pspec = powspec(C) + return paddle.sqrt(pspec + eps) + + +def logspec(C: Tensor, eps=1e-10) -> Tensor: + """Compute log-spectrum 20log10∣X_k∣. + + Args: + C (Tensor): [description] + eps ([type], optional): [description]. Defaults to 1e-10. + + Returns: + Tensor: [description] + """ + spec = magspec(C) + return 20 * paddle.log10(spec + eps) + diff --git a/third_party/paddle_audio/frontend/kaldi_test.py b/third_party/paddle_audio/frontend/kaldi_test.py new file mode 100644 index 00000000..34ff413c --- /dev/null +++ b/third_party/paddle_audio/frontend/kaldi_test.py @@ -0,0 +1,533 @@ +from typing import Tuple +import numpy as np +import paddle +import unittest + +import decimal +import numpy +import math +import logging +from pathlib import Path + +from scipy.fftpack import dct + +from third_party.paddle_audio.frontend import kaldi + +def round_half_up(number): + return int(decimal.Decimal(number).quantize(decimal.Decimal('1'), rounding=decimal.ROUND_HALF_UP)) + +def rolling_window(a, window, step=1): + # http://ellisvalentiner.com/post/2017-03-21-np-strides-trick + shape = a.shape[:-1] + (a.shape[-1] - window + 1, window) + strides = a.strides + (a.strides[-1],) + return numpy.lib.stride_tricks.as_strided(a, shape=shape, strides=strides)[::step] + + +def do_dither(signal, dither_value=1.0): + signal += numpy.random.normal(size=signal.shape) * dither_value + return signal + +def do_remove_dc_offset(signal): + signal -= numpy.mean(signal) + return signal + +def do_preemphasis(signal, coeff=0.97): + """perform preemphasis on the input signal. + + :param signal: The signal to filter. + :param coeff: The preemphasis coefficient. 0 is no filter, default is 0.95. + :returns: the filtered signal. + """ + return numpy.append((1-coeff)*signal[0], signal[1:] - coeff * signal[:-1]) + + +def framesig(sig, frame_len, frame_step, dither=1.0, preemph=0.97, remove_dc_offset=True, wintype='hamming', stride_trick=True): + """Frame a signal into overlapping frames. + + :param sig: the audio signal to frame. + :param frame_len: length of each frame measured in samples. + :param frame_step: number of samples after the start of the previous frame that the next frame should begin. + :param winfunc: the analysis window to apply to each frame. By default no window is applied. + :param stride_trick: use stride trick to compute the rolling window and window multiplication faster + :returns: an array of frames. Size is NUMFRAMES by frame_len. + """ + slen = len(sig) + frame_len = int(round_half_up(frame_len)) + frame_step = int(round_half_up(frame_step)) + if slen <= frame_len: + numframes = 1 + else: + numframes = 1 + (( slen - frame_len) // frame_step) + + # check kaldi/src/feat/feature-window.h + padsignal = sig[:(numframes-1)*frame_step+frame_len] + if wintype is 'povey': + win = numpy.empty(frame_len) + for i in range(frame_len): + win[i] = (0.5-0.5*numpy.cos(2*numpy.pi/(frame_len-1)*i))**0.85 + else: # the hamming window + win = numpy.hamming(frame_len) + + if stride_trick: + frames = rolling_window(padsignal, window=frame_len, step=frame_step) + else: + indices = numpy.tile(numpy.arange(0, frame_len), (numframes, 1)) + numpy.tile( + numpy.arange(0, numframes * frame_step, frame_step), (frame_len, 1)).T + indices = numpy.array(indices, dtype=numpy.int32) + frames = padsignal[indices] + win = numpy.tile(win, (numframes, 1)) + + frames = frames.astype(numpy.float32) + raw_frames = numpy.zeros(frames.shape) + for frm in range(frames.shape[0]): + frames[frm,:] = do_dither(frames[frm,:], dither) # dither + frames[frm,:] = do_remove_dc_offset(frames[frm,:]) # remove dc offset + raw_frames[frm,:] = frames[frm,:] + frames[frm,:] = do_preemphasis(frames[frm,:], preemph) # preemphasize + + return frames * win, raw_frames + + +def magspec(frames, NFFT): + """Compute the magnitude spectrum of each frame in frames. If frames is an NxD matrix, output will be Nx(NFFT/2+1). + + :param frames: the array of frames. Each row is a frame. + :param NFFT: the FFT length to use. If NFFT > frame_len, the frames are zero-padded. + :returns: If frames is an NxD matrix, output will be Nx(NFFT/2+1). Each row will be the magnitude spectrum of the corresponding frame. + """ + if numpy.shape(frames)[1] > NFFT: + logging.warn( + 'frame length (%d) is greater than FFT size (%d), frame will be truncated. Increase NFFT to avoid.', + numpy.shape(frames)[1], NFFT) + complex_spec = numpy.fft.rfft(frames, NFFT) + return numpy.absolute(complex_spec) + + +def powspec(frames, NFFT): + """Compute the power spectrum of each frame in frames. If frames is an NxD matrix, output will be Nx(NFFT/2+1). + + :param frames: the array of frames. Each row is a frame. + :param NFFT: the FFT length to use. If NFFT > frame_len, the frames are zero-padded. + :returns: If frames is an NxD matrix, output will be Nx(NFFT/2+1). Each row will be the power spectrum of the corresponding frame. + """ + return numpy.square(magspec(frames, NFFT)) + + + +def mfcc(signal,samplerate=16000,winlen=0.025,winstep=0.01,numcep=13, + nfilt=23,nfft=512,lowfreq=20,highfreq=None,dither=1.0,remove_dc_offset=True,preemph=0.97, + ceplifter=22,useEnergy=True,wintype='povey'): + """Compute MFCC features from an audio signal. + + :param signal: the audio signal from which to compute features. Should be an N*1 array + :param samplerate: the samplerate of the signal we are working with. + :param winlen: the length of the analysis window in seconds. Default is 0.025s (25 milliseconds) + :param winstep: the step between successive windows in seconds. Default is 0.01s (10 milliseconds) + :param numcep: the number of cepstrum to return, default 13 + :param nfilt: the number of filters in the filterbank, default 26. + :param nfft: the FFT size. Default is 512. + :param lowfreq: lowest band edge of mel filters. In Hz, default is 0. + :param highfreq: highest band edge of mel filters. In Hz, default is samplerate/2 + :param preemph: apply preemphasis filter with preemph as coefficient. 0 is no filter. Default is 0.97. + :param ceplifter: apply a lifter to final cepstral coefficients. 0 is no lifter. Default is 22. + :param appendEnergy: if this is true, the zeroth cepstral coefficient is replaced with the log of the total frame energy. + :param winfunc: the analysis window to apply to each frame. By default no window is applied. You can use numpy window functions here e.g. winfunc=numpy.hamming + :returns: A numpy array of size (NUMFRAMES by numcep) containing features. Each row holds 1 feature vector. + """ + feat,energy = fbank(signal,samplerate,winlen,winstep,nfilt,nfft,lowfreq,highfreq,dither,remove_dc_offset,preemph,wintype) + feat = numpy.log(feat) + feat = dct(feat, type=2, axis=1, norm='ortho')[:,:numcep] + feat = lifter(feat,ceplifter) + if useEnergy: feat[:,0] = numpy.log(energy) # replace first cepstral coefficient with log of frame energy + return feat + +def fbank(signal,samplerate=16000,winlen=0.025,winstep=0.01, + nfilt=40,nfft=512,lowfreq=0,highfreq=None,dither=1.0,remove_dc_offset=True, preemph=0.97, + wintype='hamming'): + """Compute Mel-filterbank energy features from an audio signal. + + :param signal: the audio signal from which to compute features. Should be an N*1 array + :param samplerate: the samplerate of the signal we are working with. + :param winlen: the length of the analysis window in seconds. Default is 0.025s (25 milliseconds) + :param winstep: the step between successive windows in seconds. Default is 0.01s (10 milliseconds) + :param nfilt: the number of filters in the filterbank, default 26. + :param nfft: the FFT size. Default is 512. + :param lowfreq: lowest band edge of mel filters. In Hz, default is 0. + :param highfreq: highest band edge of mel filters. In Hz, default is samplerate/2 + :param preemph: apply preemphasis filter with preemph as coefficient. 0 is no filter. Default is 0.97. + :param winfunc: the analysis window to apply to each frame. By default no window is applied. You can use numpy window functions here e.g. winfunc=numpy.hamming + winfunc=lambda x:numpy.ones((x,)) + :returns: 2 values. The first is a numpy array of size (NUMFRAMES by nfilt) containing features. Each row holds 1 feature vector. The + second return value is the energy in each frame (total energy, unwindowed) + """ + highfreq= highfreq or samplerate/2 + frames,raw_frames = sigproc.framesig(signal, winlen*samplerate, winstep*samplerate, dither, preemph, remove_dc_offset, wintype) + pspec = sigproc.powspec(frames,nfft) # nearly the same until this part + energy = numpy.sum(raw_frames**2,1) # this stores the raw energy in each frame + energy = numpy.where(energy == 0,numpy.finfo(float).eps,energy) # if energy is zero, we get problems with log + + fb = get_filterbanks(nfilt,nfft,samplerate,lowfreq,highfreq) + feat = numpy.dot(pspec,fb.T) # compute the filterbank energies + feat = numpy.where(feat == 0,numpy.finfo(float).eps,feat) # if feat is zero, we get problems with log + + return feat,energy + +def logfbank(signal,samplerate=16000,winlen=0.025,winstep=0.01, + nfilt=40,nfft=512,lowfreq=64,highfreq=None,dither=1.0,remove_dc_offset=True,preemph=0.97,wintype='hamming'): + """Compute log Mel-filterbank energy features from an audio signal. + + :param signal: the audio signal from which to compute features. Should be an N*1 array + :param samplerate: the samplerate of the signal we are working with. + :param winlen: the length of the analysis window in seconds. Default is 0.025s (25 milliseconds) + :param winstep: the step between successive windows in seconds. Default is 0.01s (10 milliseconds) + :param nfilt: the number of filters in the filterbank, default 26. + :param nfft: the FFT size. Default is 512. + :param lowfreq: lowest band edge of mel filters. In Hz, default is 0. + :param highfreq: highest band edge of mel filters. In Hz, default is samplerate/2 + :param preemph: apply preemphasis filter with preemph as coefficient. 0 is no filter. Default is 0.97. + :returns: A numpy array of size (NUMFRAMES by nfilt) containing features. Each row holds 1 feature vector. + """ + feat,energy = fbank(signal,samplerate,winlen,winstep,nfilt,nfft,lowfreq,highfreq,dither, remove_dc_offset,preemph,wintype) + return numpy.log(feat) + +def hz2mel(hz): + """Convert a value in Hertz to Mels + + :param hz: a value in Hz. This can also be a numpy array, conversion proceeds element-wise. + :returns: a value in Mels. If an array was passed in, an identical sized array is returned. + """ + return 1127 * numpy.log(1+hz/700.0) + +def mel2hz(mel): + """Convert a value in Mels to Hertz + + :param mel: a value in Mels. This can also be a numpy array, conversion proceeds element-wise. + :returns: a value in Hertz. If an array was passed in, an identical sized array is returned. + """ + return 700 * (numpy.exp(mel/1127.0)-1) + +def get_filterbanks(nfilt=26,nfft=512,samplerate=16000,lowfreq=0,highfreq=None): + """Compute a Mel-filterbank. The filters are stored in the rows, the columns correspond + to fft bins. The filters are returned as an array of size nfilt * (nfft/2 + 1) + + :param nfilt: the number of filters in the filterbank, default 20. + :param nfft: the FFT size. Default is 512. + :param samplerate: the samplerate of the signal we are working with. Affects mel spacing. + :param lowfreq: lowest band edge of mel filters, default 0 Hz + :param highfreq: highest band edge of mel filters, default samplerate/2 + :returns: A numpy array of size nfilt * (nfft/2 + 1) containing filterbank. Each row holds 1 filter. + """ + highfreq= highfreq or samplerate/2 + assert highfreq <= samplerate/2, "highfreq is greater than samplerate/2" + + # compute points evenly spaced in mels + lowmel = hz2mel(lowfreq) + highmel = hz2mel(highfreq) + + # check kaldi/src/feat/Mel-computations.h + fbank = numpy.zeros([nfilt,nfft//2+1]) + mel_freq_delta = (highmel-lowmel)/(nfilt+1) + for j in range(0,nfilt): + leftmel = lowmel+j*mel_freq_delta + centermel = lowmel+(j+1)*mel_freq_delta + rightmel = lowmel+(j+2)*mel_freq_delta + for i in range(0,nfft//2): + mel=hz2mel(i*samplerate/nfft) + if mel>leftmel and mel 0: + nframes,ncoeff = numpy.shape(cepstra) + n = numpy.arange(ncoeff) + lift = 1 + (L/2.)*numpy.sin(numpy.pi*n/L) + return lift*cepstra + else: + # values of L <= 0, do nothing + return cepstra + +def delta(feat, N): + """Compute delta features from a feature vector sequence. + + :param feat: A numpy array of size (NUMFRAMES by number of features) containing features. Each row holds 1 feature vector. + :param N: For each frame, calculate delta features based on preceding and following N frames + :returns: A numpy array of size (NUMFRAMES by number of features) containing delta features. Each row holds 1 delta feature vector. + """ + if N < 1: + raise ValueError('N must be an integer >= 1') + NUMFRAMES = len(feat) + denominator = 2 * sum([i**2 for i in range(1, N+1)]) + delta_feat = numpy.empty_like(feat) + padded = numpy.pad(feat, ((N, N), (0, 0)), mode='edge') # padded version of feat + for t in range(NUMFRAMES): + delta_feat[t] = numpy.dot(numpy.arange(-N, N+1), padded[t : t+2*N+1]) / denominator # [t : t+2*N+1] == [(N+t)-N : (N+t)+N+1] + return delta_feat + +##### modify for test ###### + +def framesig_without_dither_dc_preemphasize(sig, frame_len, frame_step, wintype='hamming', stride_trick=True): + """Frame a signal into overlapping frames. + + :param sig: the audio signal to frame. + :param frame_len: length of each frame measured in samples. + :param frame_step: number of samples after the start of the previous frame that the next frame should begin. + :param winfunc: the analysis window to apply to each frame. By default no window is applied. + :param stride_trick: use stride trick to compute the rolling window and window multiplication faster + :returns: an array of frames. Size is NUMFRAMES by frame_len. + """ + slen = len(sig) + frame_len = int(round_half_up(frame_len)) + frame_step = int(round_half_up(frame_step)) + if slen <= frame_len: + numframes = 1 + else: + numframes = 1 + (( slen - frame_len) // frame_step) + + # check kaldi/src/feat/feature-window.h + padsignal = sig[:(numframes-1)*frame_step+frame_len] + + if wintype is 'povey': + win = numpy.empty(frame_len) + for i in range(frame_len): + win[i] = (0.5-0.5*numpy.cos(2*numpy.pi/(frame_len-1)*i))**0.85 + elif wintype == '': + win = numpy.ones(frame_len) + elif wintype == 'hann': + win = numpy.hanning(frame_len) + else: # the hamming window + win = numpy.hamming(frame_len) + + if stride_trick: + frames = rolling_window(padsignal, window=frame_len, step=frame_step) + else: + indices = numpy.tile(numpy.arange(0, frame_len), (numframes, 1)) + numpy.tile( + numpy.arange(0, numframes * frame_step, frame_step), (frame_len, 1)).T + indices = numpy.array(indices, dtype=numpy.int32) + frames = padsignal[indices] + win = numpy.tile(win, (numframes, 1)) + + frames = frames.astype(numpy.float32) + raw_frames = frames + return frames * win, raw_frames + + +def frames(signal,samplerate=16000,winlen=0.025,winstep=0.01, + nfilt=40,nfft=512,lowfreq=0,highfreq=None, wintype='hamming'): + frames_with_win, raw_frames = framesig_without_dither_dc_preemphasize(signal, winlen*samplerate, winstep*samplerate, wintype) + return frames_with_win, raw_frames + + +def complexspec(frames, NFFT): + """Compute the magnitude spectrum of each frame in frames. If frames is an NxD matrix, output will be Nx(NFFT/2+1). + + :param frames: the array of frames. Each row is a frame. + :param NFFT: the FFT length to use. If NFFT > frame_len, the frames are zero-padded. + :returns: If frames is an NxD matrix, output will be Nx(NFFT/2+1). Each row will be the magnitude spectrum of the corresponding frame. + """ + if numpy.shape(frames)[1] > NFFT: + logging.warn( + 'frame length (%d) is greater than FFT size (%d), frame will be truncated. Increase NFFT to avoid.', + numpy.shape(frames)[1], NFFT) + complex_spec = numpy.fft.rfft(frames, NFFT) + return complex_spec + + +def stft_with_window(signal,samplerate=16000,winlen=0.025,winstep=0.01, + nfilt=40,nfft=512,lowfreq=0,highfreq=None,dither=1.0,remove_dc_offset=True, preemph=0.97, + wintype='hamming'): + frames_with_win, raw_frames = framesig_without_dither_dc_preemphasize(signal, winlen*samplerate, winstep*samplerate, wintype) + + spec = magspec(frames_with_win, nfft) # nearly the same until this part + scomplex = complexspec(frames_with_win, nfft) + + rspec = magspec(raw_frames, nfft) + rcomplex = complexspec(raw_frames, nfft) + return spec, scomplex, rspec, rcomplex + + +class TestKaldiFE(unittest.TestCase): + def setUp(self): + self. this_dir = Path(__file__).parent + + self.wavpath = str(self.this_dir / 'english.wav') + self.winlen=0.025 # ms + self.winstep=0.01 # ms + self.nfft=512 + self.lowfreq = 0 + self.highfreq = None + self.wintype='hamm' + self.nfilt=40 + + paddle.set_device('cpu') + + + def test_read(self): + import scipy.io.wavfile as wav + rate, sig = wav.read(self.wavpath) + sr, wav = kaldi.read(self.wavpath) + wav = wav[:, 0] + self.assertTrue(np.all(sig == wav)) + self.assertEqual(rate, sr) + + def test_frames(self): + sr, wav = kaldi.read(self.wavpath) + wav = wav[:, 0] + _, fs = frames(wav, samplerate=sr, + winlen=self.winlen, winstep=self.winstep, + nfilt=self.nfilt, nfft=self.nfft, + lowfreq=self.lowfreq, highfreq=self.highfreq, + wintype=self.wintype) + + t_wav = paddle.to_tensor([wav], dtype='float32') + t_wavlen = paddle.to_tensor([len(wav)]) + t_fs, t_nframe = kaldi.frames(t_wav, t_wavlen, sr, self.winlen, self.winstep, clip=False) + t_fs = t_fs.astype(fs.dtype)[0] + + self.assertEqual(t_nframe.item(), fs.shape[0]) + self.assertTrue(np.allclose(t_fs.numpy(), fs)) + + + def test_stft(self): + sr, wav = kaldi.read(self.wavpath) + wav = wav[:, 0] + + for wintype in ['', 'hamm', 'hann', 'povey']: + self.wintype=wintype + _, stft_c_win, _, _ = stft_with_window(wav, samplerate=sr, + winlen=self.winlen, winstep=self.winstep, + nfilt=self.nfilt, nfft=self.nfft, + lowfreq=self.lowfreq, highfreq=self.highfreq, + wintype=self.wintype) + + t_wav = paddle.to_tensor([wav], dtype='float32') + t_wavlen = paddle.to_tensor([len(wav)]) + + stft_class = kaldi.STFT(self.nfft, sr, self.winlen, self.winstep, window_type=self.wintype, dither=0.0, preemph_coeff=0.0, remove_dc_offset=False, clip=False) + t_stft, t_nframe = stft_class(t_wav, t_wavlen) + t_stft = t_stft.astype(stft_c_win.real.dtype)[0] + t_real = t_stft[:, :, 0] + t_imag = t_stft[:, :, 1] + + self.assertEqual(t_nframe.item(), stft_c_win.real.shape[0]) + + self.assertLess(np.sum(t_real.numpy()) - np.sum(stft_c_win.real), 1) + self.assertTrue(np.allclose(t_real.numpy(), stft_c_win.real, atol=1e-1)) + + self.assertLess(np.sum(t_imag.numpy()) - np.sum(stft_c_win.imag), 1) + self.assertTrue(np.allclose(t_imag.numpy(), stft_c_win.imag, atol=1e-1)) + + + def test_magspec(self): + sr, wav = kaldi.read(self.wavpath) + wav = wav[:, 0] + for wintype in ['', 'hamm', 'hann', 'povey']: + self.wintype=wintype + stft_win, _, _, _ = stft_with_window(wav, samplerate=sr, + winlen=self.winlen, winstep=self.winstep, + nfilt=self.nfilt, nfft=self.nfft, + lowfreq=self.lowfreq, highfreq=self.highfreq, + wintype=self.wintype) + + t_wav = paddle.to_tensor([wav], dtype='float32') + t_wavlen = paddle.to_tensor([len(wav)]) + + stft_class = kaldi.STFT(self.nfft, sr, self.winlen, self.winstep, window_type=self.wintype, dither=0.0, preemph_coeff=0.0, remove_dc_offset=False, clip=False) + t_stft, t_nframe = stft_class(t_wav, t_wavlen) + t_stft = t_stft.astype(stft_win.dtype) + t_spec = kaldi.magspec(t_stft)[0] + + self.assertEqual(t_nframe.item(), stft_win.shape[0]) + + self.assertLess(np.sum(t_spec.numpy()) - np.sum(stft_win), 1) + self.assertTrue(np.allclose(t_spec.numpy(), stft_win, atol=1e-1)) + + + def test_magsepc_winprocess(self): + sr, wav = kaldi.read(self.wavpath) + wav = wav[:, 0] + fs, _= framesig(wav, self.winlen*sr, self.winstep*sr, + dither=0.0, preemph=0.97, remove_dc_offset=True, wintype='povey', stride_trick=True) + spec = magspec(fs, self.nfft) # nearly the same until this part + + t_wav = paddle.to_tensor([wav], dtype='float32') + t_wavlen = paddle.to_tensor([len(wav)]) + stft_class = kaldi.STFT( + self.nfft, sr, self.winlen, self.winstep, + window_type='povey', dither=0.0, preemph_coeff=0.97, remove_dc_offset=True, clip=False) + t_stft, t_nframe = stft_class(t_wav, t_wavlen) + t_stft = t_stft.astype(spec.dtype) + t_spec = kaldi.magspec(t_stft)[0] + + self.assertEqual(t_nframe.item(), fs.shape[0]) + + self.assertLess(np.sum(t_spec.numpy()) - np.sum(spec), 1) + self.assertTrue(np.allclose(t_spec.numpy(), spec, atol=1e-1)) + + + def test_powspec(self): + sr, wav = kaldi.read(self.wavpath) + wav = wav[:, 0] + for wintype in ['', 'hamm', 'hann', 'povey']: + self.wintype=wintype + stft_win, _, _, _ = stft_with_window(wav, samplerate=sr, + winlen=self.winlen, winstep=self.winstep, + nfilt=self.nfilt, nfft=self.nfft, + lowfreq=self.lowfreq, highfreq=self.highfreq, + wintype=self.wintype) + stft_win = np.square(stft_win) + + t_wav = paddle.to_tensor([wav], dtype='float32') + t_wavlen = paddle.to_tensor([len(wav)]) + + stft_class = kaldi.STFT(self.nfft, sr, self.winlen, self.winstep, window_type=self.wintype, dither=0.0, preemph_coeff=0.0, remove_dc_offset=False, clip=False) + t_stft, t_nframe = stft_class(t_wav, t_wavlen) + t_stft = t_stft.astype(stft_win.dtype) + t_spec = kaldi.powspec(t_stft)[0] + + self.assertEqual(t_nframe.item(), stft_win.shape[0]) + + self.assertLess(np.sum(t_spec.numpy() - stft_win), 5e4) + self.assertTrue(np.allclose(t_spec.numpy(), stft_win, atol=1e2)) + + +# from python_speech_features import mfcc +# from python_speech_features import delta +# from python_speech_features import logfbank +# import scipy.io.wavfile as wav + +# (rate,sig) = wav.read("english.wav") + +# # note that generally nfilt=40 is used for speech recognition +# fbank_feat = logfbank(sig,nfilt=23,lowfreq=20,dither=0,wintype='povey') + +# # the computed fbank coefficents of english.wav with dimension [110,23] +# # [ 12.2865 12.6906 13.1765 15.714 16.064 15.7553 16.5746 16.9205 16.6472 16.1302 16.4576 16.7326 16.8864 17.7215 18.88 19.1377 19.1495 18.6683 18.3886 20.3506 20.2772 18.8248 18.1899 +# # 11.9198 13.146 14.7215 15.8642 17.4288 16.394 16.8238 16.1095 16.4297 16.6331 16.3163 16.5093 17.4981 18.3429 19.6555 19.6263 19.8435 19.0534 19.001 20.0287 19.7707 19.5852 19.1112 +# # ... +# # ... +# # the same with that using kaldi commands: compute-fbank-feats --dither=0.0 + + +# mfcc_feat = mfcc(sig,dither=0,useEnergy=True,wintype='povey') + +# # the computed mfcc coefficents of english.wav with dimension [110,13] +# # [ 17.1337 -23.3651 -7.41751 -7.73686 -21.3682 -8.93884 -3.70843 4.68346 -16.0676 12.782 -7.24054 8.25089 10.7292 +# # 17.1692 -23.3028 -5.61872 -4.0075 -23.287 -20.6101 -5.51584 -6.15273 -14.4333 8.13052 -0.0345329 2.06274 -0.564298 +# # ... +# # ... +# # the same with that using kaldi commands: compute-mfcc-feats --dither=0.0 + + + +if __name__ == '__main__': + unittest.main() \ No newline at end of file