AWGN信道中BPSK调制系统的BER仿真计算解读

序号（学号）：
学生实验报告书
2014 年 4 月27 日
实验一：AWGN 信道中BPSK 调制系统的
BER 仿真计算
一、实验目的
1．掌握二相BPSK 调制的工作原理
2．掌握利用MATLAB 进行误比特率测试BER 的方法
3．掌握AWGN信道中BPSK调制系统的BER仿真计算方法
二．实验内容
利用仿真程序在MATLAB 环境下完成AWGN信道中BPSK调制系统的BER仿真计算，得到仿真结果，写出实验小结，完成实验报告。

三．实验仪器：
计算机
matlab软件
四、实验原理
在数字领域进行的最多的仿真任务是进行调制解调器的误比特率测试，在相同的条件下进行比较的话，接收器的误比特率性能是一个十分重要的指标。

误比特率的测试需要一个发送器、一个接收器和一条信道。

首先需要产生一个长的随机比特序列作为发送器的输入，发送器将这些比特调制成某种形式的信号以便传送到仿真信道，我们在传输信道上加上一定的可调制噪声，这些噪声信号会变成接收器的输入，接收器解调信号然后恢复比特序列，最后比较接收到的比特和传送的比特并计算错误。

误比特率性能常能描述成二维图像。

纵
坐标是归一化的信噪比，即每个比特的能量除以噪声的单边功率谱密度，单位为分贝。

横坐标为误比特率，没有量纲。

五．实验步骤
①运行发生器：通过发送器将伪随机序列变成数字化的调制信号。

②设定信噪比：假定 SNR 为 m dB,则 Eb/N0=10,用 MATLAB 假设SNR 单位为分贝。

③确定Eb ④计算N0 ⑤计算噪声的方差σ n ⑥产生噪声：因为噪声具有零均值，所以其功率和方差相等。

我们产生一个和信号长度相同的噪声向量，且该向量方差为σ n 。

⑦加上噪声，运行接收器⑧确定时间延迟⑨产生误差向量⑩统计错误比特：误差向量“err”中的每一个非零元素对应着一个错误的比特。

最后计算误比特率 BER：每运行一次误比特率仿真，就需要传输和接收固定数量的比特，然后确定接收到的比特中有多少错误的。

使用 MATLAB 计算BER: ber=te/length(tx)。

六．实验结果及分析
MATLAB程序：
%Simulation of bpskAWGN
Max_SNR=10;
N_trials=1000;
N=200;
Eb=1;
ber_m=0;
for trial=1:1:N_trials
trial
msg=round(rand(1,N)); % 1, 0 sequence
s=1-msg.*2; %0-->1, 1-->1
n=randn(1,N)+j.*randn(1,N); %generate guass white noise
ber_v=[];
for snr_dB=1:2:Max_SNR
snr=10.^(snr_dB./10); %snr(db)-->snr(decimal)
N0=Eb./snr;
sgma=sqrt(N0./2);
y=sqrt(Eb).*s+sgma.*n;
y1=sign(real(y));
y2=(1-y1)./2; % 1, 0 sequence
error=sum(abs( msg- y2)); %error bits
ber_snr=error./N; %ber
ber_v=[ber_v,ber_snr];
end%for snr
ber_m=ber_m+ber_v;
end
ber=ber_m./N_trials;
ber_theory=[];
for snr_db=1:2:Max_SNR
snr=10.^(snr_db./10);
snr_1=Qfunct(sqrt(2*snr));
ber_theory=[ber_theory,snr_1];
end
i=1:2:Max_SNR;
semilogy(i,ber,'-r',i,ber_theory,'*b');
xlabel('E_b/N_0 (dB)')
ylabel('BER')
legend('Monte Carlo', 'Theoretic')
程序分析：
做1000次试验，每次试验取200个抽样点，求出每次试验的误比特率，然后对1000次试验的误比特率取平均值，即得仿真误比特率ber，然后将此误比特率与理论值ber_theory进行比较。

程序运行结果如图3所示。

结果分析：
由图3可知，仿真结果与理论值基本相符，只是在信噪比较大时，仿真误码率与理论值存在误差，这是因为在高信噪比情况下，误比特率会呈下降趋势，若要仿真实际的误比特率，就需要进行更多次试验，而且每次试验取的抽样点数也应该增加。

%Simulation of qpskAWGN
N_trials=1000;
N_number=100;
N_snr=10;
Es=1;
BER_m=0;
SER_m=0;
for trials=1:N_trials;
trials
s10=round(rand(1,N_number));
S=(s10*2-1)./sqrt(2);
S1=S(1:2:N_number);
S2=S(2:2:N_number);
Sc=S1+j.*S2; %generate qpsk signal
niose=randn(1,N_number/2)+j.*randn(1,N_number/2); %generate noise SER_v=[]; %Symbol error rate
BER_v=[]; %Bit error rate
for snr_db=0:1:N_snr;
sgma=(1/2)*sqrt(10.^(-snr_db./10));
Y=Sc+sgma.*niose;
Y_r=sign(real(Y))./sqrt(2);
Y_i=sign(imag(Y))./sqrt(2);
Y_bit=[];
for k=1:length(Y_r);
Y_bit=[Y_bit,[Y_r(k),Y_i(k)]];
end;
Y_symbol=Y_r+j*Y_i;
X_b=S-Y_bit;
X_s=Sc-Y_symbol;
ber_snr=0;
for k=1:N_number
if X_b(k)~=0;
ber_snr=ber_snr+1;
end;
end;
ser_snr=0;
for k=1:N_number/2;
if X_s(k)~=0;
ser_snr=ser_snr+1;
end;
end;
BER_v=[BER_v,ber_snr./N_number];
SER_v=[SER_v,ser_snr./(N_number./2)];
end;%for SNR
BER_m=BER_m+BER_v;
SER_m=SER_m+SER_v;
end% for trials
BER=BER_m./N_trials;
SER=SER_m./N_trials;
BER_T=[];
SER_T=[];
for snr_db=0:1:N_snr;
snr=10.^(snr_db./10);
BER_THEORY=Qfunct(sqrt(2.*snr));
SER_THEORY=1-(1-(1/2).*erfc(sqrt(snr))).^2;
BER_T=[BER_T,BER_THEORY];
SER_T=[SER_T,SER_THEORY];
end;
figure
i=0:1:N_snr;
semilogy(i,BER,'-r',i,BER_T,'*b');
legend('BER-simulation','BER-theory');
xlabel('Eb/N0(db)');
ylabel('BER');
figure
i=0:1:N_snr;
semilogy(i,SER,'-g',i,SER_T,'*y');
legend('SER-simulation','SER-theory');
xlabel('E_b/N_0(db)');
ylabel('SER');
程序分析：
与bpskAWGN类似，做1000次试验，每次试验取100个点，然后每两个点形成一个qpsk符号，首先计算仿真误比特率和误码率，然后与理论值进行比较。

程序运行结果如图4和图5所示，其中图4是ber图，图5是ser图。

结果分析：
由图4和图5可知，在低信噪比的情况下，BER和SER的仿真结果与理论值完全相符，在高信噪比情况下，两者存在差异，这是因为随着信噪比的增加，误比特率和误码率也会增加，这就需要更多的试验次数和更多的抽样点数。

将BER和SER仿真结果进行比较可知，在相同信噪比的情况下，BER比SER小，这是因为只要一个比特错，相应的qpsk 符号就会出错，而一个qpsk要正确接收，就必须保证它的两个比特全部正确接收。

图3 SER仿真结果与理论值
实验二：无线通信信道建模
一、实验目的
1.掌握无线移动通信信道的特点，根据相关模型进行仿真。

2.掌握MATLAB 语言对上述参数进行仿真。

二．实验内容
利用仿真程序在MATLAB 环境下完成移动信道建模的仿真分析，得到仿真结果，写出实验小结，完成实验报告。

三．实验仪器
计算机
matlab软件
四、实验原理
当移动台在一个较小的范围（小于20个工作波长）运动时，引起接收信号的幅度、相位和到达角等的快速变化，这种变化称为小尺度衰落。

典型的小尺度衰落有Rayleigh、Rician衰落，因为当信号在传播过程中经过许多反射路径后，接收到的信号幅度可以用Rayleigh
或Rician概率密度函数来描述。

在接受信号有直达信号LOS的情况下，幅度的衰落呈现Rician分布，而当在接收端没有直达信号的情况下，幅度的衰落呈现Rayleigh分布。

采用小尺度衰落模型的信道，衰落幅度是服从Rician或Rayleigh分布的随机变量，这些变量将会影响到接收信号的幅度和功率。

五．实验步骤
1．选择路径数
2．按均匀分布产生各条路径的延迟
3．按功率时延谱确定对应的各径的功率
4．按Jake 模型产生各径的瑞利衰落系数
5．对瑞利衰落系数进行统计分析并与理论值相比较
说明：1．路径数目2-4 自己确定，或采用某个国际标准 2．每条路径时间延迟满足（0，Tmax）范围内均匀分布，Tmax 为自己选择的最大采样步长数200-600 间比较合适，或采用国际标准 3．功率可以按时延迟谱求得，也可用国际标准测量值。

功率延迟谱：①若采用等功率分配产生功率：P i =P t /M；②采用指数分布的功率延迟谱产生功率：P=1/6*exp(-t/6)
六．实验结果及分析
MATLAB程序：
% Simulation Of Jakes Model
clear all;
f_max = 30;
M = 8; % # of low frequency oscillators - 1
N = 4*M+2;
Ts=1.024e-04;
sq = 2/sqrt(N);
sigma = 1/sqrt(2);
theta = 0; % Fixed Phase
count = 0;
t0=0.001;
for t = 0:Ts:0.5 % Varying time
count = count + 1
g(count) = 0;
for n = 1 : M+1,
if n <= M
c_q(count,n) = 2*sigma*sin(pi*n/M); % Gain associated with quadrature component
c_i(count,n) = 2*sigma*cos(pi*n/M); % Gain associated with inphase component
f_i(count,n) = f_max*cos(2*pi*n/N); % Discrete doppler frequencies of inphase component
f_q(count,n) = f_max*cos(2*pi*n/N); % Discrete doppler frequencies of quadrature component
else
c_i(count,n) = sqrt(2)*cos(pi/4);
c_q(count,n) = sqrt(2)*sin(pi/4);
f_i(count,n) = f_max;
f_q(count,n) = f_max;
end; % end if
g_i(count,n) = c_i(count,n)*cos(2*pi*f_i(count,n)*(t-t0) + theta); % Inphase component for one oscillator
g_q(count,n) = c_q(count,n)*cos(2*pi*f_q(count,n)*(t-t0) + theta); % Quadrature component for one oscillator
end; %end n
tp(count) = sq*sum(g_i(count,1:M+1));% Total Inphase component
tp1(count) = sq*sum(g_q(count,1:M+1));% Total quadrature component
end; % end count no n again
envelope=sqrt(tp.^2+tp1.^2); %rayleigh envelope
rmsenv=sqrt(sum(envelope.^2)/count); %root mean square envelpe
[auto_i,lag_i] = xcorr(tp,'coeff') ; % Auto-correlation associated with inphase component [auto_q,lag_q] = xcorr(tp1,'coeff'); % Auto-correlation associated with quadrature component
len=length(lag_i);
[corrx2,lag2] = xcorr(tp,tp1,'coeff');% Cross Correlation between inphase and quadrature components
aa=-(len-1)/2:1:(len-1)/2;%total duration for lag
bb=(len-2001)./2;%mid ... points for drawing figures
cc=bb+1:1:bb+2001;%for getting the mid-values
dd=-1000:1:1000;
%_______________
tdd=dd*Ts;
z=2.*pi.*f_max*tdd;
sigma0=1;
T_bessel=sigma0.^2.*besselj(0,z);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%
figure;
plot(tdd,auto_i(cc),'-',tdd,T_bessel,'*');%in-phase
xlabel('t(Second)');
ylabel('Auto-correlation');
legend('In-component')
figure;
plot(tdd,auto_q(cc),'-',tdd,T_bessel,'*');%quadrature
xlabel('t(Second)');
ylabel('Auto-correlation');
legend('Q-component')
figure
co1=1:1000;
semilogy(co1*Ts,envelope(1:1000));
xlabel('t(Second)');
ylabel('Rayleigh Coef.');
%%__________________________________________
length_r=length(envelope);
%____________________________
pdf_env=zeros(1,501);
count=0;
temp=round(100.*envelope);
for k=1:length_r
if temp(k)<=500
count=count+1;
pdf_env(1,temp(k)+1)=pdf_env(1,temp(k)+1)+1;
end
end
count
pdf_env=pdf_env./count./0.01; %simulation rayleigh pdf
sgma2=0.5;
x=[0:0.01:5];
pdf_theory=(x./sgma2).*exp(-1.*x.^2./(2.*sgma2)); %theory rayleigh pdf
figure
plot(x,pdf_env,'-',x,pdf_theory,'*');
legend('Simulated','Theoretic');
xlabel('r');
ylabel('PDF of r');
程序分析：
首先选定jake模型得到参数M=8、N=4*M+2，然后计算各路的同相与正交分量的增益和多普勒频率，以及它们的表达式，该程序假设各路的延时相等，均为t0=0.001，各路相位均为0。

然后分别求同相分量和正交分量的自相关和互相关，分别作图，并分别与理论值相比较。

接着作出瑞利衰落系数图形，最后作出仿真的和理论的瑞利分布的概率密度函数，并进行比较。

图6. 瑞利衰落系数图形
结果分析：
图4和图5的自相关和互相关函数图形仿真结果与理论结果大致相同，说明jake模型很好地分离出了接收信号同相分量和正交分量，图6是瑞利衰落系数的时域图形，图7中瑞利分布概率密度函数的仿真结果的大致变化趋势与理论值相同，只不过在概率密度较大的区域两者有较大差异，这是因为抽样点数太少，无法满足理论计算的要求，而且jake模型的M、N取值也会影响概率密度函数图形的仿真结果。

实验三： CDMA通信仿真
一、实验目的
1. CDMA通信具有很多通信特点，不仅被IS-95移动通信系统使用，目前已成为3G的主要技术。

2. 通过实验：
①掌握直接序列扩频发射机与接收机的组成与仿真；
②仿真验证AWGN信道下单用户直接序列扩频系统的BER性能；
③仿真验证平坦瑞利信道下单用户直接序列扩频系统的BER性能；
④观察存在干扰用户时的系统性能变化。

二．实验内容
设计一个CDMA系统，用MATLAB进行仿真，统计BER或SER随信噪比的关系，绘出曲线。

对统计试验的结果与单用户的理论值进行比较，对仿真结果进行分析。

三．实验仪器
计算机
matlab软件
四．实验原理
仿真基带直接序列扩频系统：
1.采用BPSK或QPSK映射。

2.扩频序列可以是随机产生，可以是m序列，也可以是Gold码，长
度自选。

3.最后对BER或SER随信噪比变化画图与理论单用户的结果比较，
并对仿真结果进行分析。

五．实验步骤
1.确定用户数目、信道特征以及调制方式。

2.确定基带扩频仿真系统的原理结构图，按照框图设计一个CDMA
系统，并进行仿真。

3.用MATLAB进行仿真，统计BER或SER随信噪比的关系，绘出
曲线。

4.对统计试验的结果与单用户的理论值进行比较。

5.对仿真结果进行分析。

六．实验结果及分析
MATLAB程序：
%main_IS95_forward.m
%此函数用于IS-95前向链路系统的仿真，包括扩
%频调制，匹配滤波，RAKE接收等相关通信模块。

%仿真环境: 加性高斯白噪声信道.
%数据速率= 9600 KBps
%
clear all
close all
clc
disp('--------------start-------------------');
global Zi Zq Zs show R Gi Gq
clear j;
show = 0; %控制程序运行中的显示
SD = 0; % 选择软/硬判决接收
%-------------------主要的仿真参数设置------------------
BitRate = 9600; %比特率
ChipRate = 1228800; %码片速率
N = 184; %源数据数
MFType = 1; % 匹配滤波器类型--升余弦
R = 5;
%+++++++++++++++++++Viterbi生成多项式++++++++++++++++++
G_Vit = [1 1 1 1 0 1 0 1 1; 1 0 1 1 1 0 0 0 1];%Viterbi生成多项式矩阵
K = size(G_Vit, 2); %列数
L = size(G_Vit, 1); %行数
%++++++++++++++++++++++++++++++++++++++++++++++++++++++
%++++++++++++++++++++++Walsh矩阵++++++++++++++++++++++++ WLen = 64; %walsh码的长度
Walsh = reshape([1;0]*ones(1, WLen/2), WLen , 1); %32个1 0行
%Walsh = zeros(WLen ,1);
%++++++++++++++++++++++++++++++++++++++++++++++++++++++
%++++++++++++++++++扩频调制PN码的生成多项式++++++++++++++
%Gi = [ 1 0 1 0 0 0 1 1 1 0 1 0 0 0 0 1]';
%Gq = [ 1 0 0 1 1 1 0 0 0 1 1 1 1 0 0 1]';
Gi_ind = [15, 13, 9, 8, 7, 5, 0]'; %i路PN码生成多项式参数
Gq_ind = [15, 12, 11, 10, 6, 5, 4, 3, 0]'; %q路PN码生成多项式参数
Gi = zeros(16, 1); %16×1的0矩阵
Gi(16-Gi_ind) = ones(size(Gi_ind));%根据Gi_ind配置i路PN码生成多项式
Zi = [zeros(length(Gi)-1, 1); 1];
% I路信道PN码生成器的初始状态
Gq = zeros(16, 1); %16×1的0矩阵
Gq(16-Gq_ind) = ones(size(Gq_ind)); %根据Gq_ind配置q路PN码生成多项式Zq = [zeros(length(Gq)-1, 1); 1];
% Q路信道PN码生成器的初始状态
%++++++++++++++++++++++++++++++++++++++++++++++++++++++
%+++++++++++++++++++扰码生成多项式++++++++++++++++++++++
Gs_ind = [42, 35, 33, 31, 27, 26, 25, 22, 21, 19, 18, 17, 16, 10, 7, 6, 5, 3, 2, 1, 0]'; Gs = zeros(43, 1); %43×1的0矩阵
Gs(43-Gs_ind) = ones(size(Gs_ind)); %根据Gs_ind配置扰码生成多项式
Zs = [zeros(length(Gs)-1, 1); 1];
% 长序列生成器的初始状态
%++++++++++++++++++++++++++++++++++++++++++++++++++++++
%++++++++++++++++++++++AWGN信道++++++++++++++++++++++++ EbEc = 10*log10(ChipRate/BitRate);%处理增益
EbEcVit = 10*log10(L);
EbNo = [-1: 0.5 : 1]; %仿真信噪比范围(dB)
%EbNo = [-2 : 0.5 : -1.5];
%++++++++++++++++++++++++++++++++++++++++++++++++++++++
%------------------------------------------------------
%-------------------------主程序-------------------------
ErrorsB = []; ErrorsC = []; NN = [];
if (SD == 1) % 判断软/硬判决接收
fprintf('\n SOFT Decision Viterbi Decoder\n\n');
else
fprintf('\n HARD Decision Viterbi Decoder\n\n');
end
for i=1:length(EbNo) %根据EbNo多次运行
fprintf('\nProcessing %1.1f (dB)', EbNo(i));%输出当前EbNo值
iter = 0; ErrB = 0; ErrC = 0;
while (ErrB <300) & (iter <150)
drawnow;
%++++++++++++++++++++++发射机+++++++++++++++++++++++
TxData = (randn(N, 1)>0);%生成源数据
% 速率为19.2Kcps
[TxChips, Scrambler] = PacketBuilder(TxData, G_Vit, Gs); %产生IS-95前向链路系统的发送数据包
% 速率为1.2288Mcps
[x PN MF] = Modulator(TxChips, MFType, Walsh);%实现IS-95前向链路系统的数据调制
%++++++++++++++++++++++++++++++++++++++++++++++++++++++
%++++++++++++++++++++++++信道+++++++++++++++++++++++++++
noise = 1/sqrt(2)*sqrt(R/2)*( randn(size(x)) + j*randn(size(x)))*10^(-(EbNo(i) - EbEc)/20);%生成噪声序列
r = x+noise;%加入噪声
%++++++++++++++++++++++++++++++++++++++++++++++++++++++
%+++++++++++++++++++++++++接收机++++++++++++++++++++++++
RxSD = Demodulator(r, PN, MF, Walsh); %软判决，速率为19.2 Kcps
RxHD = (RxSD>0); % 定义接收码片的硬判决
if (SD)
[RxData Metric]= ReceiverSD(RxSD, G_Vit, Scrambler); %软判决
else
[RxData Metric]= ReceiverHD(RxHD, G_Vit, Scrambler); %硬判决
end
%++++++++++++++++++++++++++++++++++++++++++++++++++++++
if(show)
subplot(311); plot(RxSD, '-o'); title('Soft Decisions'); %软判决结果图
subplot(312); plot(xor(TxChips, RxHD), '-o'); title('Chip Errors');%RAKE接收机输入符号与发送码相比出错的码
subplot(313); plot(xor(TxData, RxData), '-o'); %硬判决接收机与发送数据相比的出错码
title(['Data Bit Errors. Metric = ', num2str(Metric)]);
pause;
end
if(mod(iter, 50)==0) %每50次保存一次
fprintf('.');
save TempResults ErrB ErrC N iter%保存结果
end
ErrB = ErrB + sum(xor(RxData, TxData));%求出错比特数
ErrC = ErrC + sum(xor(RxHD, TxChips)); %求出错码数
iter = iter+ 1;%迭代次数
end
ErrorsB = [ErrorsB; ErrB]; %存储各EbNo值下的出错比特数
ErrorsC = [ErrorsC; ErrC]; %存储各EbNo值下的出错码数
NN = [NN; N*iter]; %存储各EbNo值下的总数据码数目
save SimData*%保存当前迭代的数据
end
%+++++++++++++++++++++++++误码率计算++++++++++++++++++++++++ PerrB = ErrorsB./NN; %出错比特比例
%PerrB1 = ErrorsB1./NN1;
PerrC = ErrorsC./NN; %出错码比例
Pbpsk= 1/2*erfc(sqrt(10.^(EbNo/10))); %EbNo的余误差
PcVit= 1/2*erfc(sqrt(10.^((EbNo-EbEcVit)/10)));%EbNo-EbEcVit的余误差
Pc = 1/2*erfc(sqrt(10.^((EbNo-EbEc)/10)));%EbNo-EbEc的余误差
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ %%+++++++++++++++++++++++++性能仿真显示++++++++++++++++++++++ figure;
semilogy(EbNo(1:length(PerrB)), PerrB, 'b-*'); hold on;%信噪比误码率图
% %semilogy(EbNo(1:length(PerrB1)), PerrB1, 'k-o'); hold on;
% semilogy(EbNo(1:length(PerrC)), PerrC, 'b-o'); grid on;
% semilogy(EbNo, Pbpsk, 'b-.^');
% %semilogy(EbNo, PcVit, 'k-.x'); ylabel('BER');
% semilogy(EbNo, Pc, 'b-.x');
xlabel('信噪比/dB');
ylabel('误码率');
grid on;
% legend('Pb of System (HD)', 'Pb of System (SD)', 'Pc before Viterbi of System',
% ... 'Pb of BPSK with no Viterbi (theory)', 'Pc on Receiver (theory)');
%
% legend('Pb of System', 'Pc before Viterbi of System', ...
%'Pb of BPSK with no Viterbi (theory)',
%'Pc before Viterbi (theory)', 'Pc on Receiver (theory)');
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ disp('--------------end-------------------');
仿真结果如下：
结果分析：从图中可以看出，信噪比增大时，误码率先增大后减小，且误码率小于0.4。