数字调制2ASK误码率分析matlab实现

写在开头：

·关于香农三大定理，以个人学习情况和个人理解以学习总结的形式给出。

香农三大贡献（定理）

一 香农第一定理（可变长无失真信源编码定理）

本课程第一节课，老师给出了信息量的度量方式。信息量的度量应该与消息的种类，重要程度无关，与种类无关是比较容易接受的，如果信息量与消息的重要程度失去关联这是开始让人难以接受的，但是重要程度这一概念是掺杂有很大的主观因素的。香农跳脱开了信息重要性的束缚，转而投靠了“概率”这一概念，一般地我们认为，越小概率发生的时间包含有更多的信息，而越确定发生的事件对我们来说越是“废话”，包含的信息量更少；除此之外，我们希望信息量之间的叠加对应其概率的相乘（相互独立的条件下）。以上是我们对信息量度量的绝大部分性质的期待，恰好对数函数拥有这个性质。 信息熵对应平均信息量的概念，也即符号信息量的期望。
之所以要对“信息”这一对外行来说抽象的概念进行度量，其目的是为了communication中信息的压缩。举个课上曾举过的例子，运输一箱灯泡的过程中，保险安全起见我们会用塑料泡沫包裹灯泡，同样运输一大箱同型号的灯泡，但是因为其中泡沫的体积占比不同，导致我们的运输效率（同体积箱子中灯泡个数）不同，那么我们是否可以通过某种方式，得到灯泡的极限运输速率？即保证运输灯泡的安全性的前提下，无限压缩泡沫的体积。
香农第一定理指明了
·一段信息的信息量是固定的，这称为这段信息的信息熵（H）
·无论怎么压缩，信息熵是无失真信源编码的极限值
·若编码的平均码长小于信息熵值，必然发生差错（也就是有损）

感悟思考：从莫尔斯码到香农编码到霍夫曼编码

莫尔斯码中常用的字母多用比较简单的短线表示，使用频率低的字母长线较多，我认为这体现了一种资源配置的思想，使用频率或者讲出现概率越高的代码，我们给它分配更快速更便捷实现的方法，以期望平均的工作效率提升。香农将这一思想再次量化，将信源符号按照出现的概率进行了排列，出现频率高的采用更短的码长，同时这里暗含了一个编码压缩的极限，也就是之后的霍夫曼编码方式，霍夫曼编码充分利用短码，也即将短码全部占用完成后，再扩充码长，从而达到了比香农编码理论中的压缩极限。

个人认为关于香农第一贡献，就是天才般的将信息量这一抽象的概念与概率建立了联系并且将其量化，有了这个度量工具，进而才能指标性的优化压缩方法，对压缩数据提供了方向和评判标准。

二 香农第二定理（有噪信道编码定理）

如果说第一定理指出了压缩编码（无损）的极限，那么第二定律则指出了信号传输速率的极限。公式2-1与我们的直观认知是相符的，首先是信道容量与信道带宽呈正相关；而后与信噪比也呈正相关，当噪声占比大时，传输越来越无效。因此也指出了提高传输速率的两个途径，即增加信道带宽与提高信噪比。从了解到资料中提及了一个例子，即网速与距离路由器远近之间有什么关系？为什么距离路由器远网速越慢？解释比较容易理解，在信息传输过程中会叠加许多噪声，随着距离增加，噪声功率增加，而信号由于信道对其衰减作用功率也会下降，造成信噪比恶化，进而导致信道容量减小，传输速率上限降低，网速变慢。
这个公式给出了传输速率上限信道容量的影响因素，对于实际的意义还不仅于此。这让我想到了信息传输可靠性与有效性之间的博弈，第二定理启发性如下：

·有噪信道编码定理指出，尽管噪声会干扰通信信道，但还是有可能在信息传输速率小于信道容量的前提下，以任意低的错误概率传送数据信息。

牺牲了信息传输速率，则有效性降低，换取的是传输更加可靠。

三 香农第三定理（保失真度准则下的有失真信源编码定理）

以上定理：
·给出了信息压缩的极限：

·R(D)在实际工程中可以作为衡量各种压缩编码方法性能优劣的一种标尺
·量化、数模转换、频带压缩和数据压缩的理论基础

关于第三定理的理解不太深入，大概意思理解如下：比如玩而你画我猜游戏，这个游戏的本质即信息的传递，我们假定表达者和接受者之间的认知差异为噪声，比如传递信息“三长两短”，如果我们一下子表达清楚这个概念是相对困难的，于是我们可以和同伴事先规定一个“通信协议”或者说“校验”，我们第一个表达动作代表字符个数，然后依次传输每一个字符信息，如此，信息的传递会变得更加可靠。回到原理本身，第三定理指明了在满足第一定理压缩极限，第二定理传输速率极限的前提下，如何通信不出错。
参考：
如何理解「香农定理」，包含哪些内容，它的发现有什么意义？
https://www.zhihu.com/question/447161080/answer/1759158425
三种编码方式的比较
https://zhuanlan.zhihu.com/p/52234093
如何连贯地理解香农三大定理？
https://www.zhihu.com/question/39296849/answer/1467262363

仿真实验部分：

·1关于Huffman压缩编码

·Huffman编码思路
霍夫曼编码使用一种特别的方法为信号源中的每个符号设定二进制码。出现频率更大的符号将获得更短
的比特，出现频率更小的符号将被分配更长的比特，以此来提高数据压缩率，提高传输效率。

算法步骤：

1.对字符出现的频率降序排列，进行步骤2
2.查找频率位于最底层的两个字符（字符组），并将其组合成新的字符组，字符组出现的频率等于内部字符（字符组出现频率之和），如果只剩下两个字符（字符与字符组；字符组与字符组）进行步骤4；否则进入步骤3
3.将步骤2中组合后的字符与字符组重新按照出现频率降序排列，重复步骤2
4.如果按照字符组（字符）的排列顺序进行编码，按照左孩子为1，右孩子为0（或者相反）的顺序进行编码，一直到叶子节点，完成编码

实现如图1-1，选取信号符号为’aaabbiieubaddyss’，最终压缩比为1.3333。

图1-1 Huffman编码
·2关于Haming信道编码
规则：
1.设分组码(n,k)中，k个为有效信息码，n为编码后长度，r = n-k为监督码
2.r个监督码元构造r个监督关系式指出一位错码的n种可能位置，r需要满足：
3.2^r-1 =n 满足此线性分组的编码方式为汉明码----能够纠正1位错码

为了避免重复造轮子，并且考虑到汉明编码和解码实现较为复杂，此处引用matlab内置的函数。

·3 2ASK
思路：
·采用与模拟调制相同的方法，直接将信号与载波相乘
·需要注意：
周期性：一个电平周期应该是载波信号周期的整数倍

问题及解决：
·在滤波的时候，应用lowpass，结果发现总是滤不掉,查阅资料，可以设置滤波器的陡峭特性，重新设置为0.95左右，work

·在进行抽样判决的时候，起初混淆了理论中判别门限A/2（0 1 概率相同的前提下，虽然0 1 频率不相等，方便起见假设为0）A的含义，误当作了解调后的幅值，导致无论怎么加大噪声，最终误码的概率非常之低，经过排查发现了该错误。

·关于2ASK采用相干解调误码率的极限为多少：许多教材以及相关资料给出了图3-2中的图，该图表明误码率的极限应该为0.5，即随机“瞎猜”也应该为0.5，这与直观的认知是相符的。

·4 2FSK
调制：
·2FSK可以视为2个不同载频的2ASK叠加，故仍然可以采用模拟调频法，不过这里需要对电平范围进行逻辑转换
·原电平∈{1,0}，可以将原电平整体逻辑取反，然后对取反前后的信号分别做2ask在w1,w2的调制，再叠加

解调：（包络检波法较为简单，此处采用相干解调法）
分别通过两路带通滤波器，然后经过相乘器，然后低通滤波，抽样判决
流程：

  |---w1 带通滤波器 --->相乘器(载波w1)---低通-----|	

S2fsk -| |----抽样判决----out
|— w2 带通滤波器 —>相乘器(载波w2)—低通----|

2FSK仿真省略。。。。。。

·5 误码率测试问题及解决
误码率分辨精度过低：
由于所给信号长度有限，因此误码率的分辨精度（这里我指错一个码元引起的误码率变化）不高，即每错一个代码，就会引起误码率波动很大，导致最终做出的图有很多“尖刺”，如图，这是由于添加高斯噪声时的随机性引起的，

解决方法：
1.可以在每一个信噪比点上，多次重复实验，取误差率的平均值，相当于在每一个点上进行均值滤波，可以有效抵抗误码率分辨精度不高以及awgn函数随机性引起的“噪声”干扰。
2.可以将仿真中psr的仿真步长设置的小一些，取得最终的误差率为相邻几个误差率的平均值，也相当于一个均值滤波。
最佳判别门限：
由于仿真时，代码数量少，导致01出现的频率并不相同，经过计算发现这种误差不可忽略。因此重新计算判别门限，而非简单取a/2作为判别门限。


最终仿真结果如图：

参考：
https://blog.csdn.net/qq_45654781/article/details/106579495
https://blog.csdn.net/hxxjxw/article/details/82628565

代码

%% 信号初始化，为了与真实情况对照，适当增加了字符，调整0-1比例大致为0.5
close all;clear;clc;
signal = 'aaabbiieubaddsssyusongsong';% 为提高分辨率精度，又添加了几个字符
%% 信号Huffman编码 + haming(7,4)编码
[dict,huffman_code,H_ratio,num] = En_Huffuman(signal);fprintf('霍夫曼编码完成\n 压缩比为:%f\n',H_ratio);global haming_code;
haming_code = encode(huffman_code,7,4);fprintf('汉明编码完成\n 编码效率为:%f\n',4/7);
%% 2ASK调制与解调 --- 时间大约0.5min左右
tic;[err,BStar,reHaming_code,snr] = deCode(huffman_code,num);toc;
figure(1);draw(err,snr);

function [reHaming_code,bStar] = ASK(snr,lnP)
% 功能说明:实现汉明码载波调制+信道传输+相干解调+抽样判决+输出单次误码率
% input:
% snr:信噪比,单位为dB
% lnp:计算最佳判别门限
% 
% output:
% reHaming_code:解调后的汉明码
% bStar:最佳判别门限

global haming_code;
reHaming_code=zeros(size(haming_code));
numCode = length(haming_code);  % 码元个数
singleT = 100;                  % 单个信号,每秒钟采样数 fs = 100Hz
numTime = singleT*numCode;      % 总的信号序列数目
t=linspace(0,100,numTime);      % 生成时间序列，方便做图
fc=10000;                        % 载波频率
st=t;                           % 相当于采样后的信号
% 信号离散化--->调制
for n=1:numCode
    if haming_code(n)==0
        for m=numTime/numCode*(n-1)+1:numTime/numCode*n
            st(m)=0;
        end
    else
        for m=numTime/numCode*(n-1)+1:numTime/numCode*n
            st(m)=1;
        end
    end
end

% 载波调制 + 信道加噪
s = cos(2*pi*fc*t);% 载波信号
s_2ask = st.*s;
s_2ask = awgn(s_2ask,snr);

% 相干解调
ss_2ask = s_2ask.*s;
ss_2ask = lowpass(ss_2ask,0.95,singleT,'ImpulseResponse','iir','Steepness',0.99);%低通滤波

% 最佳判别门限计算
% a = max(ss_2ask);
% a = min(1,a);
a = 1;
bStar = a/2+a/2/(10^(snr/20))*lnP;

% 抽样判决
for m=0:numCode-1                    
    if ss_2ask(1,m*singleT+singleT/2) < bStar
        reHaming_code(m+1)=0;
    else
        reHaming_code(m+1)=1;
    end
end

end

function [err,BStar,reHaming_code,snr] = deCode(huffman_code,num)
% 功能说明:汉明解码 + 误码率分析
% input:
% huffman_code:传输之前的霍夫曼码
% num:霍夫曼代码长度，为了减少运算复杂度，作为参数传入
% output:
% err:误码率
% BStar:最佳判别门限，方便理解
% reHaming_code:解调及抽样判决后得到的汉明码
% snr:信噪比序列，方便做图

global haming_code;
snr = -6:1:18;
L = length(snr);
err = zeros(1,L);
BStar = zeros(1,L);
epoch = 16;

% 便于计算最佳判别门限计算参数
lnP1 = sum(haming_code)/length(haming_code);
lnP = log((1-lnP1)/lnP1);
% numHaming = length(haming_code);%调试用

for i = 1:L
    for j = 1:epoch
        [reHaming_code,bStar] = ASK(snr(i),lnP);
        rehuffman_code = decode(reHaming_code,7,4);
        Rehuffman_code = rehuffman_code(1,1:num);
        err(i) = err(i) + sum(Rehuffman_code~=huffman_code)/(num);
    end
    
    BStar(i) = bStar;
    err(i) = err(i)/epoch;
end
% 对err进行移动平均滤波，进一步消除偶然误差（由于码数少，导致单位误差率大）
windowSize = 4;
b = (1/windowSize)*ones(1,windowSize);
err = filter(b,1,err);

err = err(1,windowSize+1:end);
snr = snr(1,windowSize+1:end);

end

function draw(err,snr)
% 做出仿真误差率曲线以及理论曲线

% 仿真误差率
% Err = log10(err);
% plot(snr,Err);
semilogy(snr,err);
xlabel('r/dB');ylabel('Pe');

title('2ASK数字调制系统误码率与信噪比的关系');
hold on;

% 理论误差率
r = 10.^(snr/20);
Perr = 0.5*erfc(sqrt(r/4));
% PErr = log10(Perr);
% plot(snr,PErr);
semilogy(snr,Perr);

end

function [dict,code,H_ratio,num] = En_Huffuman(signal)
% 功能说明:对signal中的字母进行霍夫曼编码 
% 参数说明：
% input:
% signal:需要转换的字符串对象，即所需要传输的信息
% output:
% dict:霍夫曼字典
% code:对信号的编码
% H_ratio: 利用Huffman编码的压缩比
% num:总码长

Signal = tabulate(signal');% Signal 为 k*3维元胞组,k为字符类数;第一列对应字符，第三列为字符对应的频率

% 计算等长码的平均码长b
k = size(Signal,1);     % k个字符采用等长码需要的b(bit)满足  2^(b) > k > 2^(b-1);
if log2(k) ~= floor(log2(k))
    b =floor(log2(k))+1;% 反解出  log2(k) <= b < log2(k) + 1
else
    b =floor(log2(k));
end

% huffman编码
p = (cell2mat(Signal(:,3))/100)';
symbols = Signal(:,1)';
dict = huffmandict(symbols,p);
code = huffmanenco(signal,dict);

% 计算压缩后的平均码长 L
L = 0;
for i = 1 : k
    L = L + size(cell2mat(dict(i,2)))*p(i);
end

% 总码长，用于计算误码率
num = length(code);% 提高代码效率，预先计算代码长度

% 计算压缩比
H_ratio = b/L(2);

end