hoo2/Multimedia_AAC_Project
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Level 3: small fixes and report

Browse Source
This commit is contained in:
hoo2
2026-02-22 01:39:12 +02:00
parent 6b1ebc3cfe
commit 17228f8539
36 changed files with 4660 additions and 154 deletions
Download Patch File Download Diff File Expand all files Collapse all files
source/level_2/core/aac_coder.py
+2 -21
View File
@@ -111,21 +111,6 @@ def _thresholds_from_smr(
return T
def _normalize_global_gain(G: GlobalGain) -> float | FloatArray:
"""
Normalize GlobalGain to match AACChannelFrameF3["G"] type:
- long: return float
- ESH: return float64 ndarray of shape (1, 8)
"""
if np.isscalar(G):
return float(G)
G_arr = np.asarray(G)
if G_arr.size == 1:
return float(G_arr.reshape(-1)[0])
return np.asarray(G_arr, dtype=np.float64)
# -----------------------------------------------------------------------------
# Public helpers (useful for level_x demo wrappers)
# -----------------------------------------------------------------------------
@@ -476,10 +461,6 @@ def aac_coder_3(
S_L, sfc_L, G_L = aac_quantizer(chl_f_tns, frame_type, SMR_L)
S_R, sfc_R, G_R = aac_quantizer(chr_f_tns, frame_type, SMR_R)
# Normalize G types for AACSeq3 schema (float | float64 ndarray).
G_Ln = _normalize_global_gain(G_L)
G_Rn = _normalize_global_gain(G_R)
# Huffman-code ONLY the DPCM differences for b>0.
# sfc[0] corresponds to alpha(0)=G and is stored separately in the frame.
sfc_L_dpcm = np.asarray(sfc_L, dtype=np.int64)[1:, ...]
@@ -503,7 +484,7 @@ def aac_coder_3(
"chl": {
"tns_coeffs": np.asarray(chl_tns_coeffs, dtype=np.float64),
"T": np.asarray(T_L, dtype=np.float64),
"G": G_Ln,
"G": G_L,
"sfc": sfc_L_stream,
"stream": mdct_L_stream,
"codebook": int(cb_L),
@@ -511,7 +492,7 @@ def aac_coder_3(
"chr": {
"tns_coeffs": np.asarray(chr_tns_coeffs, dtype=np.float64),
"T": np.asarray(T_R, dtype=np.float64),
"G": G_Rn,
"G": G_R,
"sfc": sfc_R_stream,
"stream": mdct_R_stream,
"codebook": int(cb_R),
source/level_2/core/aac_decoder.py
+3 -3
View File
@@ -42,7 +42,7 @@ def _nbands(frame_type: FrameType) -> int:
# Public helpers
# -----------------------------------------------------------------------------
def aac_unpack_seq_channels_to_frame_f(frame_type: FrameType, chl_f: FrameChannelF, chr_f: FrameChannelF) -> FrameF:
def aac_unpack_seq_channels(frame_type: FrameType, chl_f: FrameChannelF, chr_f: FrameChannelF) -> FrameF:
"""
Re-pack per-channel spectra from the Level-1 AACSeq1 schema into the stereo
FrameF container expected by aac_i_filter_bank().
@@ -167,7 +167,7 @@ def aac_decoder_1(
chl_f = np.asarray(fr["chl"]["frame_F"], dtype=np.float64)
chr_f = np.asarray(fr["chr"]["frame_F"], dtype=np.float64)
frame_f: FrameF = aac_unpack_seq_channels_to_frame_f(frame_type, chl_f, chr_f)
frame_f: FrameF = aac_unpack_seq_channels(frame_type, chl_f, chr_f)
frame_t_hat: FrameT = aac_i_filter_bank(frame_f, frame_type, win_type) # (2048, 2)
start = i * hop
@@ -427,7 +427,7 @@ def aac_decoder_3(
X_R = aac_i_tns(Xq_R, frame_type, tns_R)
# Re-pack to stereo container and inverse filterbank
frame_f = aac_unpack_seq_channels_to_frame_f(frame_type, np.asarray(X_L), np.asarray(X_R))
frame_f = aac_unpack_seq_channels(frame_type, np.asarray(X_L), np.asarray(X_R))
frame_t_hat: FrameT = aac_i_filter_bank(frame_f, frame_type, win_type)
start = i * hop
source/level_2/core/aac_huffman.py
+112
View File
@@ -0,0 +1,112 @@
# ------------------------------------------------------------
# AAC Coder/Decoder - Huffman wrappers (Level 3)
#
# Multimedia course at Aristotle University of
# Thessaloniki (AUTh)
#
# Author:
# Christos Choutouridis (ΑΕΜ 8997)
# cchoutou@ece.auth.gr
#
# Description:
# Thin wrappers around the provided Huffman utilities (material/huff_utils.py)
# so that the API matches the assignment text.
#
# Exposed API (assignment):
# huff_sec, huff_codebook = aac_encode_huff(coeff_sec, huff_LUT_list, force_codebook)
# dec_coeffs = aac_decode_huff(huff_sec, huff_codebook, huff_LUT_list)
#
# Notes:
# - Huffman coding operates on tuples. Therefore, decode(encode(x)) may return
# extra trailing symbols due to tuple padding. The AAC decoder knows the
# true section length from side information (band limits) and truncates.
# ------------------------------------------------------------
from __future__ import annotations
from typing import Any
import numpy as np
from material.huff_utils import encode_huff, decode_huff
def aac_encode_huff(
coeff_sec: np.ndarray,
huff_LUT_list: list[dict[str, Any]],
force_codebook: int | None = None,
) -> tuple[str, int]:
"""
Huffman-encode a section of coefficients (MDCT symbols or scalefactors).
Parameters
----------
coeff_sec : np.ndarray
Coefficient section to be encoded. Any shape is accepted; the input
is flattened and treated as a 1-D sequence of int64 symbols.
huff_LUT_list : list[dict[str, Any]]
List of Huffman Look-Up Tables (LUTs) as returned by material.load_LUT().
Index corresponds to codebook id (typically 1..11, with 0 reserved).
force_codebook : int | None
If provided, forces the use of this Huffman codebook. In the assignment,
scalefactors are encoded with codebook 11. For MDCT coefficients, this
argument is usually omitted (auto-selection).
Returns
-------
tuple[str, int]
(huff_sec, huff_codebook)
- huff_sec: bitstream as a string of '0'/'1'
- huff_codebook: codebook id used by the encoder
"""
coeff_sec_arr = np.asarray(coeff_sec, dtype=np.int64).reshape(-1)
if force_codebook is None:
# Provided utility returns (bitstream, codebook) in the auto-selection case.
huff_sec, huff_codebook = encode_huff(coeff_sec_arr, huff_LUT_list)
return str(huff_sec), int(huff_codebook)
# Provided utility returns ONLY the bitstream when force_codebook is set.
cb = int(force_codebook)
huff_sec = encode_huff(coeff_sec_arr, huff_LUT_list, force_codebook=cb)
return str(huff_sec), cb
def aac_decode_huff(
huff_sec: str | np.ndarray,
huff_codebook: int,
huff_LUT: list[dict[str, Any]],
) -> np.ndarray:
"""
Huffman-decode a bitstream using the specified codebook.
Parameters
----------
huff_sec : str | np.ndarray
Huffman bitstream. Typically a string of '0'/'1'. If an array is provided,
it is passed through to the provided decoder.
huff_codebook : int
Codebook id that was returned by aac_encode_huff.
Codebook 0 represents an all-zero section.
huff_LUT : list[dict[str, Any]]
Huffman LUT list as returned by material.load_LUT().
Returns
-------
np.ndarray
Decoded coefficients as a 1-D np.int64 array.
Note: Due to tuple coding, the decoded array may contain extra trailing
padding symbols. The caller must truncate to the known section length.
"""
cb = int(huff_codebook)
if cb == 0:
# Codebook 0 represents an all-zero section. The decoded length is not
# recoverable from the bitstream alone; the caller must expand/truncate.
return np.zeros((0,), dtype=np.int64)
if cb < 0 or cb >= len(huff_LUT):
raise ValueError(f"Invalid Huffman codebook index: {cb}")
lut = huff_LUT[cb]
dec = decode_huff(huff_sec, lut)
return np.asarray(dec, dtype=np.int64).reshape(-1)
source/level_2/core/aac_psycho.py
+441
View File
@@ -0,0 +1,441 @@
# ------------------------------------------------------------
# AAC Coder/Decoder - Psychoacoustic Model
#
# Multimedia course at Aristotle University of
# Thessaloniki (AUTh)
#
# Author:
# Christos Choutouridis (ΑΕΜ 8997)
# cchoutou@ece.auth.gr
#
# Description:
# Psychoacoustic model for ONE channel, based on the assignment notes (Section 2.4).
#
# Public API:
# SMR = aac_psycho(frame_T, frame_type, frame_T_prev_1, frame_T_prev_2)
#
# Output:
# - For long frames ("OLS", "LSS", "LPS"): SMR has shape (69,)
# - For short frames ("ESH"): SMR has shape (42, 8) (one column per subframe)
#
# Notes:
# - Uses Bark band tables from material/TableB219.mat:
# * B219a for long windows (69 bands, N=2048 FFT, N/2=1024 bins)
# * B219b for short windows (42 bands, N=256 FFT, N/2=128 bins)
# - Applies a Hann window in time domain before FFT magnitude/phase extraction.
# - Implements:
# spreading function -> band spreading -> tonality index -> masking thresholds -> SMR.
# ------------------------------------------------------------
from __future__ import annotations
import numpy as np
from core.aac_utils import band_limits, get_table
from core.aac_configuration import NMT_DB, TMN_DB
from core.aac_types import *
# -----------------------------------------------------------------------------
# Spreading function
# -----------------------------------------------------------------------------
def _spreading_matrix(bval: BandValueArray) -> FloatArray:
"""
Compute the spreading function matrix between psychoacoustic bands.
The spreading function describes how energy in one critical band masks
nearby bands. The formula follows the assignment pseudo-code.
Parameters
----------
bval : BandValueArray
Bark value per band, shape (B,).
Returns
-------
FloatArray
Spreading matrix S of shape (B, B), where:
S[bb, b] quantifies the contribution of band bb masking band b.
"""
bval = np.asarray(bval, dtype=np.float64).reshape(-1)
B = int(bval.shape[0])
spread = np.zeros((B, B), dtype=np.float64)
for b in range(B):
for bb in range(B):
# tmpx depends on direction (asymmetric spreading)
if bb >= b:
tmpx = 3.0 * (bval[bb] - bval[b])
else:
tmpx = 1.5 * (bval[bb] - bval[b])
# tmpz uses the "min(..., 0)" nonlinearity exactly as in the notes
tmpz = 8.0 * min((tmpx - 0.5) ** 2 - 2.0 * (tmpx - 0.5), 0.0)
tmpy = 15.811389 + 7.5 * (tmpx + 0.474) - 17.5 * np.sqrt(1.0 + (tmpx + 0.474) ** 2)
# Clamp very small values (below -100 dB) to 0 contribution
if tmpy < -100.0:
spread[bb, b] = 0.0
else:
spread[bb, b] = 10.0 ** ((tmpz + tmpy) / 10.0)
return spread
# -----------------------------------------------------------------------------
# Windowing + FFT feature extraction
# -----------------------------------------------------------------------------
def _hann_window(N: int) -> FloatArray:
"""
Hann window as specified in the notes:
w[n] = 0.5 - 0.5*cos(2*pi*(n + 0.5)/N)
Parameters
----------
N : int
Window length.
Returns
-------
FloatArray
1-D array of shape (N,), dtype float64.
"""
n = np.arange(N, dtype=np.float64)
return 0.5 - 0.5 * np.cos((2.0 * np.pi / N) * (n + 0.5))
def _r_phi_from_time(x: FrameChannelT, N: int) -> tuple[FloatArray, FloatArray]:
"""
Compute FFT magnitude r(w) and phase phi(w) for bins w = 0 .. N/2-1.
Processing:
1) Apply Hann window in time domain.
2) Compute N-point FFT.
3) Keep only the positive-frequency bins [0 .. N/2-1].
Parameters
----------
x : FrameChannelT
Time-domain samples, shape (N,).
N : int
FFT size (2048 or 256).
Returns
-------
r : FloatArray
Magnitude spectrum for bins 0 .. N/2-1, shape (N/2,).
phi : FloatArray
Phase spectrum for bins 0 .. N/2-1, shape (N/2,).
"""
x = np.asarray(x, dtype=np.float64).reshape(-1)
if x.shape[0] != N:
raise ValueError(f"Expected time vector of length {N}, got {x.shape[0]}.")
w = _hann_window(N)
X = np.fft.fft(x * w, n=N)
Xp = X[: N // 2]
r = np.abs(Xp).astype(np.float64, copy=False)
phi = np.angle(Xp).astype(np.float64, copy=False)
return r, phi
def _predictability(
r: FloatArray,
phi: FloatArray,
r_m1: FloatArray,
phi_m1: FloatArray,
r_m2: FloatArray,
phi_m2: FloatArray,
) -> FloatArray:
"""
Compute predictability c(w) per spectral bin.
The notes define:
r_pred(w) = 2*r_{-1}(w) - r_{-2}(w)
phi_pred(w) = 2*phi_{-1}(w) - phi_{-2}(w)
c(w) = |X(w) - X_pred(w)| / (r(w) + |r_pred(w)|)
where X(w) is represented in polar form using r(w), phi(w).
Parameters
----------
r, phi : FloatArray
Current magnitude and phase, shape (N/2,).
r_m1, phi_m1 : FloatArray
Previous magnitude and phase, shape (N/2,).
r_m2, phi_m2 : FloatArray
Pre-previous magnitude and phase, shape (N/2,).
Returns
-------
FloatArray
Predictability c(w), shape (N/2,).
"""
r_pred = 2.0 * r_m1 - r_m2
phi_pred = 2.0 * phi_m1 - phi_m2
num = np.sqrt(
(r * np.cos(phi) - r_pred * np.cos(phi_pred)) ** 2
+ (r * np.sin(phi) - r_pred * np.sin(phi_pred)) ** 2
)
den = r + np.abs(r_pred) + 1e-12 # avoid division-by-zero without altering behavior
return (num / den).astype(np.float64, copy=False)
# -----------------------------------------------------------------------------
# Band-domain aggregation
# -----------------------------------------------------------------------------
def _band_energy_and_pred(
r: FloatArray,
c: FloatArray,
wlow: BandIndexArray,
whigh: BandIndexArray,
) -> tuple[FloatArray, FloatArray]:
"""
Aggregate spectral bin quantities into psychoacoustic bands.
Definitions (notes):
e(b) = sum_{w=wlow(b)..whigh(b)} r(w)^2
c_num(b) = sum_{w=wlow(b)..whigh(b)} c(w) * r(w)^2
The band predictability c(b) is later computed after spreading as:
cb(b) = ct(b) / ecb(b)
Parameters
----------
r : FloatArray
Magnitude spectrum, shape (N/2,).
c : FloatArray
Predictability per bin, shape (N/2,).
wlow, whigh : BandIndexArray
Band limits (inclusive indices), shape (B,).
Returns
-------
e_b : FloatArray
Band energies e(b), shape (B,).
c_num_b : FloatArray
Weighted predictability numerators c_num(b), shape (B,).
"""
r2 = (r * r).astype(np.float64, copy=False)
B = int(wlow.shape[0])
e_b = np.zeros(B, dtype=np.float64)
c_num_b = np.zeros(B, dtype=np.float64)
for b in range(B):
a = int(wlow[b])
z = int(whigh[b])
seg_r2 = r2[a : z + 1]
e_b[b] = float(np.sum(seg_r2))
c_num_b[b] = float(np.sum(c[a : z + 1] * seg_r2))
return e_b, c_num_b
def _psycho_window(
time_x: FrameChannelT,
prev1_x: FrameChannelT,
prev2_x: FrameChannelT,
*,
N: int,
table: BarkTable,
) -> FloatArray:
"""
Compute SMR for one FFT analysis window (N=2048 for long, N=256 for short).
This implements the pipeline described in the notes:
- FFT magnitude/phase
- predictability per bin
- band energies and predictability
- band spreading
- tonality index tb(b)
- masking threshold (noise + threshold in quiet)
- SMR(b) = e(b) / np(b)
Parameters
----------
time_x : FrameChannelT
Current time-domain samples, shape (N,).
prev1_x : FrameChannelT
Previous time-domain samples, shape (N,).
prev2_x : FrameChannelT
Pre-previous time-domain samples, shape (N,).
N : int
FFT size.
table : BarkTable
Psychoacoustic band table (B219a or B219b).
Returns
-------
FloatArray
SMR per band, shape (B,).
"""
wlow, whigh, bval, qthr_db = band_limits(table)
spread = _spreading_matrix(bval)
# FFT features for current and history windows
r, phi = _r_phi_from_time(time_x, N)
r_m1, phi_m1 = _r_phi_from_time(prev1_x, N)
r_m2, phi_m2 = _r_phi_from_time(prev2_x, N)
# Predictability per bin
c_w = _predictability(r, phi, r_m1, phi_m1, r_m2, phi_m2)
# Aggregate into psycho bands
e_b, c_num_b = _band_energy_and_pred(r, c_w, wlow, whigh)
# Spread energies and predictability across bands:
# ecb(b) = sum_bb e(bb) * S(bb, b)
# ct(b) = sum_bb c_num(bb) * S(bb, b)
ecb = spread.T @ e_b
ct = spread.T @ c_num_b
# Band predictability after spreading: cb(b) = ct(b) / ecb(b)
cb = ct / (ecb + 1e-12)
# Normalized energy term:
# en(b) = ecb(b) / sum_bb S(bb, b)
spread_colsum = np.sum(spread, axis=0)
en = ecb / (spread_colsum + 1e-12)
# Tonality index (clamped to [0, 1])
tb = -0.299 - 0.43 * np.log(np.maximum(cb, 1e-12))
tb = np.clip(tb, 0.0, 1.0)
# Required SNR per band (dB): interpolate between TMN and NMT
snr_b = tb * TMN_DB + (1.0 - tb) * NMT_DB
bc = 10.0 ** (-snr_b / 10.0)
# Noise masking threshold estimate (power domain)
nb = en * bc
# Threshold in quiet (convert from dB to power domain):
# qthr_power = eps * (N/2) * 10^(qthr_db/10)
qthr_power = np.finfo('float').eps * (N / 2.0) * (10.0 ** (qthr_db / 10.0))
# Final masking threshold per band:
# np(b) = max(nb(b), qthr(b))
npart = np.maximum(nb, qthr_power)
# Signal-to-mask ratio:
# SMR(b) = e(b) / np(b)
smr = e_b / (npart + 1e-12)
return smr.astype(np.float64, copy=False)
# -----------------------------------------------------------------------------
# ESH window slicing (match filterbank conventions)
# -----------------------------------------------------------------------------
def _esh_subframes(x_2048: FrameChannelT) -> list[FrameChannelT]:
"""
Extract the 8 overlapping 256-sample short windows used by AAC ESH.
The project convention (matching the filterbank) is:
start_j = 448 + 128*j, for j = 0..7
subframe_j = x[start_j : start_j + 256]
This selects the central 1152-sample region [448, 1600) and produces
8 windows with 50% overlap.
Parameters
----------
x_2048 : FrameChannelT
Time-domain channel frame, shape (2048,).
Returns
-------
list[FrameChannelT]
List of 8 subframes, each of shape (256,).
"""
x_2048 = np.asarray(x_2048, dtype=np.float64).reshape(-1)
if x_2048.shape[0] != 2048:
raise ValueError("ESH requires 2048-sample input frames.")
subs: list[FrameChannelT] = []
for j in range(8):
start = 448 + 128 * j
subs.append(x_2048[start : start + 256])
return subs
# -----------------------------------------------------------------------------
# Public API
# -----------------------------------------------------------------------------
def aac_psycho(
frame_T: FrameChannelT,
frame_type: FrameType,
frame_T_prev_1: FrameChannelT,
frame_T_prev_2: FrameChannelT,
) -> FloatArray:
"""
Psychoacoustic model for ONE channel.
Parameters
----------
frame_T : FrameChannelT
Current time-domain channel frame, shape (2048,).
For "ESH", the 8 short windows are derived internally.
frame_type : FrameType
AAC frame type ("OLS", "LSS", "ESH", "LPS").
frame_T_prev_1 : FrameChannelT
Previous time-domain channel frame, shape (2048,).
frame_T_prev_2 : FrameChannelT
Pre-previous time-domain channel frame, shape (2048,).
Returns
-------
FloatArray
Signal-to-Mask Ratio (SMR), per psychoacoustic band.
- If frame_type == "ESH": shape (42, 8)
- Else: shape (69,)
"""
frame_T = np.asarray(frame_T, dtype=np.float64).reshape(-1)
frame_T_prev_1 = np.asarray(frame_T_prev_1, dtype=np.float64).reshape(-1)
frame_T_prev_2 = np.asarray(frame_T_prev_2, dtype=np.float64).reshape(-1)
if frame_T.shape[0] != 2048 or frame_T_prev_1.shape[0] != 2048 or frame_T_prev_2.shape[0] != 2048:
raise ValueError("aac_psycho expects 2048-sample frames for current/prev1/prev2.")
table, N = get_table(frame_type)
# Long frame types: compute one SMR vector (69 bands)
if frame_type != "ESH":
return _psycho_window(frame_T, frame_T_prev_1, frame_T_prev_2, N=N, table=table)
# ESH: compute 8 SMR vectors (42 bands each), one per short subframe.
#
# The notes use short-window history for predictability:
# - For j=0: use previous frame's subframes (7, 6)
# - For j=1: use current subframe 0 and previous frame's subframe 7
# - For j>=2: use current subframes (j-1, j-2)
#
# This matches the "within-frame history" convention commonly used in
# simplified psycho models for ESH.
cur_subs = _esh_subframes(frame_T)
prev1_subs = _esh_subframes(frame_T_prev_1)
B = int(table.shape[0]) # expected 42
smr_out = np.zeros((B, 8), dtype=np.float64)
for j in range(8):
if j == 0:
x_m1 = prev1_subs[7]
x_m2 = prev1_subs[6]
elif j == 1:
x_m1 = cur_subs[0]
x_m2 = prev1_subs[7]
else:
x_m1 = cur_subs[j - 1]
x_m2 = cur_subs[j - 2]
smr_out[:, j] = _psycho_window(cur_subs[j], x_m1, x_m2, N=256, table=table)
return smr_out
source/level_2/core/aac_quantizer.py
+600
View File
@@ -0,0 +1,600 @@
# ------------------------------------------------------------
# AAC Coder/Decoder - Quantizer / iQuantizer (Level 3)
#
# Multimedia course at Aristotle University of
# Thessaloniki (AUTh)
#
# Author:
# Christos Choutouridis (ΑΕΜ 8997)
# cchoutou@ece.auth.gr
#
# Description:
# Implements AAC quantizer and inverse quantizer for one channel.
# Based on assignment section 2.6 (Eq. 12-15).
#
# Notes:
# - Bit reservoir is not implemented (assignment simplification).
# - Scalefactor bands are assumed equal to psychoacoustic bands
# (Table B.2.1.9a / B.2.1.9b from TableB219.mat).
# ------------------------------------------------------------
from __future__ import annotations
import numpy as np
from core.aac_utils import get_table, band_limits
from core.aac_types import *
# -----------------------------------------------------------------------------
# Constants (assignment)
# -----------------------------------------------------------------------------
MAGIC_NUMBER: float = 0.4054
EPS: float = 1e-12
MAX_SF_DELTA:int = 60
# -----------------------------------------------------------------------------
# Helpers: ESH packing/unpacking (128x8 <-> 1024x1)
# -----------------------------------------------------------------------------
def _esh_pack(x_128x8: FloatArray) -> FloatArray:
"""
Pack ESH coefficients (128 x 8) into a single long vector (1024 x 1).
Packing order:
Columns are concatenated in subframe order (0..7), column-major.
Parameters
----------
x_128x8 : FloatArray
ESH coefficients, shape (128, 8).
Returns
-------
FloatArray
Packed coefficients, shape (1024, 1).
"""
x_128x8 = np.asarray(x_128x8, dtype=np.float64)
if x_128x8.shape != (128, 8):
raise ValueError("ESH pack expects shape (128, 8).")
return x_128x8.reshape(1024, 1, order="F")
def _esh_unpack(x_1024x1: FloatArray) -> FloatArray:
"""
Unpack a packed ESH vector (1024 elements) back to shape (128, 8).
Parameters
----------
x_1024x1 : FloatArray
Packed ESH vector, shape (1024,) or (1024, 1) after flattening.
Returns
-------
FloatArray
Unpacked ESH coefficients, shape (128, 8).
"""
x_1024x1 = np.asarray(x_1024x1, dtype=np.float64).reshape(-1)
if x_1024x1.shape[0] != 1024:
raise ValueError("ESH unpack expects 1024 elements.")
return x_1024x1.reshape(128, 8, order="F")
# -----------------------------------------------------------------------------
# Core quantizer formulas (Eq. 12, Eq. 13)
# -----------------------------------------------------------------------------
def _quantize_symbol(x: FloatArray, alpha: float) -> QuantizedSymbols:
"""
Quantize MDCT coefficients to integer symbols S(k).
Implements Eq. (12):
S(k) = sgn(X(k)) * int( (|X(k)| * 2^(-alpha/4))^(3/4) + MAGIC_NUMBER )
Parameters
----------
x : FloatArray
MDCT coefficients for a contiguous set of spectral lines.
Shape: (N,)
alpha : float
Scalefactor gain for the corresponding scalefactor band.
Returns
-------
QuantizedSymbols
Quantized symbols S(k) as int64, shape (N,).
"""
x = np.asarray(x, dtype=np.float64)
scale = 2.0 ** (-0.25 * float(alpha))
ax = np.abs(x) * scale
y = np.power(ax, 0.75, dtype=np.float64)
# "int" in the assignment corresponds to truncation.
q = np.floor(y + MAGIC_NUMBER).astype(np.int64)
return (np.sign(x).astype(np.int64) * q).astype(np.int64)
def _dequantize_symbol(S: QuantizedSymbols, alpha: float) -> FloatArray:
"""
Inverse quantizer (dequantization of symbols).
Implements Eq. (13):
Xhat(k) = sgn(S(k)) * |S(k)|^(4/3) * 2^(alpha/4)
Parameters
----------
S : QuantizedSymbols
Quantized symbols S(k), int64, shape (N,).
alpha : float
Scalefactor gain for the corresponding scalefactor band.
Returns
-------
FloatArray
Reconstructed MDCT coefficients Xhat(k), float64, shape (N,).
"""
S = np.asarray(S, dtype=np.int64)
scale = 2.0 ** (0.25 * float(alpha))
aS = np.abs(S).astype(np.float64)
y = np.power(aS, 4.0 / 3.0, dtype=np.float64)
return (np.sign(S).astype(np.float64) * y * scale).astype(np.float64)
# -----------------------------------------------------------------------------
# Alpha initialization (Eq. 14)
# -----------------------------------------------------------------------------
def _initial_alpha_hat(X: "FloatArray", MQ: int = 8191) -> int:
"""
Compute the initial scalefactor estimate alpha_hat for a frame.
The assignment proposes the following first approximation (Equation 14):
alpha_hat = (16/3) * log2( max_k(|X(k)|)^(3/4) / MQ )
where max_k runs over all MDCT coefficients of the frame (not per band),
and MQ is the maximum quantization level parameter (2*MQ + 1 levels).
Parameters
----------
X : FloatArray
MDCT coefficients of one frame (or one ESH subframe), shape (N,).
MQ : int
Quantizer parameter (default 8191, as per assignment).
Returns
-------
int
Integer alpha_hat (rounded to nearest integer).
"""
x_max = float(np.max(np.abs(X)))
if x_max <= 0.0:
return 0
alpha_hat = (16.0 / 3.0) * np.log2((x_max ** (3.0 / 4.0)) / float(MQ))
return int(np.round(alpha_hat))
# -----------------------------------------------------------------------------
# Band utilities
# -----------------------------------------------------------------------------
def _band_slices(frame_type: FrameType) -> list[tuple[int, int]]:
"""
Return scalefactor band ranges [wlow, whigh] (inclusive) for the given frame type.
These are derived from the psychoacoustic tables (TableB219),
and map directly to MDCT indices:
- long: 0..1023
- short (ESH subframe): 0..127
Parameters
----------
frame_type : FrameType
Frame type ("OLS", "LSS", "ESH", "LPS").
Returns
-------
list[tuple[int, int]]
List of (lo, hi) inclusive index pairs for each band.
"""
table, _Nfft = get_table(frame_type)
wlow, whigh, _bval, _qthr_db = band_limits(table)
bands: list[tuple[int, int]] = []
for lo, hi in zip(wlow, whigh):
bands.append((int(lo), int(hi)))
return bands
def _band_energy(x: FloatArray, lo: int, hi: int) -> float:
"""
Compute energy of a spectral segment x[lo:hi+1].
Parameters
----------
x : FloatArray
MDCT coefficient vector.
lo, hi : int
Inclusive index range.
Returns
-------
float
Sum of squares (energy) within the band.
"""
sec = x[lo : hi + 1]
return float(np.sum(sec * sec))
def _psychoacoustic_threshold(
X: FloatArray,
SMR_col: FloatArray,
bands: list[tuple[int, int]],
) -> FloatArray:
"""
Compute psychoacoustic thresholds T(b) per band.
Uses:
P(b) = sum_{k in band} X(k)^2
T(b) = P(b) / SMR(b)
Parameters
----------
X : FloatArray
MDCT coefficients for a frame (long) or one ESH subframe (short).
SMR_col : FloatArray
SMR values for this frame/subframe, shape (NB,).
bands : list[tuple[int, int]]
Band index ranges.
Returns
-------
FloatArray
Threshold vector T(b), shape (NB,).
"""
nb = len(bands)
T = np.zeros((nb,), dtype=np.float64)
for b, (lo, hi) in enumerate(bands):
P = _band_energy(X, lo, hi)
smr = float(SMR_col[b])
if smr <= EPS:
T[b] = 0.0
else:
T[b] = P / smr
return T
# -----------------------------------------------------------------------------
# Alpha selection per band + neighbor-difference constraint
# -----------------------------------------------------------------------------
def _best_alpha_for_band(
X: "FloatArray", lo: int, hi: int, T_b: float,
alpha_hat: int, alpha_prev: int, alpha_min: int, alpha_max: int,
) -> int:
"""
Determine the band-wise scalefactor alpha(b) following the assignment.
Procedure:
- Start from a frame-wise initial estimate alpha_hat.
- Iteratively increase alpha(b) by 1 as long as the quantization error power
stays below the psychoacoustic threshold T(b): P_e(b) = sum_{k in band} ( X(k) - Xhat(k) )^2
- Stop increasing alpha(b) if the neighbor constraint would be violated: |alpha(b) - alpha(b-1)| <= 60
When processing bands sequentially (low -> high), this becomes: alpha(b) <= alpha_prev + 60
Notes:
- This function does not decrease alpha if the initial value already violates
the threshold; the assignment only specifies iterative increase.
Parameters
----------
X : FloatArray
Full MDCT vector of the current (sub)frame, shape (N,).
lo, hi : int
Band index bounds (inclusive), defining the band slice.
T_b : float
Threshold T(b) for this band.
alpha_hat : int
Initial frame-wise estimate (Equation 14).
alpha_prev : int
Previously selected alpha for band b-1 (neighbor constraint reference).
alpha_min, alpha_max : int
Safeguard bounds for alpha.
Returns
-------
int
Selected integer alpha(b).
"""
if T_b <= 0.0:
return int(alpha_hat)
Xsec = X[lo : hi + 1]
# Neighbor constraint (sequential processing): alpha(b) <= alpha_prev + 60
alpha_limit = min(int(alpha_max), int(alpha_prev) + MAX_SF_DELTA)
# Start from alpha_hat, clamped to feasible range
alpha = int(alpha_hat)
alpha = max(int(alpha_min), min(alpha, int(alpha_limit)))
# Evaluate at current alpha
Ssec = _quantize_symbol(Xsec, alpha)
Xhat = _dequantize_symbol(Ssec, alpha)
Pe = float(np.sum((Xsec - Xhat) ** 2))
# If already above threshold, return current alpha (no decrease step specified)
if Pe > T_b:
return alpha
# Increase alpha while still under threshold and within constraints
while True:
alpha_next = alpha + 1
if alpha_next > alpha_limit:
break
Ssec = _quantize_symbol(Xsec, alpha_next)
Xhat = _dequantize_symbol(Ssec, alpha_next)
Pe_next = float(np.sum((Xsec - Xhat) ** 2))
if Pe_next > T_b:
break
alpha = alpha_next
return alpha
# -----------------------------------------------------------------------------
# Public API
# -----------------------------------------------------------------------------
def aac_quantizer(
frame_F: FrameChannelF,
frame_type: FrameType,
SMR: FloatArray,
) -> tuple[QuantizedSymbols, ScaleFactors, GlobalGain]:
"""
AAC quantizer for one channel (Level 3).
Quantizes MDCT coefficients (after TNS) using band-wise scalefactors derived
from psychoacoustic thresholds computed via SMR.
The implementation follows the assignment procedure:
- Compute an initial frame-wise alpha_hat using Equation (14), based on the
maximum MDCT coefficient magnitude of the (sub)frame.
- For each band b, increase alpha(b) by 1 while the quantization error power
P_e(b) stays below the threshold T(b).
- Enforce the neighbor constraint |alpha(b) - alpha(b-1)| <= 60 during the
band-by-band search (no post-processing needed).
Parameters
----------
frame_F : FrameChannelF
MDCT coefficients after TNS, one channel.
Shapes:
- Long frames: (1024,) or (1024, 1)
- ESH: (128, 8)
frame_type : FrameType
AAC frame type ("OLS", "LSS", "ESH", "LPS").
SMR : FloatArray
Signal-to-Mask Ratio per band.
Shapes:
- Long: (NB,) or (NB, 1)
- ESH: (NB, 8)
Returns
-------
S : QuantizedSymbols
Quantized symbols S(k), packed as shape (1024, 1) for all frame types.
For ESH, the 8 subframes are packed in column-major subframe layout.
sfc : ScaleFactors
DPCM-coded scalefactors:
sfc(0) = alpha(0) = G
sfc(b) = alpha(b) - alpha(b-1), for b > 0
Shapes:
- Long: (NB, 1)
- ESH: (NB, 8)
G : GlobalGain
Global gain G = alpha(0).
- Long: scalar float
- ESH: array shape (1, 8), dtype float64
"""
bands = _band_slices(frame_type)
NB = len(bands)
X = np.asarray(frame_F, dtype=np.float64)
SMR = np.asarray(SMR, dtype=np.float64)
# -------------------------------------------------------------------------
# ESH: 8 short subframes, each of length 128
# -------------------------------------------------------------------------
if frame_type == "ESH":
if X.shape != (128, 8):
raise ValueError("For ESH, frame_F must have shape (128, 8).")
if SMR.shape != (NB, 8):
raise ValueError(f"For ESH, SMR must have shape ({NB}, 8).")
S_out: QuantizedSymbols = np.zeros((1024, 1), dtype=np.int64)
sfc: ScaleFactors = np.zeros((NB, 8), dtype=np.int64)
G_arr = np.zeros((1, 8), dtype=np.float64)
# Packed output view: (128, 8) with column-major layout
S_pack = S_out[:, 0].reshape(128, 8, order="F")
for j in range(8):
Xj = X[:, j].reshape(128)
SMRj = SMR[:, j].reshape(NB)
# Compute psychoacoustic threshold T(b) for this subframe
T = _psychoacoustic_threshold(Xj, SMRj, bands)
# Frame-wise initial estimate alpha_hat (Equation 14)
alpha_hat = _initial_alpha_hat(Xj)
# Band-wise scalefactors alpha(b)
alpha = np.zeros((NB,), dtype=np.int64)
alpha_prev = int(alpha_hat)
for b, (lo, hi) in enumerate(bands):
alpha_b = _best_alpha_for_band(
X=Xj,
lo=lo,
hi=hi,
T_b=float(T[b]),
alpha_hat=int(alpha_hat),
alpha_prev=int(alpha_prev),
alpha_min=-4096,
alpha_max=4096,
)
alpha[b] = int(alpha_b)
alpha_prev = int(alpha_b)
# DPCM-coded scalefactors
G_arr[0, j] = float(alpha[0])
sfc[0, j] = int(alpha[0])
for b in range(1, NB):
sfc[b, j] = int(alpha[b] - alpha[b - 1])
# Quantize MDCT coefficients band-by-band
Sj = np.zeros((128,), dtype=np.int64)
for b, (lo, hi) in enumerate(bands):
Sj[lo : hi + 1] = _quantize_symbol(Xj[lo : hi + 1], float(alpha[b]))
# Store subframe in packed output
S_pack[:, j] = Sj
return S_out, sfc, G_arr
# -------------------------------------------------------------------------
# Long frames: OLS / LSS / LPS, length 1024
# -------------------------------------------------------------------------
if X.shape == (1024,):
Xv = X
elif X.shape == (1024, 1):
Xv = X[:, 0]
else:
raise ValueError("For non-ESH, frame_F must have shape (1024,) or (1024, 1).")
if SMR.shape == (NB,):
SMRv = SMR
elif SMR.shape == (NB, 1):
SMRv = SMR[:, 0]
else:
raise ValueError(f"For non-ESH, SMR must have shape ({NB},) or ({NB}, 1).")
# Compute psychoacoustic threshold T(b) for the long frame
T = _psychoacoustic_threshold(Xv, SMRv, bands)
# Frame-wise initial estimate alpha_hat (Equation 14)
alpha_hat = _initial_alpha_hat(Xv)
# Band-wise scalefactors alpha(b)
alpha = np.zeros((NB,), dtype=np.int64)
alpha_prev = int(alpha_hat)
for b, (lo, hi) in enumerate(bands):
alpha_b = _best_alpha_for_band(
X=Xv,
lo=lo,
hi=hi,
T_b=float(T[b]),
alpha_hat=int(alpha_hat),
alpha_prev=int(alpha_prev),
alpha_min=-4096,
alpha_max=4096,
)
alpha[b] = int(alpha_b)
alpha_prev = int(alpha_b)
# DPCM-coded scalefactors
sfc_out: ScaleFactors = np.zeros((NB, 1), dtype=np.int64)
sfc_out[0, 0] = int(alpha[0])
for b in range(1, NB):
sfc_out[b, 0] = int(alpha[b] - alpha[b - 1])
G: float = float(alpha[0])
# Quantize MDCT coefficients band-by-band
S_vec = np.zeros((1024,), dtype=np.int64)
for b, (lo, hi) in enumerate(bands):
S_vec[lo : hi + 1] = _quantize_symbol(Xv[lo : hi + 1], float(alpha[b]))
return S_vec.reshape(1024, 1), sfc_out, G
def aac_i_quantizer(
S: QuantizedSymbols,
sfc: ScaleFactors,
G: GlobalGain,
frame_type: FrameType,
) -> FrameChannelF:
"""
Inverse quantizer (iQuantizer) for one channel.
Reconstructs MDCT coefficients from quantized symbols and DPCM scalefactors.
Parameters
----------
S : QuantizedSymbols
Quantized symbols, shape (1024, 1) (or any array with 1024 elements).
sfc : ScaleFactors
DPCM-coded scalefactors.
Shapes:
- Long: (NB, 1)
- ESH: (NB, 8)
G : GlobalGain
Global gain (not strictly required if sfc includes sfc(0)=alpha(0)).
Present for API compatibility with the assignment.
frame_type : FrameType
AAC frame type.
Returns
-------
FrameChannelF
Reconstructed MDCT coefficients:
- ESH: (128, 8)
- Long: (1024, 1)
"""
bands = _band_slices(frame_type)
NB = len(bands)
S_flat = np.asarray(S, dtype=np.int64).reshape(-1)
if S_flat.shape[0] != 1024:
raise ValueError("S must contain 1024 symbols.")
if frame_type == "ESH":
sfc = np.asarray(sfc, dtype=np.int64)
if sfc.shape != (NB, 8):
raise ValueError(f"For ESH, sfc must have shape ({NB}, 8).")
S_128x8 = _esh_unpack(S_flat)
Xrec = np.zeros((128, 8), dtype=np.float64)
for j in range(8):
alpha = np.zeros((NB,), dtype=np.int64)
alpha[0] = int(sfc[0, j])
for b in range(1, NB):
alpha[b] = int(alpha[b - 1] + sfc[b, j])
Xj = np.zeros((128,), dtype=np.float64)
for b, (lo, hi) in enumerate(bands):
Xj[lo : hi + 1] = _dequantize_symbol(S_128x8[lo : hi + 1, j].astype(np.int64), float(alpha[b]))
Xrec[:, j] = Xj
return Xrec
sfc = np.asarray(sfc, dtype=np.int64)
if sfc.shape != (NB, 1):
raise ValueError(f"For non-ESH, sfc must have shape ({NB}, 1).")
alpha = np.zeros((NB,), dtype=np.int64)
alpha[0] = int(sfc[0, 0])
for b in range(1, NB):
alpha[b] = int(alpha[b - 1] + sfc[b, 0])
Xrec = np.zeros((1024,), dtype=np.float64)
for b, (lo, hi) in enumerate(bands):
Xrec[lo : hi + 1] = _dequantize_symbol(S_flat[lo : hi + 1], float(alpha[b]))
return Xrec.reshape(1024, 1)
source/level_2/core/aac_tns.py
+6 -6
View File
@@ -39,7 +39,7 @@ from core.aac_types import *
# -----------------------------------------------------------------------------
def _band_ranges_for_kcount(k_count: int) -> BandRanges:
def _band_ranges(k_count: int) -> BandRanges:
"""
Return Bark band index ranges [start, end] (inclusive) for the given MDCT line count.
@@ -66,7 +66,7 @@ def _band_ranges_for_kcount(k_count: int) -> BandRanges:
start = tbl[:, 1].astype(int)
end = tbl[:, 2].astype(int)
ranges: list[tuple[int, int]] = [(int(s), int(e)) for s, e in zip(start, end)]
ranges: BandRanges = [(int(s), int(e)) for s, e in zip(start, end)]
for s, e in ranges:
if s < 0 or e < s or e >= k_count:
@@ -117,7 +117,7 @@ def _compute_sw(x: MdctCoeffs) -> MdctCoeffs:
x = np.asarray(x, dtype=np.float64).reshape(-1)
k_count = int(x.shape[0])
bands = _band_ranges_for_kcount(k_count)
bands = _band_ranges(k_count)
sw = np.zeros(k_count, dtype=np.float64)
for s, e in bands:
@@ -347,7 +347,7 @@ def _apply_itns_iir(y: MdctCoeffs, a_q: MdctCoeffs) -> MdctCoeffs:
return x_hat
def _tns_one_vector(x: MdctCoeffs) -> tuple[MdctCoeffs, MdctCoeffs]:
def _tns_vector(x: MdctCoeffs) -> tuple[MdctCoeffs, MdctCoeffs]:
"""
TNS for a single MDCT vector (one long frame or one short subframe).
@@ -430,7 +430,7 @@ def aac_tns(frame_F_in: FrameChannelF, frame_type: FrameType) -> Tuple[FrameChan
a_out = np.empty((PRED_ORDER, 8), dtype=np.float64)
for j in range(8):
y[:, j], a_out[:, j] = _tns_one_vector(x[:, j])
y[:, j], a_out[:, j] = _tns_vector(x[:, j])
return y, a_out
@@ -443,7 +443,7 @@ def aac_tns(frame_F_in: FrameChannelF, frame_type: FrameType) -> Tuple[FrameChan
else:
raise ValueError('For non-ESH, frame_F_in must have shape (1024,) or (1024, 1).')
y_vec, a_q = _tns_one_vector(x_vec)
y_vec, a_q = _tns_vector(x_vec)
if out_shape == (1024,):
y_out = y_vec
source/level_2/core/aac_utils.py
+36
View File
@@ -163,6 +163,42 @@ def snr_db(x_ref: StereoSignal, x_hat: StereoSignal) -> float:
return float(10.0 * np.log10(ps / pn))
def estimate_lag_mono(x_ref: TimeSignal, x_hat: TimeSignal, max_lag=4096):
"""
Estimate time lag between two mono signals.
Returns lag (positive means x_hat delayed).
"""
n = min(len(x_ref), len(x_hat))
x_ref = x_ref[:n]
x_hat = x_hat[:n]
corr = np.correlate(x_ref, x_hat, mode='full')
lags = np.arange(-n + 1, n)
center = n - 1
lo = max(0, center - max_lag)
hi = min(len(corr), center + max_lag + 1)
best = lo + int(np.argmax(corr[lo:hi]))
return int(lags[best])
def match_gain(x_ref: StereoSignal, x_hat: StereoSignal) -> float:
"""
Least-squares gain g that best maps x_hat -> x_ref.
"""
n = min(x_ref.shape[0], x_hat.shape[0])
c = min(x_ref.shape[1], x_hat.shape[1])
r = x_ref[:n, :c].reshape(-1).astype(np.float64)
h = x_hat[:n, :c].reshape(-1).astype(np.float64)
denom = float(np.dot(h, h))
if denom <= 0.0:
return 1.0
return float(np.dot(r, h) / denom)
# -----------------------------------------------------------------------------
# Psychoacoustic band tables (TableB219.mat)
# -----------------------------------------------------------------------------
source/level_2/material/huff_utils.py
+403
View File
@@ -0,0 +1,403 @@
import numpy as np
import scipy.io as sio
import os
# ------------------ LOAD LUT ------------------
def load_LUT(mat_filename=None):
"""
Loads the list of Huffman Codebooks (LUTs)
Returns:
huffLUT : list (index 1..11 used, index 0 unused)
"""
if mat_filename is None:
current_dir = os.path.dirname(os.path.abspath(__file__))
mat_filename = os.path.join(current_dir, "huffCodebooks.mat")
mat = sio.loadmat(mat_filename)
huffCodebooks_raw = mat['huffCodebooks'].squeeze()
huffCodebooks = []
for i in range(11):
huffCodebooks.append(np.array(huffCodebooks_raw[i]))
# Build inverse VLC tables
invTable = [None] * 11
for i in range(11):
h = huffCodebooks[i][:, 2].astype(int) # column 3
hlength = huffCodebooks[i][:, 1].astype(int) # column 2
hbin = []
for j in range(len(h)):
hbin.append(format(h[j], f'0{hlength[j]}b'))
invTable[i] = vlc_table(hbin)
# Build Huffman LUT dicts
huffLUT = [None] * 12 # index 0 unused
params = [
(4, 1, True),
(4, 1, True),
(4, 2, False),
(4, 2, False),
(2, 4, True),
(2, 4, True),
(2, 7, False),
(2, 7, False),
(2, 12, False),
(2, 12, False),
(2, 16, False),
]
for i, (nTupleSize, maxAbs, signed) in enumerate(params, start=1):
huffLUT[i] = {
'LUT': huffCodebooks[i-1],
'invTable': invTable[i-1],
'codebook': i,
'nTupleSize': nTupleSize,
'maxAbsCodeVal': maxAbs,
'signedValues': signed
}
return huffLUT
def vlc_table(code_array):
"""
codeArray: list of strings, each string is a Huffman codeword (e.g. '0101')
returns:
h : NumPy array of shape (num_nodes, 3)
columns:
[ next_if_0 , next_if_1 , symbol_index ]
"""
h = np.zeros((1, 3), dtype=int)
for code_index, code in enumerate(code_array, start=1):
word = [int(bit) for bit in code]
h_index = 0
for bit in word:
k = bit
next_node = h[h_index, k]
if next_node == 0:
h = np.vstack([h, [0, 0, 0]])
new_index = h.shape[0] - 1
h[h_index, k] = new_index
h_index = new_index
else:
h_index = next_node
h[h_index, 2] = code_index
return h
# ------------------ ENCODE ------------------
def encode_huff(coeff_sec, huff_LUT_list, force_codebook = None):
"""
Huffman-encode a sequence of quantized coefficients.
This function selects the appropriate Huffman codebook based on the
maximum absolute value of the input coefficients, encodes the coefficients
into a binary Huffman bitstream, and returns both the bitstream and the
selected codebook index.
This is the Python equivalent of the MATLAB `encodeHuff.m` function used
in audio/image coding (e.g., scale factor band encoding). The input
coefficient sequence is grouped into fixed-size tuples as defined by
the chosen Huffman LUT. Zero-padding may be applied internally.
Parameters
----------
coeff_sec : array_like of int
1-D array of quantized integer coefficients to encode.
Typically corresponds to a "section" or scale-factor band.
huff_LUT_list : list
List of Huffman lookup-table dictionaries as returned by `loadLUT()`.
Index 1..11 correspond to valid Huffman codebooks.
Index 0 is unused.
Returns
-------
huffSec : str
Huffman-encoded bitstream represented as a string of '0' and '1'
characters.
huffCodebook : int
Index (1..11) of the Huffman codebook used for encoding.
A value of 0 indicates a special all-zero section.
"""
if force_codebook is not None:
return huff_LUT_code_1(huff_LUT_list[force_codebook], coeff_sec)
maxAbsVal = np.max(np.abs(coeff_sec))
if maxAbsVal == 0:
huffCodebook = 0
huffSec = huff_LUT_code_0()
elif maxAbsVal == 1:
candidates = [1, 2]
huffSec1 = huff_LUT_code_1(huff_LUT_list[candidates[0]], coeff_sec)
huffSec2 = huff_LUT_code_1(huff_LUT_list[candidates[1]], coeff_sec)
if len(huffSec1) <= len(huffSec2):
huffSec = huffSec1
huffCodebook = candidates[0]
else:
huffSec = huffSec2
huffCodebook = candidates[1]
elif maxAbsVal == 2:
candidates = [3, 4]
huffSec1 = huff_LUT_code_1(huff_LUT_list[candidates[0]], coeff_sec)
huffSec2 = huff_LUT_code_1(huff_LUT_list[candidates[1]], coeff_sec)
if len(huffSec1) <= len(huffSec2):
huffSec = huffSec1
huffCodebook = candidates[0]
else:
huffSec = huffSec2
huffCodebook = candidates[1]
elif maxAbsVal in (3, 4):
candidates = [5, 6]
huffSec1 = huff_LUT_code_1(huff_LUT_list[candidates[0]], coeff_sec)
huffSec2 = huff_LUT_code_1(huff_LUT_list[candidates[1]], coeff_sec)
if len(huffSec1) <= len(huffSec2):
huffSec = huffSec1
huffCodebook = candidates[0]
else:
huffSec = huffSec2
huffCodebook = candidates[1]
elif maxAbsVal in (5, 6, 7):
candidates = [7, 8]
huffSec1 = huff_LUT_code_1(huff_LUT_list[candidates[0]], coeff_sec)
huffSec2 = huff_LUT_code_1(huff_LUT_list[candidates[1]], coeff_sec)
if len(huffSec1) <= len(huffSec2):
huffSec = huffSec1
huffCodebook = candidates[0]
else:
huffSec = huffSec2
huffCodebook = candidates[1]
elif maxAbsVal in (8, 9, 10, 11, 12):
candidates = [9, 10]
huffSec1 = huff_LUT_code_1(huff_LUT_list[candidates[0]], coeff_sec)
huffSec2 = huff_LUT_code_1(huff_LUT_list[candidates[1]], coeff_sec)
if len(huffSec1) <= len(huffSec2):
huffSec = huffSec1
huffCodebook = candidates[0]
else:
huffSec = huffSec2
huffCodebook = candidates[1]
elif maxAbsVal in (13, 14, 15):
huffCodebook = 11
huffSec = huff_LUT_code_1(huff_LUT_list[huffCodebook], coeff_sec)
else:
huffCodebook = 11
huffSec = huff_LUT_code_ESC(huff_LUT_list[huffCodebook], coeff_sec)
return huffSec, huffCodebook
def huff_LUT_code_1(huff_LUT, coeff_sec):
LUT = huff_LUT['LUT']
nTupleSize = huff_LUT['nTupleSize']
maxAbsCodeVal = huff_LUT['maxAbsCodeVal']
signedValues = huff_LUT['signedValues']
numTuples = int(np.ceil(len(coeff_sec) / nTupleSize))
if signedValues:
coeff = coeff_sec + maxAbsCodeVal
base = 2 * maxAbsCodeVal + 1
else:
coeff = coeff_sec
base = maxAbsCodeVal + 1
coeffPad = np.zeros(numTuples * nTupleSize, dtype=int)
coeffPad[:len(coeff)] = coeff
huffSec = []
powers = base ** np.arange(nTupleSize - 1, -1, -1)
for i in range(numTuples):
nTuple = coeffPad[i*nTupleSize:(i+1)*nTupleSize]
huffIndex = int(np.abs(nTuple) @ powers)
hexVal = LUT[huffIndex, 2]
huffLen = LUT[huffIndex, 1]
bits = format(int(hexVal), f'0{int(huffLen)}b')
if signedValues:
huffSec.append(bits)
else:
signBits = ''.join('1' if v < 0 else '0' for v in nTuple)
huffSec.append(bits + signBits)
return ''.join(huffSec)
def huff_LUT_code_0():
return ''
def huff_LUT_code_ESC(huff_LUT, coeff_sec):
LUT = huff_LUT['LUT']
nTupleSize = huff_LUT['nTupleSize']
maxAbsCodeVal = huff_LUT['maxAbsCodeVal']
numTuples = int(np.ceil(len(coeff_sec) / nTupleSize))
base = maxAbsCodeVal + 1
coeffPad = np.zeros(numTuples * nTupleSize, dtype=int)
coeffPad[:len(coeff_sec)] = coeff_sec
huffSec = []
powers = base ** np.arange(nTupleSize - 1, -1, -1)
for i in range(numTuples):
nTuple = coeffPad[i*nTupleSize:(i+1)*nTupleSize]
lnTuple = nTuple.astype(float)
lnTuple[lnTuple == 0] = np.finfo(float).eps
N4 = np.maximum(0, np.floor(np.log2(np.abs(lnTuple))).astype(int))
N = np.maximum(0, N4 - 4)
esc = np.abs(nTuple) > 15
nTupleESC = nTuple.copy()
nTupleESC[esc] = np.sign(nTupleESC[esc]) * 16
huffIndex = int(np.abs(nTupleESC) @ powers)
hexVal = LUT[huffIndex, 2]
huffLen = LUT[huffIndex, 1]
bits = format(int(hexVal), f'0{int(huffLen)}b')
escSeq = ''
for k in range(nTupleSize):
if esc[k]:
escSeq += '1' * N[k]
escSeq += '0'
escSeq += format(abs(nTuple[k]) - (1 << N4[k]), f'0{N4[k]}b')
signBits = ''.join('1' if v < 0 else '0' for v in nTuple)
huffSec.append(bits + signBits + escSeq)
return ''.join(huffSec)
# ------------------ DECODE ------------------
def decode_huff(huff_sec, huff_LUT):
"""
Decode a Huffman-encoded stream.
Parameters
----------
huff_sec : array-like of int or str
Huffman encoded stream as a sequence of 0 and 1 (string or list/array).
huff_LUT : dict
Huffman lookup table with keys:
- 'invTable': inverse table (numpy array)
- 'codebook': codebook number
- 'nTupleSize': tuple size
- 'maxAbsCodeVal': maximum absolute code value
- 'signedValues': True/False
Returns
-------
decCoeffs : list of int
Decoded quantized coefficients.
"""
h = huff_LUT['invTable']
huffCodebook = huff_LUT['codebook']
nTupleSize = huff_LUT['nTupleSize']
maxAbsCodeVal = huff_LUT['maxAbsCodeVal']
signedValues = huff_LUT['signedValues']
# Convert string to array of ints
if isinstance(huff_sec, str):
huff_sec = np.array([int(b) for b in huff_sec])
eos = False
decCoeffs = []
streamIndex = 0
while not eos:
wordbit = 0
r = 0 # start at root
# Decode Huffman word using inverse table
while True:
b = huff_sec[streamIndex + wordbit]
wordbit += 1
rOld = r
r = h[rOld, b]
if h[r, 0] == 0 and h[r, 1] == 0:
symbolIndex = h[r, 2] - 1 # zero-based
streamIndex += wordbit
break
# Decode n-tuple magnitudes
if signedValues:
base = 2 * maxAbsCodeVal + 1
nTupleDec = []
tmp = symbolIndex
for p in reversed(range(nTupleSize)):
val = tmp // (base ** p)
nTupleDec.append(val - maxAbsCodeVal)
tmp = tmp % (base ** p)
nTupleDec = np.array(nTupleDec)
else:
base = maxAbsCodeVal + 1
nTupleDec = []
tmp = symbolIndex
for p in reversed(range(nTupleSize)):
val = tmp // (base ** p)
nTupleDec.append(val)
tmp = tmp % (base ** p)
nTupleDec = np.array(nTupleDec)
# Apply sign bits
nTupleSignBits = huff_sec[streamIndex:streamIndex + nTupleSize]
nTupleSign = -(np.sign(nTupleSignBits - 0.5))
streamIndex += nTupleSize
nTupleDec = nTupleDec * nTupleSign
# Handle escape sequences
escIndex = np.where(np.abs(nTupleDec) == 16)[0]
if huffCodebook == 11 and escIndex.size > 0:
for idx in escIndex:
N = 0
b = huff_sec[streamIndex]
while b:
N += 1
b = huff_sec[streamIndex + N]
# Skip the N leading '1' bits AND the terminating '0' delimiter.
# The encoder writes: '1'*N + '0' + <N4 bits>
streamIndex += N +1
N4 = N + 4
escape_word = huff_sec[streamIndex:streamIndex + N4]
escape_value = 2 ** N4 + int("".join(map(str, escape_word)), 2)
nTupleDec[idx] = escape_value
# We already consumed the delimiter above; now consume only N4 bits.
streamIndex += N4
# Apply signs again
nTupleDec[escIndex] *= nTupleSign[escIndex]
decCoeffs.extend(nTupleDec.tolist())
if streamIndex >= len(huff_sec):
eos = True
return decCoeffs