数据压缩原理与应用 | 实验五 | 基于DCT量化的音频编解码

目录

一、音频编解码算法设计

二、结果分析

数据集一览

实验结果

实验分析

三、经验总结

四、碎碎念

五、实验代码

一、音频编解码算法设计

本次实验采用了与JPEG编码类似的方式，即将音频信号分段后进行DCT变换并舍弃一定数量的高频系数，并且采取了各声道分别压缩的方式。受限于标准文档价格、CPP与Python的语言差异及多相滤波及解调的复杂性，本次实验没有采用512点滤波器，而是在尽量保留流程基础之上直接进行分段DCT并保留前32个系数，使用384点分帧；对每个块进行正弦窗加窗，并保留50%块重叠，以便能够完美复原。在DCT变换后，将每个块的DCT系数作为一行，形成矩阵格式，然后以列为单位开始后续编码：采取对数缩放并归一化，然后进行静态量化。解码时逆向复原。

二、结果分析

数据集一览

图2.1 数据集一览

数据集包含两段音频，较短的一段为较长一段的前10%长度内容。

实验结果

图2.2 量化索引概率分布

图2.3 重建信号对比

实验采取静态量化列表，对低频系数予以较高bit，而对高频系数予以较低bit，实验中从低到高比特位数从8逐渐减少到3。由图3.2可以看出，当低频也给予8bit时，已经出现了量化溢出现象，256值索引的概率密度显著增大，虽然SNR只有9.32，但实际听感差别很小，可能是由于人声的主观感受差别在几千Hz内更显著。

实验分析

本次实验在分块进行64点DCT变换后，只保留了前32点，但对于听感几乎没有影响，实际影响在于正弦窗不能完美复原，以及滑动窗口只在分帧内部实现，不同帧之间没有衔接。代码中采取的文件写入方式是单个DCT系数直接写入字节，而非位写入，会造成大量的位数浪费，所以无论怎样降低高频部分的量化位数也不会改变压缩率。

实验中发现，当把量化位数限制在8以内时，压缩率为整数，在压缩QQ音乐导出的音乐文件时压缩率为2，压缩某品牌手机的WAV文件时压缩率为3，由此可见，本次实验的压缩效果是写入位数不同带来的，因为没有采用位写入，所以压缩率实际上没有参考价值。此外，DCT保留位数也会对压缩率产生影响，但代码更改繁琐，所以没有对这方面进行探究。

三、经验总结

其实这次实验本来是让实现一个多相滤波的编解码器，但是我找不到具体的公式讲解，实现不了子带合成，老师给的参考项目又都是用CPP写的成熟项目，还是两个不同的项目，编解码器不兼容，里面的512点滤波器数值还对不上，实在是看不懂。AI只能给出用DCT编解码的可用代码，在跟AI缠斗十几个小时后我放弃了挣扎。

四、碎碎念

老师真是觉得AI无所不能了，居然直接让我们实现这种成熟复杂的编解码器，太高看我们了......

五、实验代码

import numpy as np import soundfile as sf import matplotlib.pyplot as plt import struct import os from scipy.fftpack import dct, idct plt.rcParams['font.sans-serif'] = ['SimHei'] plt.rcParams['axes.unicode_minus'] = False class MPEG1AudioLayer1PQMF: def __init__(self, bit_allocation_strategy='static', bits_per_sample=8): self.bit_allocation_strategy = bit_allocation_strategy self.bits_per_sample = bits_per_sample self.block_size = 64 self.hop_size = 32 self.frame_size = 384 self.analysis_window = self.create_prototype_window() self.synthesis_window = self.analysis_window.copy() self._init_bit_allocation() def create_prototype_window(self): n = np.arange(self.block_size) window = np.sin(np.pi / self.block_size * (n + 0.5)) return window def _init_bit_allocation(self): if self.bit_allocation_strategy == 'uniform': self.bit_allocation = np.full(32, self.bits_per_sample, dtype=int) else: self.bit_allocation = np.array([ 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 7, 7, 7, 7, 6, 6, 6, 6, 5, 5, 5, 5, 4, 4, 4, 4, 3, 3, 3, 3 ]) def analysis_filter_bank(self, audio_frame): subband_samples = [] for start in range(0, len(audio_frame) - self.block_size + 1, self.hop_size): block = audio_frame[start: start + self.block_size] block = block * self.analysis_window coeffs = dct(block, type=2, norm='ortho')[:32] subband_samples.append(coeffs) return np.array(subband_samples) def synthesis_filter_bank(self, subband_samples): num_blocks = len(subband_samples) output_length = (num_blocks - 1) * self.hop_size + self.block_size reconstructed = np.zeros(output_length) for i, coeffs in enumerate(subband_samples): full_coeffs = np.zeros(self.block_size) full_coeffs[:32] = coeffs block = idct(full_coeffs, type=2, norm='ortho') block *= self.synthesis_window start = i * self.hop_size reconstructed[start: start + self.block_size] += block return reconstructed def encode_scale_factor(self, scale_factor): if scale_factor < 1e-10: return 0 log_sf = np.log2(scale_factor) sf_index = int(np.clip(round((log_sf + 32) * 2), 0, 63)) return sf_index def decode_scale_factor(self, sf_index): log_sf = sf_index / 2 - 32 return 2 ** log_sf def uniform_quantizer(self, samples, bits): if bits == 0: return np.zeros_like(samples, dtype=np.uint16) levels = 2 ** bits normalized = np.clip(samples, -1, 1) q = ((normalized + 1) * 0.5 * (levels - 1)) return np.round(q).astype(np.uint16) def inverse_uniform_quantizer(self, q, scale_factor, bits): if bits == 0: return np.zeros_like(q, dtype=np.float64) levels = 2 ** bits normalized = q / (levels - 1) normalized = normalized * 2 - 1 return normalized * scale_factor def encode_frame(self, audio_frame): subbands = self.analysis_filter_bank(audio_frame) scale_factors = [] quantized = [] all_indices = [] for sb in range(32): samples = subbands[:, sb] # 取某一列 sf = np.max(np.abs(samples)) sf = max(sf, 1e-10) sf_index = self.encode_scale_factor(sf) scale_factors.append(sf_index) decoded_sf = self.decode_scale_factor(sf_index) normalized = samples / decoded_sf bits = self.bit_allocation[sb] q = self.uniform_quantizer(normalized, bits) quantized.append(q) all_indices.extend(q.tolist()) return { 'scale_factors': np.array(scale_factors, dtype=np.uint8), 'quantized_indices': np.array(quantized, dtype=object), 'all_indices': np.array(all_indices) } def decode_frame(self, encoded): num_blocks = len(encoded['quantized_indices'][0]) # 从列len获取行数量 reconstructed_subbands = [] for block in range(num_blocks): reconstructed_subbands.append(np.zeros(32)) reconstructed_subbands = np.array(reconstructed_subbands) for sb in range(32): sf = self.decode_scale_factor(encoded['scale_factors'][sb]) q = encoded['quantized_indices'][sb] bits = self.bit_allocation[sb] reconstructed = self.inverse_uniform_quantizer(q, sf, bits) reconstructed_subbands[:, sb] = reconstructed return self.synthesis_filter_bank(reconstructed_subbands) def save_compressed(self, filename, channels_data, sample_rate, num_samples, num_channels): """保存多声道压缩数据""" with open(filename, 'wb') as f: # 写入文件头 f.write(b'dct1') f.write(struct.pack('<I', sample_rate)) f.write(struct.pack('<I', num_samples)) f.write(struct.pack('<I', num_channels)) # 声道数 f.write(struct.pack('<I', len(channels_data[0]))) # 帧数（所有声道相同） # 写入每个声道的数据 for channel_idx, frames_data in enumerate(channels_data): # 写入声道标识 f.write(struct.pack('<I', channel_idx)) for frame in frames_data: # 写入缩放因子 f.write(frame['scale_factors'].astype(np.uint8).tobytes()) # 写入量化索引 for sb in range(32): bits = self.bit_allocation[sb] q = frame['quantized_indices'][sb] if bits <= 8: f.write(np.array(q, dtype=np.uint8).tobytes()) else: f.write(np.array(q, dtype=np.uint16).tobytes()) def load_compressed(self, filename): """加载多声道压缩数据""" with open(filename, 'rb') as f: if f.read(4) != b'dct1': raise ValueError("非法文件") sample_rate = struct.unpack('<I', f.read(4))[0] num_samples = struct.unpack('<I', f.read(4))[0] num_channels = struct.unpack('<I', f.read(4))[0] num_frames = struct.unpack('<I', f.read(4))[0] # 计算每个帧的时间块数 num_blocks = 11 # (frame_size - block_size) / hop_size + 1 channels_data = [] for ch in range(num_channels): # 读取声道标识 channel_idx = struct.unpack('<I', f.read(4))[0] frames_data = [] for _ in range(num_frames): # 读取缩放因子 scale_factors = np.frombuffer(f.read(32), dtype=np.uint8) # 读取量化索引 quantized = [] for sb in range(32): bits = self.bit_allocation[sb] if bits <= 8: q = np.frombuffer(f.read(num_blocks), dtype=np.uint8) else: q = np.frombuffer(f.read(num_blocks * 2), dtype=np.uint16) quantized.append(q) frames_data.append({ 'scale_factors': scale_factors, 'quantized_indices': np.array(quantized, dtype=object) }) channels_data.append(frames_data) return sample_rate, num_samples, num_channels, channels_data def compress_to_file(self, input_wav, compressed_file): """压缩多声道音频文件""" audio, sr = sf.read(input_wav) # 确保是2D数组 [样本数, 声道数] if audio.ndim == 1: audio = audio.reshape(-1, 1) num_channels = audio.shape[1] num_samples = audio.shape[0] num_frames = int(np.ceil(num_samples / self.frame_size)) # 填充音频 padded = np.zeros((num_frames * self.frame_size, num_channels)) padded[:num_samples, :] = audio # 为每个声道创建独立的数据 channels_data = [] all_channel_indices = [] for ch in range(num_channels): print(f"压缩声道 {ch + 1}/{num_channels}...") channel_audio = padded[:, ch].astype(np.float64) frames_data = [] channel_indices = [] for i in range(num_frames): start = i * self.frame_size frame = channel_audio[start:start + self.frame_size] encoded = self.encode_frame(frame) frames_data.append(encoded) channel_indices.extend(encoded['all_indices'].tolist()) channels_data.append(frames_data) all_channel_indices.append(np.array(channel_indices)) # 保存压缩文件 self.save_compressed(compressed_file, channels_data, sr, num_samples, num_channels) return all_channel_indices def decompress_file(self, compressed_file, output_wav): """解压多声道音频文件""" sr, num_samples, num_channels, channels_data = self.load_compressed(compressed_file) print(f"解压 {num_channels} 个声道...") # 为每个声道解码 reconstructed_channels = [] for ch in range(num_channels): print(f"解压声道 {ch + 1}/{num_channels}...") frames_data = channels_data[ch] reconstructed = [] for frame in frames_data: rec = self.decode_frame(frame) reconstructed.extend(rec) reconstructed = np.array(reconstructed) reconstructed = reconstructed[:num_samples] # 截断到原始长度 # 归一化 maxv = np.max(np.abs(reconstructed)) if maxv > 1e-10: reconstructed /= maxv reconstructed *= 0.9 reconstructed_channels.append(reconstructed) # 合并多声道 if num_channels == 1: output_audio = reconstructed_channels[0] else: output_audio = np.column_stack(reconstructed_channels) # 保存为WAV文件 sf.write(output_wav, output_audio.astype(np.float32), sr) return output_audio def calculate_compression_ratio(self, original_file, compressed_file): original_size = os.path.getsize(original_file) compressed_size = os.path.getsize(compressed_file) return original_size / compressed_size def seg_snr(x, y): """计算分段信噪比""" # 确保维度匹配 if x.ndim != y.ndim: if x.ndim == 1 and y.ndim == 2: y = y[:, 0] elif y.ndim == 1 and x.ndim == 2: x = x[:, 0] if x.ndim > 1: snrs = [] for ch in range(min(x.shape[1], y.shape[1])): snr = seg_snr(x[:, ch], y[:, ch]) snrs.append(snr) return np.mean(snrs) # 单声道处理 frame_len = 256 hop = 128 snrs = [] min_len = min(len(x), len(y)) x = x[:min_len] y = y[:min_len] for i in range(0, len(x) - frame_len, hop): a = x[i:i + frame_len] b = y[i:i + frame_len] noise = a - b p_signal = np.sum(a ** 2) p_noise = np.sum(noise ** 2) if p_noise < 1e-10: continue snr = 10 * np.log10(p_signal / p_noise) snrs.append(snr) return np.mean(snrs) if snrs else float('inf') def plot_histogram(all_indices): """绘制直方图（支持多声道）""" if isinstance(all_indices, list): # 多声道情况，合并所有声道的数据 all_data = np.concatenate([indices.flatten() for indices in all_indices if len(indices) > 0]) else: all_data = all_indices.flatten() plt.figure(figsize=(8, 4)) plt.hist(all_data, bins=50, density=True, alpha=0.7) plt.xlabel("量化索引") plt.ylabel("概率密度") plt.title("量化索引概率分布（多声道合并）") plt.grid(True, alpha=0.3) plt.show() def plot_reconstruction(x, y, fs): """绘制重建信号对比""" # 确保都是2D或都是1D if x.ndim == 1 and y.ndim == 2: y = y[:, 0] elif y.ndim == 1 and x.ndim == 2: x = x[:, 0] elif x.ndim == 2 and y.ndim == 2: # 都取第一声道 x = x[:, 0] y = y[:, 0] n = min(len(x), len(y)) x = x[:n] y = y[:n] t = np.arange(n) / fs show_len = int(min(5 * fs, n)) plt.figure(figsize=(12, 8)) # 原始信号 plt.subplot(3, 1, 1) plt.plot(t[:show_len], x[:show_len]) plt.title("原始信号（第一声道）") plt.ylabel("幅度") plt.grid(True, alpha=0.3) # 重建信号 plt.subplot(3, 1, 2) plt.plot(t[:show_len], y[:show_len]) plt.title("重建信号（第一声道）") plt.ylabel("幅度") plt.grid(True, alpha=0.3) # 重建误差 error = x - y plt.subplot(3, 1, 3) plt.plot(t[:show_len], error[:show_len]) plt.title("重建误差") plt.xlabel("时间 (秒)") plt.ylabel("误差幅度") plt.grid(True, alpha=0.3) plt.tight_layout() plt.show() def main(): codec = MPEG1AudioLayer1PQMF(bit_allocation_strategy='static') input_file = "传奇爱尔兰哨笛.wav" compressed_file = "compressed.dct1" reconstructed_file = "reconstructed.wav" print("开始压缩（多声道独立编码）...") all_indices = codec.compress_to_file(input_file, compressed_file) print("压缩完成") print("开始解压...") reconstructed = codec.decompress_file(compressed_file, reconstructed_file) print("解压完成") # 读取原始文件用于评估 original, sr = sf.read(input_file) # 计算分段信噪比 segsnr = seg_snr(original, reconstructed) # 计算压缩比 ratio = codec.calculate_compression_ratio(input_file, compressed_file) print("\n======================") print(f"SegSNR = {segsnr:.2f} dB") print(f"Compression Ratio = {ratio:.2f}") print(f"Compressed Size = {os.path.getsize(compressed_file) / 1024:.2f} KB") # 显示声道信息 if hasattr(original, 'shape') and original.ndim > 1: print(f"Original channels: {original.shape[1]}") if hasattr(reconstructed, 'shape') and reconstructed.ndim > 1: print(f"Reconstructed channels: {reconstructed.shape[1]}") # 绘制直方图 plot_histogram(all_indices) # 绘制重建对比图 plot_reconstruction(original, reconstructed, sr) if __name__ == "__main__": main()