武汉理工大学移动通信实验报告

管我坏不坏,又没叫你爱!
590次浏览
2020年01月01日 22:41
最佳经验
本文由作者推荐
武汉开发移动app
序号(学号):

实验报告书
实验课程名称
移动通信系统实验
开课学院
信息工程学院
指导老师姓名
学生姓名
学生专业班级
实验课程名称:移动通信系统实验
实验项目名称
AWGN信道中BPSK调制系统的BER仿真计算
实验成绩
实 验 者
专业班级
组    别
同 组 者
/
实验日期
2014年4月11日
1.掌握二相BPSK调制的工作原理;
2.掌握利用MATLAB进行误比特率测试BER的方法;
3.掌握AWGN信道中BPSK调制系统的BER仿真计算方法。
二.实验仪器
1.计算器及操作系统
软件
三.实验原理
1. 仿真概述及原理
在数字领域进行的最多的仿真任务是进行调制解调器的误比特率测试,在相同的条件下
进行比较的话,接收器的误比特率性能是一个十分重要的指标。误比特率的测试需要一个发送器、一个接收器和一条信道。首先需要产生一个长的随机比特序列作为发送器的输入,发送器将这些比特调制成某种形式的信号以便传送到仿真信道,我们在传输信道上加上一定的可调制噪声,这些噪声信号会变成接收器的输入,接收器解调信号然后恢复比特序列,最后比较接收到的比特和传送的比特并计算错误。
误比特率性能常能描述成二维图像。纵坐标是归一化的信噪比,即每个比特的能量除以噪声的单边功率谱密度,单位为分贝。横坐标为误比特率,没有量纲。
2. 仿真过程及计算
1 运行发生器:通过发送器将伪随机序列变成数字化的调制信号。
2 设定信噪比:假定SNR为m dB,则Eb/N0=10,用MATLAB假设SNR单位为分贝。
3 确定Eb
4 计算N0
5 计算噪声的方差σn
6 产生噪声:因为噪声具有零均值,所以其功率和方差相等。我们产生一个和信号长度相同的噪声向量,且该向量方差为σn。
7 加上噪声,运行接收器
8 确定时间延迟
9 产生误差向量
10 统计错误比特:误差向量“err”中的每一个非零元素对应着一个错误的比特。
最后计算误比特率BER:每运行一次误比特率仿真,就需要传输和接收固定数量的比特,然后确定接收到的比特中有多少错误的。使用MATLAB计算BER: ber=te/length(tx)。
四.实验内容
1. 实验程序a
% Simulation of BPSK AWGN
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=qfunc(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')
2. 实验程序b
%Simulation of QPSK AWGN
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
noise=randn(1,N_number/2)+j.*randn(1,N_number/2);
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.*noise;
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=qfunc(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/NO(db)');
ylabel('BER');
figure
i=0:1:N_snr;
semilogy(i,SER,'-r',i,SER_T,'*b');
legend('SER-simulation','SER-theory');
xlabel('Eb/NO(db)');
ylabel('SER');
五.仿真结果
1.实验程序a

图a
3. 实验程序b

图b1

图b2
五.实验小结
通过本次实验,掌握了二相BPSK调制的工作原理及利用MATLAB进行误比特率测试BER的方法,学会了AWGN信道中BPSK调制系统的BER仿真计算方法。在实验过程中我通过不断的调试与学习,对本次实验的内容有了整体的把握,对MATLAB的使用也更加熟练,达到了预期的效果,收获很大。
实验课程名称:移动通信系统实验
实验项目名称
移动信道建模的仿真分析
实验成绩
实 验 者
专业班级
组    别
同 组 者
/
实验日期
一、实验目的
1. 无线通信信道的建模与仿真是实现移动通信系统仿真与分析的基础,宽带无线通与移动通信信道属频率选择性瑞利衰落信道模型。
2. 通过信道设计实验
1 掌握频率选择性信道模型的仿真建模方法
2 掌握模型中瑞利衰落系数的设计方法
3 掌握多径数目、功率和时延参数的设计
4 学会采用MATLAB语言对上述参数进行仿真。
二、实验仪器
1.计算器及操作系统
软件
三、实验方案和技术路线
1. 选择路径数
2. 按均匀分布产生各条路径的延迟
3. 按功率时延谱确定对应的各径的功率
4. 按Jake模型产生各径的瑞利衰落系数
5. 对瑞利衰落系数进行统计分析并与理论值相比较
说明:
1. 路径数目2-4自己确定,或采用某个国际标准
2. 每条路径时间延迟满足(0,Tmax)范围内均匀分布,Tmax为自己选择的最大采样步长数200-600间比较合适,或采用国际标准
3. 功率可以按时延迟谱求得,也可用国际标准测量值。功率延迟谱:①若采用等功率分配产生功率:Pi=Pt/M;②采用指数分布的功率延迟谱产生功率:P=1/6*exp(-t/6)
四、实验内容
实验程序如下:
% Simulation of Jakes Model
clear all;
f_max=30;
M=8; 
N=4*M+2;
Ts=1.024e-04;
sq=2/sqrt(N);
sigma=1/sqrt(2);
theta=0;
count=0;
t0=0.001;
for t=0:Ts:0.5
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 componentforoneoscillator
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 nagain
envelope=sqrt(tp.^2+tp1.^2);
rmsenv=sqrt(sum(envelope.^2)/count);
[auto_i,lag_i]=xcorr(tp,'coeff'); 
%Auto-correlation associated with inphase component
[auto_q,lag_q]=xcorr(tp,'coeff'); 
%Auto-correlation associated with quadrature component
len=length(lag_i);
[corrx2,lag2]=xcorr(tp,tp1,'coeff');
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;
sgma2=0.5;
x=[0:0.01:5];
pdf_theory=(x./sgma2).*exp(-1.*x.^2./(2.*sgma2));
figure;
plot(x,pdf_env,'-',x,pdf_theory,'*');
legend('Simulated','Theoretic');
xlabel('r');
ylabel('PDF of r');
五.仿真结果

图1

图2

图3

图4
六.实验小结
通过本次实验,我进行了无线通信信道的建模与仿真,认识到它是实现移动通信系统仿真与分析的基础,宽带无线通信与移动通信信道属频率选择性瑞利衰落信道模型。通过信道设计实验,基本掌握了频率选择性信道模型的仿真建模方法以及模型中瑞利衰落系数的设计方法,学会了多径数目、功率和时延参数的设计和采用MATLAB语言对上述参数进行仿真。实验过程中,我通过不断调试与修改,终于成功完成了任务,加深对MATLAB的应用的同时掌握了相关的知识。   
实验项目名称
CDMA通信系统仿真
实验成绩
实 验 者
专业班级
组    别
同 组 者
/
实验日期
一、实验目的
1. CDMA通信具有很多通信特点,不仅被IS-95移动通信系统使用,目前已成为3G的主要技术。
2. 通过实验:
(1)掌握直接序列扩频发射机与接收机的组成与仿真;
(2)仿真验证AWGN信道下单用户直接序列扩频系统的BER性能;
(3)仿真验证平坦瑞利信道下单用户直接序列扩频系统的BER性能;
(4)观察存在干扰用户时的系统性能变化。
二、实验仪器
1.计算器及操作系统
软件
三、实验原理
仿真基带直接序列扩频系统:
1. 采用BPSK或QPSK映射
2. 扩频序列可以是随机产生,可以是m序列,也可以是Gold码,长度自选
3. 最后对BER或SER随信噪比变化画图与理论单用户的结果比较,并对仿真结果进行分析.
四、实验方案与技术路线
1. 确定用户数目、信道特征以及调制方式
2. 确定基带扩频仿真系统的原理结构图,按照框图设计一个CDMA系统,并进行仿真。

3. 用MATLAB进行仿真,统计BER或SER随信噪比的关系,绘出曲线
4. 对统计试验的结果与单用户的理论值进行比较
5. 对仿真结果进行分析
五、仿真结果
1.仿真主程序
%main_IS95_forward.m
%此函数用于IS-95前向链路系统的仿真,包括扩
%频调制,匹配滤波,RAKE接收等相关通信模块。
%仿真环境: 加性高斯白噪声信道.
%数据速率 = 9600 KBps
%
clear all
close all
clc
disp('--------------start-------------------');
global ZiZqZs show R GiGq
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 TempResultsErrBErrC 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-------------------');
调用程序:
%VitEnc.m
function y = VitEnc(G, x);
% 此函数根据生成多项式进行Viterbi编码       
%+++++++++++++++++++++++variables++++++++++++++++++++++++++++
% G    生成多项式的矩阵
% x    输入数据(二进制形式)
% y    Viterbi编码输出序列
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
K = size(G, 1); %每个数据比特的码片数   
L = length(x); %输入数据的长度
yy = conv2(G, x'); %二维卷积
yy = yy(:, 1:L); %根据L重新设定yy长度
y = reshape(yy,K*L, 1);%矩阵变形
y = mod(y, 2); %模二运算
%SoftVitDec.m
function [xx, BestMetric] = SoftVitDec(G, y, ZeroTail);
%
% 此函数是实现软判决输入的Viterbi译码
%+++++++++++++++++++++++variables++++++++++++++++++++++++++++
% G            生成多项式的矩阵
% y            输入的待译码序列
% ZeroTail    判断是否包含‘0’尾
% xx          Viterbi译码输出序列
% BestMetric  最后的最佳度量
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++         
L = size(G, 1);    % 输出码片数
K= size(G, 2);      % 生成多项式的长度
N = 2^(K-1);        % 状态数
T = length(y)/L;    % 最大栅格深度
OutMtrx = zeros(N, 2*L); %输出矩阵的定义
for s = 1:N
in0 = ones(L, 1)*[0, (dec2bin((s-1), (K-1))-'0')];
in1 = ones(L, 1)*[1, (dec2bin((s-1), (K-1))-'0')];
out0 = mod(sum((G.*in0)'), 2);
out1 = mod(sum((G.*in1)'), 2);
OutMtrx(s, :) = [out0, out1]; %生成输出矩阵
end
OutMtrx = sign(OutMtrx-1/2);
PathMet = [100; zeros((N-1), 1)];      % 初始状态 = 100
PathMetTemp = PathMet(:,1); %副本
Trellis = zeros(N, T); %栅格的矩阵
Trellis(:,1) = [0 : (N-1)]';%给第一列赋值
y = reshape(y, L, length(y)/L);%矩阵按输出码片数变形
for t = 1:T %主栅格计算循环
yy = y(:, t); %取出y的第t列
for s = 0:N/2-1
[B0 ind0] = max(  PathMet(1+[2*s, 2*s+1]) + [OutMtrx(1+2*s, 0+[1:L]) * yy; OutMtrx(1+(2*s+1), 0+[1:L])*yy] );
[B1 ind1] = max(  PathMet(1+[2*s, 2*s+1]) + [OutMtrx(1+2*s, L+[1:L]) * yy; OutMtrx(1+(2*s+1), L+[1:L]) * yy] );
PathMetTemp(1+[s, s+N/2]) =  [B0; B1]; %改变状态
Trellis(1+[s, s+N/2], t+1) = [2*s+(ind0-1); 2*s + (ind1-1)];%生成栅格矩阵         
end
PathMet = PathMetTemp;%赋状态值
end
xx = zeros(T, 1);%生成单列0矩阵,输出变量
if (ZeroTail)  %确定最佳度量
BestInd = 1;
else
[Mycop, BestInd]  = max(PathMet); %非‘0’尾,取最大值所在位置
end
BestMetric = PathMet(BestInd); %得到最后的最佳度量
xx(T) = floor((BestInd-1)/(N/2)); %赋值xx最后一个数
NextState = Trellis(BestInd, (T+1)); %从栅格矩阵获得初态
for t=T:-1:2
xx(t-1) = floor(NextState/(N/2));%倒序生成xx
NextState = Trellis( (NextState+1), t); %从栅格矩阵获得次态
end
if (ZeroTail)
xx = xx(1:end-K+1);%限定译码输出序列长度
end
%ReceiverSD.m
function [DataOut, Metric] = ReceiverSD(SDchips, G, Scrambler);
% 此函数用于实现基于Viterbi译码的发送数据的恢复
%+++++++++++++++++++++++variables++++++++++++++++++++++++++++
% SDchips          软判决RAKE接收机输入符号
% G                Viterbi编码生成多项式矩阵
% Scrambler        扰码序列
% DataOut          接收数据(二进制形式)
% Metric            Viterbi译码最佳度量
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
if (nargin == 1)%判断只有SDchips传入时在此生成Viterbi编码生成多项式矩阵
G = [1 1 1 1 0 1 0 1 1; 1 0 1 1 1 0 0 0 1];
end
% 速率=19.2 KBps 
SDchips = SDchips.*sign(1/2-Scrambler);%解扰
INTERL = reshape(SDchips, 16, 24);%解交织       
SDchips = reshape(INTERL', length(SDchips), 1);  % 速率=19.2 KBps
[DataOut Metric] = SoftVitDec(G, SDchips, 1);%实现软判决输入的Viterbi译码
%ReceiverHD.m
function [DataOut, Metric] = ReceiverHD(HDchips, G, Scrambler);
% 此函数用于实现基于Viterbi译码的硬判决接收机
%+++++++++++++++++++++++variables++++++++++++++++++++++++++++
%HDchips          硬判决RAKE接收机输入符号
% G                Viterbi编码生成多项式矩阵
% Scrambler        扰码序列
% DataOut          接收数据(二进制形式)
% Metric            Viterbi译码最佳度量
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
if (nargin == 1) %判断只有HDchips传入时在此生成Viterbi编码生成多项式矩阵
G = [1 1 1 1 0 1 0 1 1; 1 0 1 1 1 0 0 0 1];
end
% 速率=19.2 KBps 
HDchips = xor(HDchips, Scrambler);%解扰
INTERL = reshape(HDchips, 16, 24);%解交织       
HDchips = reshape(INTERL', length(HDchips), 1);%速率=19.2 KBps 
[DataOut Metric] = VitDec(G, HDchips, 1);%维特比解码
%************************end of file***********************************
% ************************beginning of file*****************************
%VitDec.m
function [xx, BestMetric] = VitDec(G, y, ZeroTail);
%
% 此函数是实现硬判决输入的Viterbi译码
%+++++++++++++++++++++++variables++++++++++++++++++++++++++++
% G            生成多项式的矩阵
% y            输入的待译码序列
% ZeroTail    判断是否包含‘0’尾
% xx          Viterbi译码输出序列
% BestMetric  最后的最佳度量
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
L = size(G, 1);    % 输出码片数
K= size(G, 2);      % 生成多项式长度
N = 2^(K-1);        % 状态数
T = length(y)/L;    % 最大栅格深度
OutMtrx = zeros(N, 2*L);%输出矩阵的定义
for s = 1:N
in0 = ones(L, 1)*[0, (dec2bin((s-1), (K-1))-'0')];
in1 = ones(L, 1)*[1, (dec2bin((s-1), (K-1))-'0')];
out0 = mod(sum((G.*in0)'), 2);
out1 = mod(sum((G.*in1)'), 2);
OutMtrx(s, :) = [out0, out1];%生成输出矩阵
end
PathMet = [0; 100*ones((N-1), 1)];%初始状态为0
PathMetTemp = PathMet(:,1);%副本
Trellis = zeros(N, T);%栅格的矩阵
Trellis(:,1) = [0 : (N-1)]';%给第一列赋值
y = reshape(y, L, length(y)/L);%矩阵按输出码片数变形
for t = 1:T %主栅格计算循环
yy = y(:, t)';%取出y的第t列
for s = 0:N/2-1
[B0 ind0] = min(  PathMet(1+[2*s, 2*s+1]) + [sum(abs(OutMtrx(1+2*s, 0+[1:L]) - yy).^2); sum(abs(OutMtrx(1+(2*s+1), 0+[1:L]) - yy).^2)] );
[B1 ind1] = min(  PathMet(1+[2*s, 2*s+1]) + [sum(abs(OutMtrx(1+2*s, L+[1:L]) - yy).^2); sum(abs(OutMtrx(1+(2*s+1), L+[1:L]) - yy).^2)] );
PathMetTemp(1+[s, s+N/2]) =  [B0; B1];%改变状态
Trellis(1+[s, s+N/2], t+1) = [2*s+(ind0-1); 2*s + (ind1-1)];%生成栅格矩阵       
end
PathMet = PathMetTemp;%赋状态值 
end
xx = zeros(T, 1);%生成单列0矩阵,输出变量
if (ZeroTail) %确定最佳度量
BestInd = 1;
else
[Mycop, BestInd]  = min(PathMet);%非‘0’尾,取最小值所在位置
end
BestMetric = PathMet(BestInd);%得到最后的最佳度量
xx(T) = floor((BestInd-1)/(N/2));%赋值xx最后一个数
NextState = Trellis(BestInd, (T+1));%从栅格矩阵获得初态
for t=T:-1:2
xx(t-1) = floor(NextState/(N/2)); %倒序生成xx
NextState = Trellis( (NextState+1), t);%从栅格矩阵获得次态 
end
if (ZeroTail)
xx = xx(1:end-K+1); %限定译码输出序列长度
end
%PNGen.m
function [y, Z] = PNGen(G, Zin, N);
%
% 此函数是根据生成多项式和输入状态产生长度为N的伪随机序列
%+++++++++++++++++++++++variables++++++++++++++++++++++++++++
% G            生成多项式
% Zin          移位寄存器初始化
% N            PN序列长度
% y            生成的PN码序列
% Z            移位寄存器的输出状态
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++         
L = length(G);%扰码生成多项式长度
Z = Zin;    % 移位寄存器的初始化
y = zeros(N, 1);%N*1的0矩阵
for i=1:N
y(i) = Z(L); %获取当前状态输出值(移位寄存器的最后一位输出)
Z = xor(G*Z(L), Z); %生成移位寄存器次态
Z = [Z(L); Z(1:L-1)]; %移位寄存器后移1位
end   
%yy = filter(1, flipud(G), [1; zeros(N-1, 1)]);
%yy = mod(yy, 2);
%PacketBuilder.m
function [ChipsOut, Scrambler] = PacketBuilder(DataBits, G, Gs);
%此函数用于产生IS-95前向链路系统的发送数据包
%+++++++++++++++++++++++variables++++++++++++++++++++++++++++
% DataBits            发送数据(二进制形式)
% G                  Viterbi编码生成多项式
% Gs                  长序列生成多项式(扰码生成多项式)
% ChipsOut            输入到调制器的码序列(二进制形式)
% Scrambler          扰码
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++         
global Zs  %扰码状态
K = size(G, 2); %维特比多项式的长度   
L = size(G, 1);  %每个数据比特的码片数
N = 64*L*(length(DataBits)+K-1);% 码片数 (9.6 Kbps -> 1.288 Mbps)
chips = VitEnc(G, [DataBits; zeros(K-1,1)]);    % Viterbi编码
% 交织编码
INTERL = reshape(chips, 24, 16);            % IN:列, OUT:行
chips = reshape(INTERL', length(chips), 1);  %速率=19.2 KBps
% 产生扰码
[LongSeq Zs] = PNGen(Gs, Zs, N);%根据生成多项式和输入状态产生长度为N的PN序列
Scrambler = LongSeq(1:64:end);%扰码   
ChipsOut = xor(chips, Scrambler); %加扰
%************************end of file***********************************
%Modulator.m
function [TxOut, PN, MF] = Modulator(chips, MFType, Walsh);
%此函数用于实现IS-95前向链路系统的数据调制
%+++++++++++++++++++++++variables++++++++++++++++++++++++++++
% chips          发送的初始数据
% MFType          成型滤波器的类型选择
% Walsh          walsh码
% TxOut          调制输出信号序列
% PN              用于扩频调制的PN码序列
% MF              匹配滤波器参数
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++           
global Zi Zq show R Gi Gq
N = length(chips)*length(Walsh); %输出数据的数目 
% 输入速率 = 19.2 KBps, 输出速率= 1.2288 Mcps
tmp = sign(Walsh-1/2)*sign(chips'-1/2); %求扩频调制码数目的中间变量
chips = reshape(tmp, prod(size(tmp)), 1);%矩阵变形
[PNi Zi] = PNGen(Gi, Zi, N);%i路PN序列生成 
[PNq Zq] = PNGen(Gq, Zq, N);%q路PN序列生成 
PN = sign(PNi-1/2) + j*sign(PNq-1/2); %i、q路以复数形式合并
chips_out = chips.*PN;%得到复数形式的码序列 
chips = [chips_out, zeros(N, R-1)];%码序列0插值   
chips = reshape(chips.' , N*R, 1);%矩阵变形 
%成型滤波器
switch (MFType) %根据MFType选择滤波器类型
case 1
%升余弦滤波器 
L = 25;
L_2 = floor(L/2);
n = [-L_2:L_2]; %升余弦滤波器点数
B = 0.7; %B越大拖尾越小
MF = sinc(n/R).*(cos(pi*B*n/R)./(1-(2*B*n/R).^2)); %升余弦滤波器形状
MF = MF/sqrt(sum(MF.^2)); %升余弦滤波器特性曲线
case 2
%矩形滤波器 
L = R;
L_2 = floor(L/2);
MF = ones(L, 1); %1->0,锐截止
MF = MF/sqrt(sum(MF.^2)); %矩形滤波器特性曲线
case 3
%汉明滤波器 
L = R;
L_2 = floor(L/2);
MF = hamming(L);%生成汉明滤波器
MF = MF/sqrt(sum(MF.^2));%汉明滤波器特性曲线
end
MF = MF(:); %转置
TxOut = sqrt(R)*conv(MF, chips)/sqrt(2);%通过成型滤波器
TxOut = TxOut(L_2+1: end - L_2); %限定序列区间
if (show)
figure;
subplot(211); plot(MF, '-o'); title('Matched Filter'); grid on;%成型滤波器特性曲线图
subplot(212); psd(TxOut, 1024, 1e3, 113); title('Spectrum'); %功率谱密度估计
end
%Demodulator.m
function [SD] = Demodulator(RxIn, PN, MF, Walsh);
% 此函数是实现基于RAKE接收机的IS-95前向信链路系统的数据包的解调
%+++++++++++++++++++++++variables++++++++++++++++++++++++++++
% RxIn            输入信号
% PN              PN码序列(用于解扩)
% MF              匹配滤波器参数
% Walsh          用于解调的walsh码
% SD              RAKE接收机的软判决输出
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% DEMODULATOR            This function performs demodulation of the forward
%                                      Channel packet, based on RAKE Receiver
%                         Block Diagram
%                        Input Signal -> [Matched Filter] -> [Sampler] -> [RAKE Receiver] -> [Walsh] -> [DeSpreading]
%                         Inputs: RxIn - input signal (I/Q) analoge
%                                    PN - PN sequence (used for De-spreading)
%                                    MF - matched filter taps
%                                    Walsh - Used row of Walsh matrix for  recovering
%
%                         Outputs: SD - Soft Decisions of RAKE receiver
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++           
global R 
N = length(RxIn)/R; %有用码的个数
L = length(MF);
L_2 = floor(L/2);
rr = conv(flipud(conj(MF)), RxIn); %通过匹配接收滤波器
rr = rr(L_2+1: end - L_2); %限定接收符号序列长度
Rx = sign(real(rr(1:R:end))) + j*sign(imag(rr(1:R:end)));%接收符号采样
Rx = reshape(Rx, 64, N/64); %列导向   
Walsh = ones(N/64, 1)*sign(Walsh'-1/2);%行导向walsh码
PN = reshape(PN, 64, N/64)';%矩阵变形
PN = PN.*Walsh;%walsh正交               
% 输入速率 = 1.2288 Mpbs, 输出速率 = 19.2 KBps
SD= PN*Rx;%解扩     
SD= real(diag(SD));%确定软判决输出
2.仿真结果

图1

图2
六、实验小结
在本次实验中,我掌握了直接序列扩频发射机与接收机的组成与仿真,利用仿真验证了AWGN信道下单用户直接序列扩频系统的BER性能。我还通过仿真验证平坦瑞利信道下单用户直接序列扩频系统的BER性能并观察存在干扰用户时的系统性能变化。通过本次实验,我加深了对CDMA通信系统的认识,对其自身所具有的良好性能有了一定的掌握,通过实验提升了自己的综合技能,受益颇多。