QPSK仿真实现要点

基于MATLAB的QPSK仿真设计与实现

一、摘要

本次方向设计根据当今现代通信技术的发展，对QPSK信号的工作原理进行了仿真分析。并结合调制解调的基本性能和通信原理的基础知识，利用MATLAB 仿真工具设计出一个QPSK仿真程序，以衡量QPSK在理想信道、高斯白噪声信道和先通过瑞利衰落信道再通过高斯信道三种方式的调制解调得到的功率谱密度、噪声曲线、星座图及误码性能,并对仿真结果进行了分析。

关键字：MATLAB仿真；QPSK调制；QPSK解调；误码率；信噪比

Based on the modern communication technology, design of oriented major has implemented a simulated analysis in regard to the principle of QPSK signal. Associating with the performance of design of oriented major and underpinning communication knowledge, a QPSK simulate program was implemented and analysed in the aim of measuring its power spectral density,curve of noise,constellation and bit error performance under ideal channel, white Gaussian noise channel and via Rayleigh fading channel and white Gaussian noise channel respectively。

Key words：MATLAB Simulate; QPSK Modulation; QPSK Demodulation; Error Rate; Signal to Noise Ratio

二、设计目的和意义

近年来，软件无线电作为解决通信体制兼容性问题的重要方法受到各方面的注意。它的中心思想是在通用的硬件平台上，用软件来实现各种功能，包括调制解调类型、数据格式、通信协议等。通过软件的增加、修改或升级就可以实现新的功能，充分体现了体制的灵活性、可扩展性等。其中软件的增加、高频谱效率的调制解调模块是移动通信系统的关键技术，它的软件化也是实现软件无线电的重要环节。

通过完成设计内容，复习QPSK调制解调的基本原理，同时也要复习通信系统的主要组成部分，了解调制解调方式中最基础的方法。了解QPSK的实现方法及数学原理。并对“通信”这个概念有个整体的理解，学习数字调制中误码率测试的标准及计算方法。同时还要复习随机信号中时域用自相关函数，频域用功率谱密度来描述平稳随机过程的特性等基础知识，来理解高斯信道中噪声的表示方法，以便在编程中使用。

理解QPSK调制解调的基本原理，并使用MATLAB编程实现QPSK信号在高斯信道和瑞利衰落信道下传输，以及该方式的误码率测试。复习MATLAB编程的基础知识和编程的常用算法以及使用MATLAB仿真系统的注意事项，并锻炼自己的编程能力，通过编程完成QPSK调制解调系统的仿真，以及误码率测试，并得出响应波形。在完成要求任务的条件下，尝试优化程序。

三、设计原理

四相相移键控信号简称“QPSK”。它分为绝对相移和相对相移两种。由于绝对相移方式存

在相位模糊问题，所以在实际中主要采用相对移相方式QDPSK。它具有一系列独特的优点，目前已经广泛应用于无线通信中，成为现代通信中一种十分重要的调制解调方式。在数字信号的调制方式中QPSK四相移键控是目前最常用的一种卫星数字信号调制方式,它具有较高的频谱利用率、较强的抗干扰性、在电路上实现也较为简单。

数字相位调制PSK是角度调制、恒定幅度数字调制的一种方式，通过改变发送波的相位来实现，除了其输入信号是数字信号以及输出的相位受限制以外，PSK与传统的相位调制相似。

QPSK信号的正弦载波有4个可能的离散相位状态，每个载波相位携带2个二进制信号。

QPSK利用载波的四种不同相位来表征数字信息。因此，对于输入的二进制数字序列应该进行分组，将每两个比特编为一组；然后用四种不同的载波相位去表征它们。例如，若输入二进制数字信息序列为10110100，则可将它们分成10,11,01,00，然后用四种不同的相位来分别表示它们。由于每一种载波相位代表2个比特信息，故每个四进制码元又被称为双比特码。

四相相移调制是利用载波的四种不同相位差来表征输入的数字信息，是四进制移相键控。QPSK是在M=4时的调相技术，它规定了四种载波相位，分别为45°，135°，225°，315°，调制器输入的数据是二进制数字序列，为了能和四进制的载波相位配合起来，则需要把二进制数据变换为四进制数据，这就是说需要把二进制数字序列中每两个比特分成一组，共有四种组合，即00，01，10，11，其中每一组称为双比特码元。每一个双比特码元是由两位二进制信息比特组成，它们分别代表四进制四个符号中的一个符号。QPSK 中每次调制可传输2个信息比特，这些信息比特是通过载波的四种相位来传递的。解调器根据星座图及接收到的载波信号的相位来判断发送端发送的信息比特。

数字调制用“星座图”来描述，星座图中定义了一种调制技术的两个基本参数：（1）信号分布；（2）与调制数字比特之间的映射关系。星座图中规定了星座点与传输比特间的对应关系，这种关系称为“映射”，一种调制技术的特性可由信号分布和映射完全定义，即可由星座图来完全定义。

同时QPSK信号可以看作两个载波正交2PSK信号的合成，下图表示QPSK正交调制器。

图1、QPSK调制系统原理图

由QPSK信号的调制可知，对它的解调可以采用与2PSK信号类似的解调方法进行解调。解调原理图如下所示，同相支路和正交支路分别采用相干解调方式解调，得到()

I t和()

Q t，经过抽样判决和并/串交换器，将上下支路得到的并行数据恢复成串行数据。

原理分析：

基本原理及系统结构

连续相位QPSK可表示为

S QPSK（t）=A

其中，为随时间连续变化的相位，f c为未调载波频率，A为已调信号幅度，由2FSK 信号正交条件可知，最小频差为

这里f1，f2分别为2FSK信号的两个频率，T S为信号码元间隔，T b为二进制信息的间隔。此时有

f c=1/2(f1+f2)

QPSK与二进制PSK一样，传输信号包含的信息都存在于相位中。的别的载波相位取四个等间隔值之一，如л/4, 3л/4,5л/4,和7л/4。相应的，可将发射信号定义为

其中，i＝1，2，2，4；E为发射信号的每个符号的能量，T为符号持续时间，载波频率f 等于nc/T，nc为固定整数。每一个可能的相位值对应于一个特定的二位组。例如，可用前述

的一组相位值来表示格雷码的一组二位组：10，00，01，11。

下面介绍QPSK 信号的产生和检测。如果a 为典型的QPSK 发射机框图。输入的二进制数据序列

首先被不归零（NRZ ）电平编码转换器转换为极性形式，即负号1和0示。接着，该二进制波形被分接器分成两个分别由输入序列的奇数位偶数位组成的彼此独立的二进制波形，这两个二进制波形分别用a1（t ），和a2（t ）表示。容易注意到，在任何一信号时间间隔内a1（t ），和a2（t ）的幅度恰好分别等于Si1和 Si2，即由发送的二位组决定。这两个二进制波形a1（t ），和a2（t ）被用来调制一对正交载波或者说正交基本函数：φ1（t ）

)c f t π，φ2（t )c f t π。这样就得到一对二进制PSK 信号。φ1（t ）和

φ2（t ）的正交性使这两个信号可以被独立地检测。最后，将这两个二进制PSK 信号相加，从而得期望的QPSK 。

图2、QPSK 解调系统原理图

四、 详细设计步骤

实验内容

1.构建一个理想信道基本QPSK 仿真系统,要求仿真结果有

a.基带输入波形及其功率谱

b. QPSK 信号及其功率谱

c. QPSK 信号星座图

2.构建一个在AWGN （高斯白噪声）信道条件下的QPSK 仿真系统，要求仿真结果有

a. QPSK 信号及其功率谱

b. QPSK 信号星座图

c. 高斯白噪声信道条件下的误码性能以及高斯白噪声的理论曲线，要求所有误码性能

曲线在同一坐标比例下绘制

3．构建一个先经过Rayleigh（瑞利衰落信道），再通过AWGN（高斯白噪声）信道条件下的条件下的QPSK仿真系统，要求仿真结果有

a. QPSK信号及其功率谱

b. 通过瑞利衰落信道之前和之后的信号星座图，前后进行比较

c. 在瑞利衰落信道和在高斯白噪声条件下的误码性能曲线，并和2.c中所要求的误码

性能曲线在同一坐标比例下绘制

QPSK调制的实现：QPSK的调制部分由电平变换，串/并变换，相乘器与相加器组成。将把原始信号按1->1,0->-1变换，之后再进行穿并变换。由于输入的串行的二进制信号，要对其进行QPSK调制，要完成的就是把串行信号变换成并行信号。信号通过平衡调制，在这里可以当作一个乘法器来进行处理。再把两路信号通过加法器叠加起来则是原信号经调制后的输出波形。

QPSK解调的实现：QPSK的解调部分由相乘器，低通滤波器，抽样判决，并/串变换组成。在解调QPSK的时候，首先会将受到的信号分为相同的两路在将这两路信号分别经过乘法器得到出去的信号波形。

抽样判决时，首先要确定门限值，即如果该抽样值大于等于门限值，则为1，否则为0。

实现串/并变换时，是将原数组的奇数位赋予数组a，偶数位赋予数组b。

实验工具：MATLAB

MATLAB是一套功能强大的工程技术数值运算和系统仿真软件，它具有数值运算和系统仿真软件，它具有数值分析、矩阵运算、数字信号处理、仿真建模、系统控制盒优化等功能。MATLAB 的编程功能简单，并且很容易扩展和创造新的命令与函数。

随着通信系统复杂性的增加，传统的手工分析与电路板试验等分析设计方法已经不能适应发展的需要，通信系统计算机模拟仿真技术日益显示出其巨大的优越性。计算机仿真是根据被研究的真实系统的模型，利用计算机进行实验研究的一种方法。它具有利用模型进行仿真的一系列优点，如费用低易于进行真实系统难于实现的各种试验，以及易于实现完全相同条件下的重复试验等。MATLAB仿真软件就是分析通信系统常用的工具之一。

MATLAB作为一种功能强大的数据分析和工程计算高级语言，已被广泛应用于现代科学技术研究和工程设计的各个领域。其信号处理工具箱可以解决通信中信号变换、调制解调、滤波、频谱估计、线性系统分析等多项功能，并且能够通过图形用户界面显示结果。MATLAB包括许多标准函数，每个函数都由完成某一特定功能的代码组成，同时，MATLAB也允许用户自行写所需的函数，其扩展名为.m，称为M文件与M函数。通过流程控制和函数语句来实现特定功能，并可利用图形显示结果。

MATLAB主要优点有很多，首先编程效率高，它是一种面向科学与工程计算的高级语言，允许用数学形式的语言编写程序。因此，MATLAB语言也可通俗的称为演算纸式科学算法语言由于它编写简单，所以编程效率高，易学易懂。其次，用户使用方便。MATLAB语言是一种解释执行的语言，它灵活、方便，其调试程序手段丰富，调试速度快，需要学习时间少，扩充能力

强，MATLAB语言有丰富的库函数，在进行复制的数学运算时可以直接调用。

MATLAB语句简单，内涵丰富，它有高效方便的矩阵和数组、运算。MATLAB语言像Basic 和C语言一样规定了矩阵的算术运算符、关系运算符、逻辑运算符、条件运算符及赋值运算符，而且这些运算符大部分可以毫无改变地照搬到数组间的运算。

程序见附录

五、设计结果及分析

根据图1和图2的流程框图设计仿真程序，得出结果并且分析如下：

5.1、理想信道下的仿真，实验结果如图3所示：

图3

实验结果分析：

如图上结果显示，完成了QPSK信号在理想信道上的调制，传输，解调的过程，由于调制过程中加进了载波，因此调制信号的功率谱密度会发生变化。并且可以看出调制解调的结果没有误码。

5.2、高斯信道下的仿真，结果如图4、5所示：

图4

图5 实验结果分析：

由图4、5可以得到高斯信道下的调制信号，高斯噪声，调制输出功率谱密度曲线和QPSK 信号的星座图。

在高斯噪声的影响下，调制信号的波形发生了明显的变化，其功率谱密度函数相对于图1中的调制信号的功率谱密度只发生了微小的变化，原因在于高斯噪声是一个均值为0的白噪声，在各个频率上其功率是均匀的，因此此结果是真确的。星座图反映可接收信号早高斯噪声的影响下发生了误码，但是大部分还是保持了原来的特性。

5.3先通过瑞利衰落信道再通过高斯信道的仿真。实验结果如图6、7、8所示：

图6

图7

图8

实验结果分析：

由图7可以得到瑞利衰落信道前后的星座图，调制信号的曲线图及其功率谱密度。最后图8显示的是高斯信道和瑞利衰落信道的误码率对比。由图可知瑞利衰落信道下的误码率比高斯信道下的误码率高。

至此，仿真实验就全部完成。

六、总结

由此次仿真实验可知，由于高斯信道和瑞利信道的影响，波形发生了明显的变化，功率谱密度也发生变化，星座图可反映在噪声影响下产生的误码。从仿真解调前与调制后的图对比可知仿真正确。通过图8可知瑞利信道下的误码率和误比特率明显高于高斯信道，并与理论值基本符合。

随着通信事业的发展，通信系统的设计也会越来越复杂，通过计算机的仿真，可以大大地降低通信过程实验成本。本文设计出一个QPSK仿真模型，以衡量QPSK在理想信道，高斯白噪声信道，以及先通过瑞利再通过高斯信道的性能，通过仿真，可以更好地了解QPSK系统的工作原理，而且为硬件的研制提供一定的参考作用。

七、体会

这次的通信专业方向设计让我把通信原理的一些内容又重新复习了一遍，中间出了很多问

题，比如用matlab7.0版本就不能把几个图整合到一起，也就是不能使用subplot函数，会出现Undefined function or variable "hanalysisparammenu"的错误.而我后来又换成了matlab7.8就可以实现该功能。我把程序发给老师后老师狠狠的教育了我，不应该依靠高版本的优化功能把在低版本不能实现的功能强制执行，所以我还有很多地方需要学习，包括对误码率进行计算的两个函数我也是借鉴的网上的程序。这次方向设计让我学会了很多，也认识到了自己还有很多方面的欠缺。

八、参考文献

[1] 李人厚、张平安.精通MATLAB，西安交通大学1997.3

[2] 肖明波，通信系统仿真原理与无线应用机械工业出版社

[3] 万永革编著，通信系统仿真原理与无线应用，北京科学出版社

[4] 郭文彬、桑林，通信原理—基于MATLAB的计算机仿真，北京邮电大学出版社，2006.1

[5]网上资料

附件：

主程序：

T=1; % 基带信号宽度，也就是频率

fc=10/T; % 载波频率

ml=2; % 调制信号类型的一个标志位

nb=100; % 传输的比特数

delta_T=T/200; % 采样间隔

fs=1/delta_T; % 采样频率

SNR=0; % 信噪比

t=0:delta_T:nb*T-delta_T; % 限定t的取值范围

N=length(t); % 采样数

%------------------------- 调制部分

% 基带信号的产生

data=randn(1,nb)>0.5;

datanrz=data.*2-1;

data1=zeros(1,nb/delta_T);

for q=1:nb

data1((q-1)/delta_T+1:q/delta_T)=datanrz(q);

end

% 串并转换，将奇偶位数据分开

idata=datanrz(1:ml:(nb-1));

qdata=datanrz(2:ml:nb);

% QPSK信号的调制

ich=zeros(1,nb/delta_T/2);

for i=1:nb/2

ich((i-1)/delta_T+1:i/delta_T)=idata(i);

end

for ii=1:N/2

a(ii)=sqrt(2/T)*cos(2*pi*fc*t(ii));

end

idata1=ich.*a;

qch=zeros(1,nb/2/delta_T);

for j1=1:nb/2

qch((j1-1)/delta_T+1:j1/delta_T)=qdata(j1);

end

for jj=1:N/2

b(jj)=sqrt(2/T)*sin(2*pi*fc*t(jj));

end

qdata1=qch.*b;

s=idata1+qdata1;

ss=abs(fft(s));

%---------------------------- 瑞利衰落信道和高斯信道

% 瑞利衰落信道

ray_ich=raylrnd(0.8,1,nb/2/delta_T);

ray_qch=raylrnd(0.8,1,nb/2/delta_T);

Ray_idata=idata1.*ray_ich;

Ray_qdata=qdata1.*ray_qch;

Ray_s=Ray_idata+Ray_qdata;

% 高斯信道

s1=awgn(s,SNR);

s11=abs(fft(s1));

s111=s1-s;

Awgn_s=awgn(Ray_s,SNR); % 通过高斯信道再通过瑞利衰落信道%----------------------------- QPSK 解调部分

% 解调部分（高斯信道）

idata2=s1.*a;

qdata2=s1.*b;

idata3=zeros(1,nb/2);

qdata3=zeros(1,nb/2);

% 抽样判决的过程，与0作比较，data>=0,则置1，否则置0

for n=1:nb/2

if sum(idata2((n-1)/delta_T+1:n/delta_T))>=0

idata3(n)=1;

else idata3(n)=0;

end

if sum(qdata2((n-1)/delta_T+1:n/delta_T))>=0

qdata3(n)=1;

else qdata3(n)=0;

end

end

% 为了显示星座图,将信号进行处理

idata4=zeros(1,nb/2);

qdata4=zeros(1,nb/2);

for n=1:nb/2

Awgn_ichsum(n)=sum(idata2((n-1)/delta_T+1:n/delta_T))*delta_T; if Awgn_ichsum(n)>=0

idata4(n)=1;

else idata4(n)=0;

end

Awgn_qchsum(n)=sum(qdata2((n-1)/delta_T+1:n/delta_T))*delta_T; if Awgn_qchsum(n)>=0

qdata4(n)=1;

else qdata4(n)=0;

end

end

% 将判决之后的数据存放进数组

demodata=zeros(1,nb);

demodata(1:ml:(nb-1))=idata3;

demodata(2:ml:nb)=qdata3;

%为了显示，将它变成波形信号（即传输一个1代表单位宽度的高电平）demodata1=zeros(1,nb/delta_T);

for q=1:nb

demodata1((q-1)/delta_T+1:q/delta_T)=demodata(q);

end

% 累计误码数

% abs(demodata-data)求接收端和发射端

% 数据差的绝对值，累计之后就是误码个数

Awgn_num_BER=sum(abs(demodata-data))

% 解调部分（瑞利+高斯）

Ray_idata2=Ray_s.*a;

Ray_qdata2=Ray_s.*b;

% 为了显示星座图,将信号进行处理

Ray_idata4=zeros(1,nb/2);

Ray_qdata4=zeros(1,nb/2);

for n=1:nb/2

Ray_ichsum(n)=sum(idata2((n-1)/delta_T+1:n/delta_T))*delta_T;

if Ray_ichsum(n)>=0

Ray_idata4(n)=1;

else Ray_idata4(n)=0;

end

Ray_qchsum(n)=sum(qdata2((n-1)/delta_T+1:n/delta_T))*delta_T; if Ray_qchsum(n)>=0

Ray_qdata4(n)=1;

else Ray_qdata4(n)=0;

end

end

% 将判决之后的数据存放进数组

Ray_demodata=zeros(1,nb);

Ray_demodata(1:ml:(nb-1))=Ray_idata4;

Ray_demodata(2:ml:nb)=Ray_qdata4;

%为了显示，将它变成波形信号（即传输一个1代表单位宽度的高电平）Ray_demodata1=zeros(1,nb/delta_T);

for q=1:nb

Ray_demodata1((q-1)/delta_T+1:q/delta_T)=Ray_demodata(q); end

% 累计误码数

% abs(demodata-data)求接收端和发射端

% 数据差的绝对值，累计之后就是误码个数

Ray_num_BER=sum(abs(Ray_demodata-data))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% 误码率计算

%% 调用了cm_sm32();和cm_sm33(）函数

%%声明：函数声明在另外俩个M文件中

%%作用： cm_sm32()用于瑞利信道误码率的计算

%% cm_sm33()用于高斯信道误码率的计算

%% ecoh on/off 作用在于决定是否显示指令内容%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

SNRindB1=0:1:6;

SNRindB2=0:0.1:6;

% 瑞利衰落信道

for i=1:length(SNRindB1),

[pb,ps]=cm_sm32(SNRindB1(i)); % 比特误码率

smld_bit_ray_err_prb(i)=pb;

smld_symbol_ray_err_prb(i)=ps;

disp([ps,pb]);

echo off;

end;

% 高斯信道

echo on;

for i=1:length(SNRindB1),

[pb1,ps1]=cm_sm33(SNRindB1(i));

smld_bit_awgn_err_prb(i)=pb1;

smld_symbol_awgn_err_prb(i)=ps1;

disp([ps1,pb1]);

echo off;

end;

% 理论曲线

echo on;

for i=1:length(SNRindB2),

SNR=exp(SNRindB2(i)*log(10)/10); theo_err_awgn_prb(i)=0.5*erfc(sqrt(SNR)); theo_err_ray_prb(i)=0.5*(1-1/sqrt(1+1/SNR)); echo off;

end;

h = spectrum.welch;

%---------------------------- 输出显示部分

% 第一部分（理想）

figure(1)

subplot(3,2,1);

plot(data0),title('基带信号');

axis([0 20000 -2 2]);

subplot(3,2,2);

psd(h,data1,'fs',fs),title('基带信号功率谱密度'); subplot(3,2,3);

plot(s),title('调制信号');

axis([0 500 -3 3]);

subplot(3,2,4);

psd(h,s,'fs',fs),title('调制信号功率谱密度');

subplot(3,2,5);

plot(demodata1),title('解调输出');

axis([0 20000 -2 2]);

subplot(3,2,6);

psd(h,demodata1,'fs',fs),title('解调输出功率谱密度');

% 通过高斯信道

figure(2)

subplot(2,2,1);

plot(s1),title('调制信号(Awgn)');

axis([0 500 -5 5]);

subplot(2,2,2);

psd(h,s1,'fs',fs),title('调制信号功率谱密度(Awgn)'); subplot(2,2,3);

plot(s111),title('高斯噪声曲线');

axis([0 2000 -5 5]);

figure(3)

for i=1:nb/2

plot(idata(i),qdata(i),'r+'),title('QPSK信号星座图（Awgn）');hold on; axis([-2 2 -2 2]);

plot(Awgn_ichsum(i),Awgn_qchsum(i),'*');hold on;

legend('理论值（发射端）','实际值（接收端）');

end

%通过高斯信道再通过瑞利衰落信道

figure(4)

subplot(2,2,1)

plot(Ray_s),title('调制信号(Ray+Awgn)');

axis([0 500 -5 5]);

subplot(2,2,2);

psd(h,Ray_s,'fs',fs),title('调制信号功率谱密度(Ray)');

figure(5)

for i=1:nb/2

plot(idata(i),qdata(i),'r+'),title('QPSK信号星座图（Awgn+Ray）');hold on; axis([-2 2 -2 2]);

plot(Ray_ichsum(i),Ray_qchsum(i),'*');hold on;

legend('理论值（发射端）','实际值（接收端）');

end

figure(6)

semilogy(SNRindB2,theo_err_awgn_prb,'r'),title('误码率曲线');hold on;

semilogy(SNRindB1,smld_bit_awgn_err_prb,'r*');hold on;

semilogy(SNRindB2,theo_err_ray_prb);hold on;

semilogy(SNRindB1,smld_bit_ray_err_prb,'+');

xlabel('Eb/No');ylabel('BER');

legend('理论AWGN','仿真AWGN','理论Rayleigh','仿真Rayleigh');

cm_sm32()用于瑞利信道误码率的计算：

function [pb,ps]=cm_sm32(snr_in_dB)

% [pb,ps]=cm_sm32(snr_in_dB)

% CM_SM3 finds the probability of bit error and symbol error for

% the given value of snr_in_dB, signal to noise ratio in dB.

N=100;

E=1; % energy per symbol numofsymbolerror=0;

numofbiterror=0;

counter=0;

snr=10^(snr_in_dB/10); % signal to noise ratio

sgma=sqrt(E/snr)/2; % noise variance

s00=[1 0]; s01=[0 1]; s11=[-1 0]; s10=[0 -1]; % signal mapping

% generation of the data source

while(numofbiterror<100)

for i=1:N,

temp=rand; % a uniform random variable between 0 and 1

if (temp<0.25), % with probability 1/4, source output is "00"

dsource1(i)=0; dsource2(i)=0;

elseif (temp<0.5), % with probability 1/4, source output is "01"

dsource1(i)=0; dsource2(i)=1;

elseif (temp<0.75), % with probability 1/4, source output is "10"

dsource1(i)=1; dsource2(i)=0;

else % with probability 1/4, source output is "11"

dsource1(i)=1; dsource2(i)=1;

end;

end;

% detection and the probability of error calculation

for i=1:N,

ray=raylrnd(0.8);

n=sgma*randn(1,2); % 2 normal distributed r.v with 0, variance sgma

if ((dsource1(i)==0) & (dsource2(i)==0)),

r=ray*s00+n;

elseif ((dsource1(i)==0) & (dsource2(i)==1)),

r=ray*s01+n;

elseif ((dsource1(i)==1) & (dsource2(i)==0)),

r=s10*ray+n;

else

r=s11*ray+n;

end;

% The correlation metrics are computed below

c00=dot(r,s00); c01=dot(r,s01); c10=dot(r,s10); c11=dot(r,s11);

% The decision on the ith symbol is made next

c_max=max([c00,c01,c10,c11]);

if (c00==c_max), decis1=0; decis2=0;

elseif (c01==c_max), decis1=0; decis2=1;

elseif (c10==c_max), decis1=1; decis2=0;

else decis1=1; decis2=1;

end;

% Increment the error counter, if the decision is not correct

symbolerror=0;

if (decis1~=dsource1(i)), numofbiterror=numofbiterror+1; symbolerror=1;

end;

if (decis2~=dsource2(i)), numofbiterror=numofbiterror+1; symbolerror=1;

end;

if (symbolerror==1), numofsymbolerror=numofsymbolerror+1;

end;

end

counter=counter+1;

end

ps=numofsymbolerror/(N*counter); % since there are totally N symbols

pb=numofbiterror/(2*N*counter); % since 2N bits are transmitted

cm_sm33()用于高斯信道误码率的计算：

function [pb1,ps1]=cm_sm33(snr_in_dB)

% [pb,ps]=cm_sm32(snr_in_dB)

% CM_SM3 finds the probability of bit error and symbol error for

% the given value of snr_in_dB, signal to noise ratio in dB.

N=100;

E=1; % energy per symbol

snr=10^(snr_in_dB/10); % signal to noise ratio

sgma=sqrt(E/snr)/2; % noise variance

s00=[1 0]; s01=[0 1]; s11=[-1 0]; s10=[0 -1]; % signal mapping

% generation of the data source

numofsymbolerror=0;

numofbiterror=0;

counter=0;

while(numofbiterror<100)

for i=1:N,

temp=rand; % a uniform random variable between 0 and 1

if (temp<0.25), % with probability 1/4, source output is "00"

dsource1(i)=0; dsource2(i)=0;

elseif (temp<0.5), % with probability 1/4, source output is "01"

dsource1(i)=0; dsource2(i)=1;

elseif (temp<0.75), % with probability 1/4, source output is "10"

dsource1(i)=1; dsource2(i)=0;

else % with probability 1/4, source output is "11"

dsource1(i)=1; dsource2(i)=1;

end;

end;

% detection and the probability of error calculation

for i=1:N,

% the received signal at the detection, for the ith symbol,is:

n=sgma*randn(1,2); % 2 normal distributed r.v with 0, variance

sgma

if ((dsource1(i)==0) & (dsource2(i)==0)),

r=s00+n;

elseif ((dsource1(i)==0) & (dsource2(i)==1)),

r=s01+n;

elseif ((dsource1(i)==1) & (dsource2(i)==0)),

r=s10+n;

else

r=s11+n;

end;

% The correlation metrics are computed below

c00=dot(r,s00); c01=dot(r,s01); c10=dot(r,s10); c11=dot(r,s11);

% The decision on the ith symbol is made next

c_max=max([c00,c01,c10,c11]);

if (c00==c_max), decis1=0; decis2=0;

elseif (c01==c_max), decis1=0; decis2=1;

elseif (c10==c_max), decis1=1; decis2=0;

else decis1=1; decis2=1;

end;

% Increment the error counter, if the decision is not correct

symbolerror=0;

if (decis1~=dsource1(i)), numofbiterror=numofbiterror+1; symbolerror=1;

end;

if (decis2~=dsource2(i)), numofbiterror=numofbiterror+1; symbolerror=1;

end;

if (symbolerror==1), numofsymbolerror=numofsymbolerror+1;

end;

end

counter=counter+1;

end

ps1=numofsymbolerror/(N*counter); % since there are totally N symbols

pb1=numofbiterror/(2*N*counter); % since 2N bits are transmitted