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

📅 2026/7/19 17:20:48
数据压缩原理与应用 | 实验五 | 基于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_strategystatic, bits_per_sample8): 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, dtypeint) 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, type2, normortho)[: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, type2, normortho) 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, dtypenp.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, dtypenp.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, dtypenp.uint8), quantized_indices: np.array(quantized, dtypeobject), 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(bdct1) 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, dtypenp.uint8).tobytes()) else: f.write(np.array(q, dtypenp.uint16).tobytes()) def load_compressed(self, filename): 加载多声道压缩数据 with open(filename, rb) as f: if f.read(4) ! bdct1: 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), dtypenp.uint8) # 读取量化索引 quantized [] for sb in range(32): bits self.bit_allocation[sb] if bits 8: q np.frombuffer(f.read(num_blocks), dtypenp.uint8) else: q np.frombuffer(f.read(num_blocks * 2), dtypenp.uint16) quantized.append(q) frames_data.append({ scale_factors: scale_factors, quantized_indices: np.array(quantized, dtypeobject) }) 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, bins50, densityTrue, alpha0.7) plt.xlabel(量化索引) plt.ylabel(概率密度) plt.title(量化索引概率分布多声道合并) plt.grid(True, alpha0.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, alpha0.3) # 重建信号 plt.subplot(3, 1, 2) plt.plot(t[:show_len], y[:show_len]) plt.title(重建信号第一声道) plt.ylabel(幅度) plt.grid(True, alpha0.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, alpha0.3) plt.tight_layout() plt.show() def main(): codec MPEG1AudioLayer1PQMF(bit_allocation_strategystatic) 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(fSegSNR {segsnr:.2f} dB) print(fCompression Ratio {ratio:.2f}) print(fCompressed Size {os.path.getsize(compressed_file) / 1024:.2f} KB) # 显示声道信息 if hasattr(original, shape) and original.ndim 1: print(fOriginal channels: {original.shape[1]}) if hasattr(reconstructed, shape) and reconstructed.ndim 1: print(fReconstructed channels: {reconstructed.shape[1]}) # 绘制直方图 plot_histogram(all_indices) # 绘制重建对比图 plot_reconstruction(original, reconstructed, sr) if __name__ __main__: main()