通信原理PCM（对数量化：A律13折线）的matlab仿真

A律量化

要点如下：

1、将归一化x的输入范围等比地划分为8个段。每段用折线近似。（共8条折线）
2、各段上采用16个电平的均匀量化，共8×16=128个电平。
3、正负两部分结合在一起，共256个电平，对应于=8个量化比特。

量化如下：

PCM编码（对应A律量化）

代码模块

% 读取音频文件->对数量化->插入同步序列->脉冲整形->AWGN信道->匹配滤波->同步->解压缩以恢复信号
clear
clc
close all
format long

%% 参数设置
fs = 8000;                % 采样频率
duration = 1;             % 信号持续时间(s)
t = 0:1/fs:duration-1/fs; % 时间向量
SNR = 5;                 % 信噪比(dB)

%% 信号产生
% % 读取音频文件
% [audio1_real, fs] = audioread('audio1.wav');
% [audio2_real, fs] = audioread('audio2.wav');

% 生成两个原始音频信号
audio1_real = 1.8*sin(2*pi*500*t);    % 500Hz正弦信号
audio2_real = 1.5*cos(2*pi*400*t);     % 400Hz余弦信号

% 归一化信号到[-1, 1]范围
audio1 = audio1_real / max(abs(audio1_real));
audio2 = audio2_real / max(abs(audio2_real));

%% A律压缩(1位符号位+7位数值位)
% 获取符号位和数值位
signs1 = sign(audio1) >= 0;  % 1表示正，0表示负
signs2 = sign(audio2) >= 0;
% 压缩幅值（7位）
audio1_compressed = max(0, min(127, round(128 * compress_Alaw(abs(audio1)))));
audio2_compressed = max(0, min(127, round(128 * compress_Alaw(abs(audio2)))));

% 生成8bit同步序列
sync_sequence = [1 0 1 0 1 0 1 0];  % 8bit同步码

% 构建数据帧
num_samples = length(audio1_compressed);
frame_data = zeros(num_samples, 24);  % 24bit每帧

for i = 1:num_samples
    % 组装音频1的8bit（符号位+7位数值）
    audio1_bits = [signs1(i), de2bi(audio1_compressed(i), 7, 'left-msb')];
    % 组装音频2的8bit（符号位+7位数值）
    audio2_bits = [signs2(i), de2bi(audio2_compressed(i), 7, 'left-msb')];
    % 组装帧：同步序列 + 音频1 + 音频2
    frame_data(i,:) = [sync_sequence audio1_bits audio2_bits];
end

% 串行化数据
serial_data = reshape(frame_data.', [], 1);

%% 信号发射
% 参数设置
rolloff = 0.25; % 滚降系数
span = 4; % 滤波器的跨度
sps = 9; % 每个符号的采样点数
% 上采样
data_upsampled = upsample(serial_data, sps);
% 平方根升余弦FIR滤波器
h = rcosdesign(rolloff, span, sps);
% 脉冲整形
data_filtered = conv(data_upsampled, h,'same');

%% 信道传输
% 通过AWGN信道传输
received_signal = awgn(data_filtered, SNR, 'measured');

%% 信号接收
% 匹配滤波
filtered_signal = conv(received_signal, fliplr(h),'same');
% 下采样
downsampled_signal = downsample(filtered_signal,sps);
% 二值化判决
decoded_bits = downsampled_signal > 0.5;

%% A律解压缩
% 重新整形为帧格式
received_frames = reshape(decoded_bits(1:floor(length(decoded_bits)/24)*24), 24, []).';

% 帧同步和数据提取（初始化）
recovered_audio1 = zeros(num_samples, 1);
recovered_audio2 = zeros(num_samples, 1);
recovered_signs1 = zeros(num_samples, 1);
recovered_signs2 = zeros(num_samples, 1);
sync_errors = 0;

for i = 1:size(received_frames, 1)
    % 提取当前帧的同步序列
    current_sync = received_frames(i, 1:8);
    % 上一帧是否正确同步
    idx=1;
    % 检查同步序列
    % 我们这里不妨在上一帧正确同步的情况下允许1位同步误差，否则重新同步，该帧舍弃
    if (sum(abs(current_sync - sync_sequence)) == 0)||...
            (sum(abs(current_sync - sync_sequence)) <= 1)&&idx
        % 分别读取两个序列
        audio1_bits = received_frames(i, 9:16);
        audio2_bits = received_frames(i, 17:24);
        
        % 分离符号位和数值位
        recovered_signs1(i) = audio1_bits(1);
        recovered_signs2(i) = audio2_bits(1);
        
        % 转换数值位回十进制
        recovered_audio1(i) = bi2de(audio1_bits(2:8), 'left-msb');
        recovered_audio2(i) = bi2de(audio2_bits(2:8), 'left-msb');
        idx=1;
    else
        sync_errors = sync_errors + 1;
        idx=0;
    end
end

%% A律解压缩并恢复符号
recovered_signs1 = 2 * recovered_signs1 - 1;  % 转换为±1
recovered_signs2 = 2 * recovered_signs2 - 1;
recovered_audio1 = max(abs(audio1_real))*expand_Alaw(recovered_audio1/128) .* recovered_signs1;
recovered_audio2 = max(abs(audio2_real))*expand_Alaw(recovered_audio2/128) .* recovered_signs2;

%% 计算结果
% 计算nmse
nmse_audio1 = calculate_nmse(audio1_real',recovered_audio1);
nmse_audio2 = calculate_nmse(audio2_real',recovered_audio2);
% 计算误码率
ber = sum(decoded_bits ~= serial_data) / length(serial_data);

%% 绘制结果
figure;
subplot(2,2,1);
plot(received_signal(1:200));
title('接收信号片段');

subplot(2,2,2);
plot(filtered_signal(1:200));
title('滤波后信号片段');

subplot(2,2,3);
plot(t(1:200), audio1_real(1:200), 'b', t(1:200), recovered_audio1(1:200), 'r--');
title('音频1: 原始信号 / 恢复信号');
legend('原始信号', '恢复信号');

subplot(2,2,4);
plot(t(1:200), audio2_real(1:200), 'b', t(1:200), recovered_audio2(1:200), 'r--');
title('音频2: 原始信号 / 恢复信号');
legend('原始信号', '恢复信号');

fprintf('同步错误次数: %d\n', sync_errors);
fprintf('音频1nmse: %.4f\n', nmse_audio1);
fprintf('音频2nmse: %.4f\n', nmse_audio2);
fprintf('误码率: %.4f\n', ber);

%% A律压缩函数(13折线近似)
function compressed = compress_Alaw(x)
    compressed = zeros(size(x));
    for i = 1:length(x)
        x_abs = abs(x(i));
        switch 1
            case x_abs<=1/128
                 compressed(i) = floor(x_abs/(1/2048))/128;
            case (x_abs<=1/64)&&(x_abs>1/128)
                 compressed(i) = 1/8+floor((x_abs-1/128)/(1/2048))/128;
            case (x_abs <= 1/32)&&(x_abs>1/64)
                 compressed(i) = 1/4 + floor((x_abs - 1/64) * 1024) / 128;
            case (x_abs <= 1/16)&&(x_abs>1/32)
                 compressed(i) = 3/8 + floor((x_abs - 1/32) * 512) / 128;
            case (x_abs <= 1/8)&&(x_abs>1/16)
                 compressed(i) = 1/2 + floor((x_abs - 1/16) * 256) / 128;
            case (x_abs <= 1/4)&&(x_abs>1/8)
                 compressed(i) = 5/8 + floor((x_abs - 1/8) * 128) / 128;
            case (x_abs <= 1/2)&&(x_abs>1/4)
                 compressed(i) = 3/4 + floor((x_abs - 1/4) * 64) / 128;
            case (x_abs <= 1)&&(x_abs>1/2)
                 compressed(i) = 7/8 + floor((x_abs - 1/2) * 32) / 128;
        end
    end 
end

%% A律解压缩函数(13折线近似)
function expanded = expand_Alaw(y)
    expanded = zeros(size(y));
    for i = 1:length(y)
        y_abs = abs(y(i));
        switch 1
            case y_abs <= 1/8
                expanded(i) = (y_abs * 128) / 2048;
            case (y_abs > 1/8) && (y_abs <= 1/4)
                expanded(i) = 1/128 + ((y_abs - 1/8) * 128) / 2048;
            case (y_abs > 1/4) && (y_abs <= 3/8)
                expanded(i) = 1/64 + ((y_abs - 1/4) * 128) / 1024;
            case (y_abs > 3/8) && (y_abs <= 1/2)
                expanded(i) = 1/32 + ((y_abs - 3/8) * 128) / 512;
            case (y_abs > 1/2) && (y_abs <= 5/8)
                expanded(i) = 1/16 + ((y_abs - 1/2) * 128) / 256;
            case (y_abs > 5/8) && (y_abs <= 3/4)
                expanded(i) = 1/8 + ((y_abs - 5/8) * 128) / 128;
            case (y_abs > 3/4) && (y_abs <= 7/8)
                expanded(i) = 1/4 + ((y_abs - 3/4) * 128) / 64;
            case (y_abs > 7/8) && (y_abs <= 1)
                expanded(i) = 1/2 + ((y_abs - 7/8) * 128) / 32;
        end
    end
end 

%% 计算nmse
function nmse_value = calculate_nmse(original, recovered)
    % 计算原始信号和恢复信号的差的平方
    diff_squared = (original - recovered) .^ 2;
    
    % 计算原始信号的平方
    original_squared = original .^ 2;
    
    % 计算NMSE
    nmse_value = sum(diff_squared) / sum(original_squared);
end

结果展示