atmos/os_lib/
aac.rs

1// aac.rs - AAC-LC デコーダ（ISO/IEC 14496-3 の GA 4.x、純粋モジュール）。
2//
3// 【2026-10-17】YouTube の音声（itag 140 = mp4a.40.2、AAC-LC）を再生するため新設。
4// ハードウェアにもグローバル状態にも依存しないので、ホスト試験で ffmpeg の復号結果と照合する。
5//
6// 対応: AAC-LC（object type 2）の raw_data_block（SCE / CPE / LFE / DSE / FIL）、長窓・短窓
7// （窓のグループ化）、サイン窓・KBD 窓、M/S ステレオ、インテンシティステレオ、PNS、TNS、
8// パルス、AudioSpecificConfig と ADTS ヘッダ。
9// 未対応（明示的にエラー）: AAC Main の予測、LTP、SSR のゲイン制御、CCE、PCE、
10// 960 サンプルのフレーム。HE-AAC（SBR/PS）は FIL の拡張として読み捨てる（LC 部分だけ出る）。
11
12use super::aac_tables::*;
13use alloc::vec;
14use alloc::vec::Vec;
15use core::f32::consts::PI;
16
17// ===================== ビット読み =====================
18
19struct Bits<'a> {
20    data: &'a [u8],
21    pos: usize,
22}
23
24impl<'a> Bits<'a> {
25    fn new(data: &'a [u8]) -> Self {
26        Self { data, pos: 0 }
27    }
28
29    fn bit(&mut self) -> Result<u32, &'static str> {
30        let b = *self.data.get(self.pos / 8).ok_or("AAC: unexpected end of data")?;
31        let v = (b >> (7 - (self.pos % 8))) & 1;
32        self.pos += 1;
33        Ok(v as u32)
34    }
35
36    fn read(&mut self, n: u32) -> Result<u32, &'static str> {
37        let mut v = 0u32;
38        for _ in 0..n.min(32) {
39            v = (v << 1) | self.bit()?;
40        }
41        Ok(v)
42    }
43
44    fn align(&mut self) {
45        self.pos = (self.pos + 7) & !7;
46    }
47
48    fn skip(&mut self, bits: usize) -> Result<(), &'static str> {
49        self.pos += bits;
50        if self.pos > self.data.len() * 8 {
51            return Err("AAC: unexpected end of data");
52        }
53        Ok(())
54    }
55}
56
57// ===================== ハフマン =====================
58
59/// 二分木のハフマン復号器。`nodes[i] = [0 の子, 1 の子]`、負の値は葉（-(記号 + 1)）、0 は未使用。
60struct Huff {
61    nodes: Vec<[i32; 2]>,
62}
63
64impl Huff {
65    fn new(codes: &[u32], bits: &[u8]) -> Self {
66        let mut nodes = vec![[0i32; 2]];
67        for (sym, (&code, &len)) in codes.iter().zip(bits.iter()).enumerate() {
68            let mut cur = 0usize;
69            for i in (0..len as u32).rev() {
70                let b = ((code >> i) & 1) as usize;
71                if i == 0 {
72                    nodes[cur][b] = -(sym as i32) - 1;
73                } else {
74                    if nodes[cur][b] <= 0 {
75                        nodes.push([0, 0]);
76                        nodes[cur][b] = (nodes.len() - 1) as i32;
77                    }
78                    cur = nodes[cur][b] as usize;
79                }
80            }
81        }
82        Self { nodes }
83    }
84
85    fn decode(&self, r: &mut Bits) -> Result<usize, &'static str> {
86        let mut cur = 0usize;
87        for _ in 0..32 {
88            let b = r.bit()? as usize;
89            let n = self.nodes.get(cur).ok_or("AAC: bad huffman node")?[b];
90            if n < 0 {
91                return Ok((-n - 1) as usize);
92            }
93            if n == 0 {
94                return Err("AAC: invalid huffman code");
95            }
96            cur = n as usize;
97        }
98        Err("AAC: huffman code too long")
99    }
100}
101
102// ===================== 設定 =====================
103
104/// 標本化周波数インデックス → Hz（表 1.18）。
105pub const SAMPLE_RATES: [u32; 13] = [96000, 88200, 64000, 48000, 44100, 32000, 24000, 22050, 16000, 12000, 11025, 8000, 7350];
106
107/// 復号に必要な設定（AudioSpecificConfig / ADTS ヘッダから）。
108#[derive(Clone, Copy, Debug, PartialEq, Eq)]
109pub struct AacConfig {
110    pub object_type: u32,
111    pub sf_index: usize,
112    pub sample_rate: u32,
113    pub channels: u32,
114}
115
116/// AudioSpecificConfig（1.6.2.1、MP4 の esds の DecoderSpecificInfo）。
117pub fn parse_audio_specific_config(data: &[u8]) -> Result<AacConfig, &'static str> {
118    let mut r = Bits::new(data);
119    let mut object_type = r.read(5)?;
120    if object_type == 31 {
121        object_type = 32 + r.read(6)?;
122    }
123    let sf_index = r.read(4)? as usize;
124    let sample_rate = if sf_index == 15 {
125        r.read(24)?
126    } else {
127        *SAMPLE_RATES.get(sf_index).ok_or("AAC: bad sampling frequency index")?
128    };
129    let channels = r.read(4)?;
130    if object_type == 2 || object_type == 4 {
131        // GASpecificConfig: frameLengthFlag（960 は未対応）。
132        if r.read(1)? == 1 {
133            return Err("AAC: 960-sample frames not supported");
134        }
135    }
136    let sf_index = if sf_index == 15 {
137        SAMPLE_RATES.iter().position(|&f| f <= sample_rate).unwrap_or(12)
138    } else {
139        sf_index
140    };
141    Ok(AacConfig { object_type, sf_index, sample_rate, channels })
142}
143
144/// ADTS ヘッダ（1.A.2.2）。戻り値は (設定, ヘッダ長, フレーム長（ヘッダ込み）)。
145pub fn parse_adts_header(data: &[u8]) -> Option<(AacConfig, usize, usize)> {
146    if data.len() < 7 || data[0] != 0xff || data[1] & 0xf6 != 0xf0 {
147        return None;
148    }
149    let protection_absent = data[1] & 1;
150    let profile = (data[2] >> 6) as u32;
151    let sf_index = ((data[2] >> 2) & 0xf) as usize;
152    let channels = (((data[2] & 1) << 2) | (data[3] >> 6)) as u32;
153    let frame_len = (((data[3] & 3) as usize) << 11) | ((data[4] as usize) << 3) | ((data[5] >> 5) as usize);
154    let hdr = if protection_absent == 1 { 7 } else { 9 };
155    let sample_rate = *SAMPLE_RATES.get(sf_index)?;
156    if frame_len < hdr {
157        return None;
158    }
159    Some((AacConfig { object_type: profile + 1, sf_index, sample_rate, channels }, hdr, frame_len))
160}
161
162// ===================== 個々のチャンネルの構文（ICS） =====================
163
164const ONLY_LONG: u8 = 0;
165const LONG_START: u8 = 1;
166const EIGHT_SHORT: u8 = 2;
167const LONG_STOP: u8 = 3;
168
169const ZERO_HCB: u8 = 0;
170const ESC_HCB: u8 = 11;
171const NOISE_HCB: u8 = 13;
172const INTENSITY_HCB2: u8 = 14;
173const INTENSITY_HCB: u8 = 15;
174
175/// ics_info（4.4.2.1）。
176#[derive(Clone, Copy)]
177struct IcsInfo {
178    window_sequence: u8,
179    window_shape: u8,
180    max_sfb: usize,
181    num_window_groups: usize,
182    group_len: [usize; 8],
183    num_windows: usize,
184    num_swb: usize,
185    swb_offset: &'static [u16],
186    tns_max_bands: usize,
187}
188
189impl IcsInfo {
190    fn parse(r: &mut Bits, sf_index: usize) -> Result<Self, &'static str> {
191        let _reserved = r.bit()?;
192        let window_sequence = r.read(2)? as u8;
193        let window_shape = r.read(1)? as u8;
194        let mut info = IcsInfo {
195            window_sequence,
196            window_shape,
197            max_sfb: 0,
198            num_window_groups: 1,
199            group_len: [1, 0, 0, 0, 0, 0, 0, 0],
200            num_windows: 1,
201            num_swb: NUM_SWB_1024[sf_index] as usize,
202            swb_offset: SWB_OFFSET_1024[sf_index],
203            tns_max_bands: TNS_MAX_BANDS_1024[sf_index] as usize,
204        };
205        if window_sequence == EIGHT_SHORT {
206            info.max_sfb = r.read(4)? as usize;
207            for _ in 0..7 {
208                if r.bit()? == 1 {
209                    info.group_len[info.num_window_groups - 1] += 1;
210                } else {
211                    info.num_window_groups += 1;
212                    info.group_len[info.num_window_groups - 1] = 1;
213                }
214            }
215            info.num_windows = 8;
216            info.num_swb = NUM_SWB_128[sf_index] as usize;
217            info.swb_offset = SWB_OFFSET_128[sf_index];
218            info.tns_max_bands = TNS_MAX_BANDS_128[sf_index] as usize;
219        } else {
220            info.max_sfb = r.read(6)? as usize;
221            if r.bit()? == 1 {
222                return Err("AAC: prediction not allowed in AAC-LC");
223            }
224        }
225        if info.max_sfb > info.num_swb {
226            return Err("AAC: max_sfb exceeds num_swb");
227        }
228        Ok(info)
229    }
230}
231
232/// TNS のフィルタ（窓ごとに最大 3 本、4.6.9）。
233#[derive(Clone, Copy, Default)]
234struct TnsFilter {
235    length: usize,
236    order: usize,
237    direction: bool,
238    coef: [f32; 12],
239}
240
241#[derive(Clone, Copy, Default)]
242struct Tns {
243    present: bool,
244    n_filt: [usize; 8],
245    filt: [[TnsFilter; 3]; 8],
246}
247
248/// 1 チャンネル分の復号結果（スペクトルは逆量子化・スケール済み）。
249struct ChannelData {
250    ics: IcsInfo,
251    coef: [f32; 1024],
252    /// [g * max_sfb + sfb]
253    band_type: [u8; 8 * 64],
254    /// スケールファクタ（通常帯域）・インテンシティの位置・ノイズのエネルギー。
255    sf: [i32; 8 * 64],
256    tns: Tns,
257}
258
259/// 1 チャンネルの出力側の状態（重ね合わせの後半と、前のフレームの窓の形）。
260#[derive(Clone)]
261struct ChannelState {
262    overlap: Vec<f32>,
263    prev_shape: u8,
264}
265
266// ===================== デコーダ =====================
267
268/// AAC-LC デコーダ。1 フレーム（raw_data_block）ごとに 1024 サンプル × チャンネル数を返す。
269pub struct AacDecoder {
270    cfg: AacConfig,
271    huff_sf: Huff,
272    huff_spec: Vec<Huff>,
273    pow43: Vec<f32>,
274    /// 窓の前半（立ち上がり）: [サイン, KBD]。長窓 1024・短窓 128。
275    win_long: [Vec<f32>; 2],
276    win_short: [Vec<f32>; 2],
277    imdct_long: Imdct,
278    imdct_short: Imdct,
279    ch: Vec<ChannelState>,
280    rand: u32,
281}
282
283impl AacDecoder {
284    pub fn new(cfg: AacConfig) -> Result<Self, &'static str> {
285        if cfg.object_type != 2 {
286            return Err("AAC: only AAC-LC (object type 2) is supported");
287        }
288        if cfg.sf_index > 12 {
289            return Err("AAC: bad sampling frequency index");
290        }
291        let sf_codes: Vec<u32> = SF_CODES.to_vec();
292        let huff_sf = Huff::new(&sf_codes, &SF_BITS);
293        let spec: [(&[u16], &[u8]); 11] = [
294            (&CODES1, &BITS1),
295            (&CODES2, &BITS2),
296            (&CODES3, &BITS3),
297            (&CODES4, &BITS4),
298            (&CODES5, &BITS5),
299            (&CODES6, &BITS6),
300            (&CODES7, &BITS7),
301            (&CODES8, &BITS8),
302            (&CODES9, &BITS9),
303            (&CODES10, &BITS10),
304            (&CODES11, &BITS11),
305        ];
306        let huff_spec = spec
307            .iter()
308            .map(|(c, b)| {
309                let codes: Vec<u32> = c.iter().map(|&v| v as u32).collect();
310                Huff::new(&codes, b)
311            })
312            .collect();
313        let pow43 = (0..8192).map(|i| libm::powf(i as f32, 4.0 / 3.0)).collect();
314        let sine = |n: usize| -> Vec<f32> { (0..n / 2).map(|i| libm::sinf(PI / n as f32 * (i as f32 + 0.5))).collect() };
315        let nch = cfg.channels.clamp(1, 2) as usize;
316        Ok(Self {
317            cfg,
318            huff_sf,
319            huff_spec,
320            pow43,
321            win_long: [sine(2048), kbd_window(2048, 4.0)],
322            win_short: [sine(256), kbd_window(256, 6.0)],
323            imdct_long: Imdct::new(1024),
324            imdct_short: Imdct::new(128),
325            ch: vec![ChannelState { overlap: vec![0.0; 1024], prev_shape: 0 }; nch],
326            rand: 0x1f2e3d4c,
327        })
328    }
329
330    pub fn config(&self) -> AacConfig {
331        self.cfg
332    }
333
334    /// raw_data_block を 1 つ復号し、-1.0..1.0 の float をチャンネル交互（1024 × チャンネル数）で返す。
335    pub fn decode_frame_f32(&mut self, data: &[u8]) -> Result<Vec<f32>, &'static str> {
336        let mut r = Bits::new(data);
337        let nch = self.ch.len();
338        let mut outs: Vec<[f32; 1024]> = Vec::new();
339        loop {
340            let id = r.read(3)?;
341            match id {
342                0 | 3 => {
343                    // SCE / LFE
344                    let _tag = r.read(4)?;
345                    let mut c = self.decode_ics(&mut r, None)?;
346                    if c.tns.present {
347                        apply_tns(&mut c);
348                    }
349                    let idx = outs.len();
350                    if idx < nch {
351                        let pcm = self.synthesize(idx, &c);
352                        outs.push(pcm);
353                    }
354                }
355                1 => {
356                    // CPE
357                    let _tag = r.read(4)?;
358                    let common = r.bit()? == 1;
359                    let mut ms_mask = [false; 8 * 64];
360                    let mut ms_present = 0;
361                    let mut info = None;
362                    if common {
363                        let i = IcsInfo::parse(&mut r, self.cfg.sf_index)?;
364                        ms_present = r.read(2)?;
365                        if ms_present == 3 {
366                            return Err("AAC: reserved ms_mask_present");
367                        }
368                        let n = i.num_window_groups * i.max_sfb;
369                        for m in ms_mask.iter_mut().take(n) {
370                            *m = ms_present == 2 || (ms_present == 1 && r.bit()? == 1);
371                        }
372                        info = Some(i);
373                    }
374                    let mut c0 = self.decode_ics(&mut r, info)?;
375                    let mut c1 = self.decode_ics(&mut r, info)?;
376                    if common {
377                        stereo_process(&mut c0, &mut c1, &ms_mask, ms_present != 0);
378                    }
379                    if c0.tns.present {
380                        apply_tns(&mut c0);
381                    }
382                    if c1.tns.present {
383                        apply_tns(&mut c1);
384                    }
385                    let idx = outs.len();
386                    if idx < nch {
387                        let pcm = self.synthesize(idx, &c0);
388                        outs.push(pcm);
389                    }
390                    let idx = outs.len();
391                    if idx < nch {
392                        let pcm = self.synthesize(idx, &c1);
393                        outs.push(pcm);
394                    }
395                }
396                2 => return Err("AAC: coupling channel element not supported"),
397                4 => {
398                    // DSE
399                    let _tag = r.read(4)?;
400                    let align = r.bit()? == 1;
401                    let mut count = r.read(8)? as usize;
402                    if count == 255 {
403                        count += r.read(8)? as usize;
404                    }
405                    if align {
406                        r.align();
407                    }
408                    r.skip(count * 8)?;
409                }
410                5 => return Err("AAC: program config element not supported"),
411                6 => {
412                    // FIL（SBR などの拡張は読み捨てる）
413                    let mut count = r.read(4)? as usize;
414                    if count == 15 {
415                        count += r.read(8)? as usize;
416                        count -= 1;
417                    }
418                    r.skip(count * 8)?;
419                }
420                _ => break,
421            }
422        }
423        if outs.is_empty() {
424            return Err("AAC: no audio element in frame");
425        }
426        // 足りないチャンネルは複製（モノラル → ステレオ指定など）。
427        while outs.len() < nch {
428            let last = outs[outs.len() - 1];
429            outs.push(last);
430        }
431        let mut out = Vec::with_capacity(1024 * nch);
432        for i in 0..1024 {
433            for o in &outs {
434                out.push(o[i] * (1.0 / 32768.0));
435            }
436        }
437        Ok(out)
438    }
439
440    /// `decode_frame_f32` の 16 ビット PCM 版。
441    pub fn decode_frame(&mut self, data: &[u8]) -> Result<Vec<i16>, &'static str> {
442        Ok(self
443            .decode_frame_f32(data)?
444            .iter()
445            .map(|&v| (v * 32768.0).clamp(-32768.0, 32767.0) as i16)
446            .collect())
447    }
448
449    fn next_rand(&mut self) -> f32 {
450        // 線形合同法（PNS の雑音。規格は乱数の種類を定めない）。
451        self.rand = self.rand.wrapping_mul(1664525).wrapping_add(1013904223);
452        (self.rand as i32) as f32
453    }
454
455    /// individual_channel_stream（4.4.2.7）。`common` は CPE の共通 ics_info。
456    fn decode_ics(&mut self, r: &mut Bits, common: Option<IcsInfo>) -> Result<ChannelData, &'static str> {
457        let global_gain = r.read(8)? as i32;
458        let ics = match common {
459            Some(i) => i,
460            None => IcsInfo::parse(r, self.cfg.sf_index)?,
461        };
462        let mut c = ChannelData {
463            ics,
464            coef: [0.0; 1024],
465            band_type: [0; 8 * 64],
466            sf: [0; 8 * 64],
467            tns: Tns::default(),
468        };
469        // section_data
470        let sect_bits = if ics.window_sequence == EIGHT_SHORT { 3 } else { 5 };
471        let esc = (1u32 << sect_bits) - 1;
472        for g in 0..ics.num_window_groups {
473            let mut k = 0;
474            while k < ics.max_sfb {
475                let cb = r.read(4)? as u8;
476                if cb == 12 {
477                    return Err("AAC: reserved codebook");
478                }
479                let mut end = k;
480                loop {
481                    let inc = r.read(sect_bits)?;
482                    end += inc as usize;
483                    if inc != esc {
484                        break;
485                    }
486                }
487                if end > ics.max_sfb {
488                    return Err("AAC: section exceeds max_sfb");
489                }
490                for sfb in k..end {
491                    c.band_type[g * ics.max_sfb + sfb] = cb;
492                }
493                k = end;
494            }
495        }
496        // scale_factor_data
497        let mut sf = global_gain;
498        let mut is_pos = 0i32;
499        let mut noise = global_gain - 90;
500        let mut noise_first = true;
501        for g in 0..ics.num_window_groups {
502            for sfb in 0..ics.max_sfb {
503                let idx = g * ics.max_sfb + sfb;
504                match c.band_type[idx] {
505                    ZERO_HCB => {}
506                    INTENSITY_HCB | INTENSITY_HCB2 => {
507                        is_pos += self.huff_sf.decode(r)? as i32 - 60;
508                        c.sf[idx] = is_pos;
509                    }
510                    NOISE_HCB => {
511                        if noise_first {
512                            noise += r.read(9)? as i32 - 256;
513                            noise_first = false;
514                        } else {
515                            noise += self.huff_sf.decode(r)? as i32 - 60;
516                        }
517                        c.sf[idx] = noise;
518                    }
519                    _ => {
520                        sf += self.huff_sf.decode(r)? as i32 - 60;
521                        if !(0..=255).contains(&sf) {
522                            return Err("AAC: scalefactor out of range");
523                        }
524                        c.sf[idx] = sf;
525                    }
526                }
527            }
528        }
529        // pulse_data
530        let mut pulses: [(usize, i32); 4] = [(0, 0); 4];
531        let mut num_pulse = 0;
532        if r.bit()? == 1 {
533            if ics.window_sequence == EIGHT_SHORT {
534                return Err("AAC: pulse data in short window");
535            }
536            num_pulse = r.read(2)? as usize + 1;
537            let start_sfb = r.read(6)? as usize;
538            if start_sfb >= ics.num_swb {
539                return Err("AAC: bad pulse_start_sfb");
540            }
541            let mut pos = ics.swb_offset[start_sfb] as usize;
542            for (i, p) in pulses.iter_mut().take(num_pulse).enumerate() {
543                let off = r.read(5)? as usize;
544                pos += off;
545                let _ = i;
546                if pos >= 1024 {
547                    return Err("AAC: bad pulse offset");
548                }
549                *p = (pos, r.read(4)? as i32);
550            }
551        }
552        // tns_data
553        if r.bit()? == 1 {
554            c.tns.present = true;
555            let short = ics.window_sequence == EIGHT_SHORT;
556            for w in 0..ics.num_windows {
557                let n_filt = r.read(if short { 1 } else { 2 })? as usize;
558                c.tns.n_filt[w] = n_filt;
559                if n_filt == 0 {
560                    continue;
561                }
562                let coef_res = r.bit()?;
563                for f in 0..n_filt {
564                    let flt = &mut c.tns.filt[w][f];
565                    flt.length = r.read(if short { 4 } else { 6 })? as usize;
566                    flt.order = r.read(if short { 3 } else { 5 })? as usize;
567                    if flt.order > if short { 7 } else { 12 } {
568                        return Err("AAC: TNS order too large");
569                    }
570                    if flt.order > 0 {
571                        flt.direction = r.bit()? == 1;
572                        let compress = r.bit()?;
573                        let len = coef_res + 3 - compress;
574                        // 4.6.9.3: 逆量子化 sin(v / iqfac)、負は iqfac_m。
575                        let half = (1u32 << (coef_res + 2)) as f32; // 2^(coef_res + 3 - 1)
576                        let iqfac = (half - 0.5) / (PI / 2.0);
577                        let iqfac_m = (half + 0.5) / (PI / 2.0);
578                        for i in 0..flt.order {
579                            let raw = r.read(len)? as i32;
580                            let v = if raw & (1 << (len - 1)) != 0 { raw - (1 << len) } else { raw };
581                            flt.coef[i] = libm::sinf(v as f32 / if v >= 0 { iqfac } else { iqfac_m });
582                        }
583                    }
584                }
585            }
586        }
587        // gain_control_data（SSR 専用）
588        if r.bit()? == 1 {
589            return Err("AAC: gain control not supported");
590        }
591        // spectral_data
592        let mut q = [0i32; 1024];
593        let mut win = 0usize;
594        for g in 0..ics.num_window_groups {
595            let glen = ics.group_len[g];
596            for sfb in 0..ics.max_sfb {
597                let cb = c.band_type[g * ics.max_sfb + sfb];
598                if cb == ZERO_HCB || cb >= NOISE_HCB {
599                    continue;
600                }
601                let start = ics.swb_offset[sfb] as usize;
602                let end = ics.swb_offset[sfb + 1] as usize;
603                for w in 0..glen {
604                    let base = (win + w) * 128;
605                    let mut k = start;
606                    while k < end {
607                        let vals = self.spectral_values(r, cb)?;
608                        let dim = if cb < 5 { 4 } else { 2 };
609                        for (j, &v) in vals.iter().take(dim).enumerate() {
610                            if let Some(slot) = q.get_mut(base + k + j) {
611                                *slot = v;
612                            }
613                        }
614                        k += dim;
615                    }
616                }
617            }
618            win += glen;
619        }
620        for &(pos, amp) in pulses.iter().take(num_pulse) {
621            q[pos] = if q[pos] > 0 { q[pos] + amp } else { q[pos] - amp };
622        }
623        // 逆量子化とスケール（4.6.1.3、4.6.2.3）、PNS（4.6.13）。
624        let mut win = 0usize;
625        for g in 0..ics.num_window_groups {
626            let glen = ics.group_len[g];
627            for sfb in 0..ics.max_sfb {
628                let idx = g * ics.max_sfb + sfb;
629                let cb = c.band_type[idx];
630                let start = ics.swb_offset[sfb] as usize;
631                let end = ics.swb_offset[sfb + 1] as usize;
632                for w in 0..glen {
633                    let base = (win + w) * 128;
634                    if cb == NOISE_HCB {
635                        let mut energy = 0.0f32;
636                        for k in start..end {
637                            let v = self.next_rand();
638                            c.coef[base + k] = v;
639                            energy += v * v;
640                        }
641                        let scale = libm::powf(2.0, 0.25 * c.sf[idx] as f32) / libm::sqrtf(energy.max(1e-30));
642                        for k in start..end {
643                            c.coef[base + k] *= scale;
644                        }
645                    } else if cb != ZERO_HCB && cb < NOISE_HCB {
646                        let gain = libm::powf(2.0, 0.25 * (c.sf[idx] - 100) as f32);
647                        for k in start..end {
648                            let v = q[base + k];
649                            let a = v.unsigned_abs() as usize;
650                            let m = if a < 8192 { self.pow43[a] } else { libm::powf(a as f32, 4.0 / 3.0) };
651                            c.coef[base + k] = if v < 0 { -m * gain } else { m * gain };
652                        }
653                    }
654                }
655            }
656            win += glen;
657        }
658        Ok(c)
659    }
660
661    /// 1 つの符号語のスペクトル値（コードブック 1〜11、4.6.3）。
662    fn spectral_values(&self, r: &mut Bits, cb: u8) -> Result<[i32; 4], &'static str> {
663        let h = self.huff_spec.get(cb as usize - 1).ok_or("AAC: bad codebook")?;
664        let idx = h.decode(r)? as i32;
665        let mut v = [0i32; 4];
666        let unsigned = !matches!(cb, 1 | 2 | 5 | 6);
667        let n = match cb {
668            1 | 2 => {
669                v = [idx / 27 - 1, (idx / 9) % 3 - 1, (idx / 3) % 3 - 1, idx % 3 - 1];
670                4
671            }
672            3 | 4 => {
673                v = [idx / 27, (idx / 9) % 3, (idx / 3) % 3, idx % 3];
674                4
675            }
676            5 | 6 => {
677                v[0] = idx / 9 - 4;
678                v[1] = idx % 9 - 4;
679                2
680            }
681            7 | 8 => {
682                v[0] = idx / 8;
683                v[1] = idx % 8;
684                2
685            }
686            9 | 10 => {
687                v[0] = idx / 13;
688                v[1] = idx % 13;
689                2
690            }
691            _ => {
692                v[0] = idx / 17;
693                v[1] = idx % 17;
694                2
695            }
696        };
697        if unsigned {
698            for x in v.iter_mut().take(n) {
699                if *x != 0 && r.bit()? == 1 {
700                    *x = -*x;
701                }
702            }
703        }
704        if cb == ESC_HCB {
705            for x in v.iter_mut().take(2) {
706                if x.abs() == 16 {
707                    let mut nbits = 4;
708                    while r.bit()? == 1 {
709                        nbits += 1;
710                        if nbits > 12 {
711                            return Err("AAC: bad escape sequence");
712                        }
713                    }
714                    let w = (1i32 << nbits) + r.read(nbits)? as i32;
715                    *x = if *x < 0 { -w } else { w };
716                }
717            }
718        }
719        Ok(v)
720    }
721
722    /// 逆 MDCT・窓掛け・重ね合わせ（4.6.11）。戻り値は 1024 サンプル（±32768 程度のスケール）。
723    fn synthesize(&mut self, ch: usize, c: &ChannelData) -> [f32; 1024] {
724        let seq = c.ics.window_sequence;
725        let shape = c.ics.window_shape as usize;
726        let prev = self.ch[ch].prev_shape as usize;
727        let mut buf = [0.0f32; 2048];
728        let (wl, ws) = (&self.win_long, &self.win_short);
729        if seq == EIGHT_SHORT {
730            let mut tmp = [0.0f32; 256];
731            for w in 0..8 {
732                self.imdct_short.run(&c.coef[w * 128..w * 128 + 128], &mut tmp);
733                let rise = if w == 0 { &ws[prev] } else { &ws[shape] };
734                let off = 448 + w * 128;
735                for i in 0..128 {
736                    buf[off + i] += tmp[i] * rise[i];
737                    buf[off + 128 + i] += tmp[128 + i] * ws[shape][127 - i];
738                }
739            }
740        } else {
741            let mut tmp = [0.0f32; 2048];
742            self.imdct_long.run(&c.coef, &mut tmp);
743            for i in 0..1024 {
744                buf[i] = tmp[i]
745                    * match seq {
746                        LONG_STOP => {
747                            if i < 448 {
748                                0.0
749                            } else if i < 576 {
750                                ws[prev][i - 448]
751                            } else {
752                                1.0
753                            }
754                        }
755                        _ => wl[prev][i],
756                    };
757                buf[1024 + i] = tmp[1024 + i]
758                    * match seq {
759                        LONG_START => {
760                            if i < 448 {
761                                1.0
762                            } else if i < 576 {
763                                ws[shape][127 - (i - 448)]
764                            } else {
765                                0.0
766                            }
767                        }
768                        _ => wl[shape][1023 - i],
769                    };
770            }
771        }
772        let st = &mut self.ch[ch];
773        let mut out = [0.0f32; 1024];
774        for i in 0..1024 {
775            out[i] = buf[i] + st.overlap[i];
776            st.overlap[i] = buf[1024 + i];
777        }
778        st.prev_shape = c.ics.window_shape;
779        out
780    }
781}
782
783/// M/S とインテンシティステレオ（4.6.8）。M/S は両方が通常の帯域のときだけ。
784fn stereo_process(c0: &mut ChannelData, c1: &mut ChannelData, ms_mask: &[bool; 8 * 64], ms: bool) {
785    let ics = c1.ics;
786    let mut win = 0usize;
787    for g in 0..ics.num_window_groups {
788        let glen = ics.group_len[g];
789        for sfb in 0..ics.max_sfb {
790            let idx = g * ics.max_sfb + sfb;
791            let start = ics.swb_offset[sfb] as usize;
792            let end = ics.swb_offset[sfb + 1] as usize;
793            let bt1 = c1.band_type[idx];
794            if bt1 == INTENSITY_HCB || bt1 == INTENSITY_HCB2 {
795                let mut sign = if bt1 == INTENSITY_HCB { 1.0 } else { -1.0 };
796                if ms && ms_mask[idx] {
797                    sign = -sign;
798                }
799                let scale = sign * libm::powf(0.5, 0.25 * c1.sf[idx] as f32);
800                for w in 0..glen {
801                    let base = (win + w) * 128;
802                    for k in start..end {
803                        c1.coef[base + k] = c0.coef[base + k] * scale;
804                    }
805                }
806            } else if ms && ms_mask[idx] && c0.band_type[idx] < NOISE_HCB && bt1 < NOISE_HCB {
807                for w in 0..glen {
808                    let base = (win + w) * 128;
809                    for k in start..end {
810                        let (l, r) = (c0.coef[base + k], c1.coef[base + k]);
811                        c0.coef[base + k] = l + r;
812                        c1.coef[base + k] = l - r;
813                    }
814                }
815            }
816        }
817        win += glen;
818    }
819}
820
821/// TNS の全極フィルタ（4.6.9.3）。反射係数を LPC に変換し、帯域の係数へ周波数方向に掛ける。
822fn apply_tns(c: &mut ChannelData) {
823    let ics = c.ics;
824    let mmm = ics.tns_max_bands.min(ics.max_sfb);
825    if mmm == 0 {
826        return;
827    }
828    for w in 0..ics.num_windows {
829        let mut bottom = ics.num_swb;
830        for f in 0..c.tns.n_filt[w] {
831            let flt = c.tns.filt[w][f];
832            let top = bottom;
833            bottom = top.saturating_sub(flt.length);
834            if flt.order == 0 {
835                continue;
836            }
837            // 反射係数 → LPC（ステップアップ）。
838            let mut lpc = [0.0f32; 12];
839            for m in 0..flt.order {
840                let k = flt.coef[m];
841                let mut tmp = [0.0f32; 12];
842                for i in 0..m {
843                    tmp[i] = lpc[i] + k * lpc[m - 1 - i];
844                }
845                tmp[m] = k;
846                lpc[..=m].copy_from_slice(&tmp[..=m]);
847            }
848            let start = ics.swb_offset[bottom.min(mmm)] as usize;
849            let end = ics.swb_offset[top.min(mmm)] as usize;
850            if end <= start {
851                continue;
852            }
853            let size = end - start;
854            let base = w * 128;
855            let mut state = [0.0f32; 12];
856            for m in 0..size {
857                let pos = if flt.direction { base + end - 1 - m } else { base + start + m };
858                let mut y = c.coef[pos];
859                for i in 0..flt.order {
860                    y -= lpc[i] * state[i];
861                }
862                for i in (1..flt.order).rev() {
863                    state[i] = state[i - 1];
864                }
865                state[0] = y;
866                c.coef[pos] = y;
867            }
868        }
869    }
870}
871
872/// KBD 窓の前半（4.6.11.3.2、`alpha` は長窓 4・短窓 6）。
873fn kbd_window(n: usize, alpha: f32) -> Vec<f32> {
874    let half = n / 2;
875    let bessel_i0 = |x: f64| -> f64 {
876        let mut sum = 1.0f64;
877        let mut term = 1.0f64;
878        for k in 1..50 {
879            term *= (x / 2.0) / k as f64;
880            sum += term * term;
881        }
882        sum
883    };
884    let a = core::f64::consts::PI * alpha as f64;
885    let w: Vec<f64> = (0..=half)
886        .map(|i| {
887            let t = (i as f64 - half as f64 / 2.0) / (half as f64 / 2.0);
888            bessel_i0(a * libm::sqrt((1.0 - t * t).max(0.0)))
889        })
890        .collect();
891    let total: f64 = w.iter().sum();
892    let mut acc = 0.0;
893    (0..half)
894        .map(|i| {
895            acc += w[i];
896            libm::sqrt(acc / total) as f32
897        })
898        .collect()
899}
900
901// ===================== 逆 MDCT（FFT による） =====================
902
903/// 逆 MDCT（M 係数 → 2M サンプル、x[n] = 2/N Σ X[k] cos(2π/N (n + n0)(k + 1/2))、N = 2M）。
904/// DCT-IV を M/2 点の複素 FFT で計算し、2M サンプルへ展開する。
905struct Imdct {
906    m: usize,
907    /// 前処理の回転 exp(-iπ(k + 1/4)/M)
908    pre: Vec<(f32, f32)>,
909    /// 後処理の回転 exp(-iπk/M)
910    post: Vec<(f32, f32)>,
911    /// FFT の回転 exp(-2πi k / (M/2))
912    tw: Vec<(f32, f32)>,
913}
914
915impl Imdct {
916    fn new(m: usize) -> Self {
917        let h = m / 2;
918        let rot = |a: f64| (libm::cos(a) as f32, libm::sin(a) as f32);
919        let pi = core::f64::consts::PI;
920        Self {
921            m,
922            pre: (0..h).map(|k| rot(-pi * (k as f64 + 0.25) / m as f64)).collect(),
923            post: (0..h).map(|k| rot(-pi * k as f64 / m as f64)).collect(),
924            tw: (0..h / 2).map(|k| rot(-2.0 * pi * k as f64 / h as f64)).collect(),
925        }
926    }
927
928    fn fft(&self, re: &mut [f32], im: &mut [f32]) {
929        let n = re.len();
930        // ビット反転
931        let mut j = 0usize;
932        for i in 1..n {
933            let mut bit = n >> 1;
934            while j & bit != 0 {
935                j ^= bit;
936                bit >>= 1;
937            }
938            j |= bit;
939            if i < j {
940                re.swap(i, j);
941                im.swap(i, j);
942            }
943        }
944        let mut len = 2;
945        while len <= n {
946            let step = n / len;
947            for s in (0..n).step_by(len) {
948                for k in 0..len / 2 {
949                    let (wr, wi) = self.tw[k * step];
950                    let (a, b) = (s + k, s + k + len / 2);
951                    let tr = re[b] * wr - im[b] * wi;
952                    let ti = re[b] * wi + im[b] * wr;
953                    re[b] = re[a] - tr;
954                    im[b] = im[a] - ti;
955                    re[a] += tr;
956                    im[a] += ti;
957                }
958            }
959            len <<= 1;
960        }
961    }
962
963    /// `x`（M 係数）→ `out`（2M サンプル）。
964    fn run(&self, x: &[f32], out: &mut [f32]) {
965        let m = self.m;
966        let h = m / 2;
967        let mut re = vec![0.0f32; h];
968        let mut im = vec![0.0f32; h];
969        for k in 0..h {
970            let (a, b) = (x[2 * k], x[m - 1 - 2 * k]);
971            let (c, s) = self.pre[k];
972            re[k] = a * c - b * s;
973            im[k] = a * s + b * c;
974        }
975        self.fft(&mut re, &mut im);
976        let mut u = vec![0.0f32; m];
977        for n in 0..h {
978            let (c, s) = self.post[n];
979            let (r, i) = (re[n] * c - im[n] * s, re[n] * s + im[n] * c);
980            u[2 * n] = r;
981            u[m - 1 - 2 * n] = -i;
982        }
983        let scale = 2.0 / (2 * m) as f32;
984        for n in 0..2 * m {
985            let v = if n < m / 2 {
986                u[n + m / 2]
987            } else if n < 3 * m / 2 {
988                -u[3 * m / 2 - 1 - n]
989            } else {
990                -u[n - 3 * m / 2]
991            };
992            out[n] = v * scale;
993        }
994    }
995}
996
997/// 試験用: 直接計算の逆 MDCT（O(N^2)）。
998pub fn imdct_reference(x: &[f32]) -> Vec<f32> {
999    let m = x.len();
1000    let n = 2 * m;
1001    let n0 = (n as f64 / 2.0 + 1.0) / 2.0;
1002    (0..n)
1003        .map(|i| {
1004            let s: f64 = (0..m)
1005                .map(|k| x[k] as f64 * libm::cos(2.0 * core::f64::consts::PI / n as f64 * (i as f64 + n0) * (k as f64 + 0.5)))
1006                .sum();
1007            (s * 2.0 / n as f64) as f32
1008        })
1009        .collect()
1010}
1011
1012/// 試験用: FFT 版の逆 MDCT。
1013pub fn imdct_fast(x: &[f32]) -> Vec<f32> {
1014    let t = Imdct::new(x.len());
1015    let mut out = vec![0.0; 2 * x.len()];
1016    t.run(x, &mut out);
1017    out
1018}
atmos/os_lib/aac.rs

atmos/os_lib/
aac.rs
