% main_DSCDMA.m
% 这个仿真程序用于实现DS-CDMA通信系统仿真
%
%+++++++++++++++++++++++准备部分+++++++++++++++++++++++++++
clear all;
clc
sr = 256000.0; % 符号速率
ml = 2; % 调制阶数
br = sr * ml; % 比特速率
nd = 100; % 符号数
SNR=-5:1:10; % Eb/No
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%+++++++++++++++++++++滤波器初值设定+++++++++++++++++++++++
irfn = 21; % 滤波器阶数
IPOINT = 8; % 过采样倍数
alfs = 0.5; % 滚降因子
[xh] = hrollfcoef(irfn,IPOINT,sr,alfs,1);
[xh2] = hrollfcoef(irfn,IPOINT,sr,alfs,0);
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%++++++++++++++++++++++++扩频码初值设定+++++++++++++++++++
user = 1; % 用户数
seq = 1; % 1:m序列 2:Gold序列 3:正交Gold序列
stage = 3; % 序列阶数
ptap1 = [1 3]; % 第一个线性移位寄存器的系数
ptap2 = [2 3]; % 第二个线性移位寄存器的系数
regi1 = [1 1 1]; % 第一个线性移位寄存器的初始化
regi2 = [1 1 1]; % 第二个线性移位寄存器的初始化
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
disp('--------------start-------------------');
%+++++++++++++++++扩频码的产生+++++++++++++++++
for ebn0=-5:1:10
switch seq
case 1 % m序列
code = mseq(stage,ptap1,regi1,user);
case 2 % Gold序列
m1 = mseq(stage,ptap1,regi1);
m2 = mseq(stage,ptap2,regi2);
code = goldseq(m1,m2,user);
case 3 % 正交Gold序列
m1 = mseq(stage,ptap1,regi1);
m2 = mseq(stage,ptap2,regi2);
code = [goldseq(m1,m2,user),zeros(user,1)];
end
code = code * 2 - 1;
clen = length(code);
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%+++++++++++++++++++++信道衰减初值设定+++++++++++++++++++++++
rfade = 0; % 瑞利衰减 0:不考虑 1:考虑
itau = [0,8]; % 延时
dlvl1 = [0.0,40.0]; % 衰减电平
n0 = [6,7]; % 用于产生衰落的波数
th1 = [0.0,0.0]; % 延迟波形的初始相位
itnd1 = [3001,4004]; % 设置衰落计数
now1 = 2; % 主径和延迟波形总数
tstp = 1 / sr / IPOINT / clen; % 时间分辨率
fd = 160; % 多普勒频移[Hz]
flat = 1; % 平坦瑞利衰落环境
itndel = nd * IPOINT * clen * 30;
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%++++++++++++++++++++++仿真运算开始++++++++++++++++++++++++
nloop = 1000; % 仿真循环次数
noe = 0;
nod = 0;
for ii=1:nloop
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%++++++++++++++++++++++++发射机++++++++++++++++++++++++++++
data = rand(user,nd*ml) > 0.5;
[ich, qch] = qpskmod(data,user,nd,ml); % QPSK 调制
[ich1,qch1] = spread(ich,qch,code); % 扩频
[ich2,qch2] = compoversamp2(ich1,qch1,IPOINT); % 过采样
[ich3,qch3] = compconv2(ich2,qch2,xh); % 滤波
if user == 1
ich4 = ich3;
qch4 = qch3;
else
ich4 = sum(ich3);
qch4 = sum(qch3);
end
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%+++++++++++++++++++++++衰减信道++++++++++++++++++++++++++
if rfade == 0
ich5 = ich4;
qch5 = qch4;
else
[ich5,qch5] = sefade(ich4,qch4,itau,dlvl1,th1,n0,itnd1, ...
now1,length(ich4),tstp,fd,flat);
itnd1 = itnd1 + itndel;
end
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%+++++++++++++++++++++++++接收机++++++++++++++++++++++++++++
spow = sum(rot90(ich3.^2 + qch3.^2)) / nd; % 衰减计算
attn = sqrt(0.5 * spow * sr / br * 10^(-ebn0/10));
[ich6,qch6] = comb2(ich5,qch5,attn); % 添加高斯白噪声(AWGN)
[ich7,qch7] = compconv2(ich6,qch6,xh2); % 滤波
sampl = irfn * IPOINT + 1;
ich8 = ich7(:,sampl:IPOINT:IPOINT*nd*clen+sampl-1);
qch8 = qch7(:,sampl:IPOINT:IPOINT*nd*clen+sampl-1);
[ich9 qch9] = despread(ich8,qch8,code); % 解扩
demodata = qpskdemod(ich9,qch9,user,nd,ml); % QPSK解调
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%++++++++++++++++++++误码率分析++++++++++++++++++++
noe2 = sum(sum(abs(data-demodata)));
nod2 = user * nd * ml;
noe = noe + noe2;
nod = nod + nod2;
% fprintf('%d\t%e\n',ii,noe2/nod2);
end
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%+++++++++++++++++++++++++数据文件++++++++++++++++++++++++++
ber = noe / nod;
fprintf('%d\t%d\t%d\t%e\n',ebn0,noe,nod,noe/nod);
fid = fopen('BER.dat','a');
fprintf(fid,'%d\t%e\t%f\t%f\t\n',ebn0,noe/nod,noe,nod);
fclose(fid);
err_rate_final(ebn0+6)=ber;
end
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%+++++++++++++++++++++++++性能仿真图++++++++++++++++++++++++++
figure
semilogy(SNR,err_rate_final,'b-*');
xlabel('信噪比/dB')
ylabel('误码率')
axis([-5,10,0,1])
grid on
disp('--------------end-------------------');
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% ************************beginning of file*****************************
% hrollfcoef.m
%
% 此函数用于产生Nyquist滤波器的系数
%
function [xh] = hrollfcoef(irfn,ipoint,sr,alfs,ncc)
%+++++++++++++++++++++++variables++++++++++++++++++++++++++++
% irfn
用于滤波的符号数
% ipoint 每个符号的抽样值数
% sr 符号速率
% alfs 衰减截止频率
% ncc 1 -- 发送端滤波 0 -- 接收端滤波
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
xi=zeros(1,irfn*ipoint+1);
xq=zeros(1,irfn*ipoint+1);
point = ipoint;
tr = sr ;
tstp = 1.0 ./ tr ./ ipoint;
n = ipoint .* irfn;
mid = ( n ./ 2 ) + 1;
sub1 = 4.0 .* alfs .* tr;
for i = 1 : n
icon = i - mid;
ym = icon;
if icon == 0.0
xt = (1.0-alfs+4.0.*alfs./pi).* tr;
else
sub2 =16.0.*alfs.*alfs.*ym.*ym./ipoint./ipoint;
if sub2 ~= 1.0
x1=sin(pi*(1.0-alfs)/ipoint*ym)./pi./(1.0-sub2)./ym./tstp;
x2=cos(pi*(1.0+alfs)/ipoint*ym)./pi.*sub1./(1.0-sub2);
xt = x1 + x2;
else
xt = alfs.*tr.*((1.0-2.0/pi).*cos(pi/4.0/alfs)+(1.0+2.0./pi).*sin(pi/4.0/alfs))./sqrt(2.0);
end
end
if ncc == 0
%当接收机情况
xh( i ) = xt ./ ipoint ./ tr;
elseif ncc == 1 %当发射机情况
xh( i ) = xt ./ tr;
else
error('ncc error');
end
end
%************************end of file**********************************
% ************************beginning of file*****************************
% mseq.m
%
% 此函数产生m序列
%
% 试举一例:
% stg = 3
% taps = [ 1 , 3 ]
% inidata = [ 1 , 1 , 1 ]
% n = 2
%
function [mout] = mseq(stg, taps, inidata, n)
%+++++++++++++++++++++++variables++++++++++++++++++++++++++++
% stg m序列阶数
% taps 线性移位寄存器的系数
% inidata 序列的初始化
% n 输出序列的数目
% mout 输出的m序列
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
if nargin < 4
n = 1;
end
mout = zeros(n,2^stg-1);
fpos = zeros(stg,1);
fpos(taps) = 1;
for ii=1:2^stg-1
mout(1,ii) = inidata(stg); % 输出数据的存储
num = mod(inidata*fpos,2); % 反馈数据的计算
inidata(2:stg) = inidata(1:stg-1); % 线形移位寄存器的一次移位
inidata(1) = num; % 返回反馈值
end
if n > 1
for ii=2:n
mout(ii,:) = shift(mout(ii-1,:),1,0);
end
end
%************************end of file**********************************
% ************************beginning of file*****************************
% shift.m
%
% 此函数用于实现线性移位寄存器的移位操作
%
function [outregi] = shift(inregi,shiftr,shiftu)
%+++++++++++++++++++++++variables++++++++++++++++++++++++++++
% inrege 向量或矩阵
% shiftr 右移量
% shiftu 顶部移位量
% outregi 寄存器的输出
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[h, v] = size(inregi);
outregi = inregi;
shiftr = rem(shiftr,v);
shiftu = rem(shiftu,h);
if shiftr > 0
outregi(:,1 :shiftr) = inregi(:,v-shiftr+1:v );
outregi(:,1+shiftr:v ) = inregi(:,1 :v-shiftr);
elseif shiftr < 0
outregi(:,1 :v+shiftr) = inregi(:,1-shiftr:v );
outregi(:,v+shiftr+1:v ) = inregi(:,1 :-shiftr);
end
inregi = outregi;
if shiftu > 0
outregi(1 :h-shiftu,:) = inregi(1+shiftu:h, :);
outregi(h-shiftu+1:h, :) = inregi(1 :shiftu,:);
elseif shiftu < 0
outregi(1 :-shiftu,:) = inregi(h+shiftu+1:h, :);
outregi(1-shiftu:h, :) = inregi(1 :h+shiftu,:);
end
%************************end of file**********************************
% ************************beginning of file*****************************
% goldseq.m
%
%此函数用于产生gold序列
%
function [gout] = goldseq(m1, m2, n)
%+++++++++++++++++++++++variables++++++++++++++++++++++++++++
% m1 第一个m序列
% m2 第二个m序列
% n 输出序列的数目
% gout 输出产生的gold序列
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
if nargin < 3
n = 1;
end
gout = zeros(n,length(m1));
for ii=1:n
gout(ii,:) = xor(m1,m2);
m2 = shift(m2,1,0);
end
%************************end of file**********************************
% ************************beginning of file*****************************
% qpskmod.m
%
% 此函数实现QPSK调制
%
function [iout,qout]=qpskmod(paradata,para,nd,ml)
%+++++++++++++++++++++++variables++++++++++++++++++++++++++++
% paradata 输入数据
% iout 输出的实部数据
% qout 输出的虚部数据
% para 并行信道数
% nd 输入数据个数
% ml 调制阶数
% %++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
m2=ml./2;
paradata2=paradata.*2-1;
count2=0;
for jj=1:nd
isi = zeros(para,1);
isq = zeros(para,1);
for ii = 1 : m2
isi = isi + 2.^( m2 - ii ) .* paradata2((1:para),ii+count2);
isq = isq + 2.^( m2 - ii ) .* paradata2((1:para),m2+ii+count2);
end
iout((1:para),jj)=isi;
qout((1:para),jj)=isq;
count2=count2+ml;
end
%************************end of file**********************************
% ************************beginning of file*****************************
% spread.m
%
%此函数用于实现扩频调制
%
function [iout, qout] = spread(idata, qdata, code1)
%+++++++++++++++++++++++variables++++++++++++++++++++++++++++
% idata 输入序列实部
% qdata 输入序列虚部
% iout 输出序列实部
% qout 输出序列虚部
% code1 扩频码序列
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
switch nargin
case { 0 , 1 }
error('lack of input argument');
case 2
code1 = qdata;
qdata = idata;
end
[hn,vn] = size(idata);
[hc,vc] = size(code1);
if hn > hc
error('lack of spread code sequences');
end
iout = zeros(hn,vn*vc);
qout = zeros(hn,vn*vc);
for ii=1:hn
iout(ii,:) = reshape(rot90(code1(ii,:),3)*idata(ii,:),1,vn*vc);
qout(ii,:) = reshape(rot90(code1(ii,:),3)*qdata(ii,:),1,vn*vc);
end
%************************end of file**********************************
% ************************beginning of file*****************************
% compoversamp2.m
%
% 此函数实现"sample"倍升采样
%
function [iout,qout] = compoversamp2(iin, qin, sample)
%+++++++++++++++++++++++variables++++++++++++++++++++++++++++
% iin 输入序列实部
% qin 输入序列虚部
% iout 输出序列实部
% qout 输出序列虚部
% sample 升采样的倍数
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[h,v] = size(iin);
iout = zeros(h,v*sample);
qout = zeros(h,v*sample);
iout(:,1:sample:1+sample*(v-1)) = iin;
qout(:,1:sample:1+sample*(v-1)) = qin;
%************************end of file**********************************
% ************************beginning of file*****************************
% compconv2.m
%
% 此函数用于实现有用信号的滤波
%
function [iout, qout] = compconv2(idata, qdata, filter)
%+++++++++++++++++++++++variables++++++++++++++++++++++++++++
% idata 输入序列实部
% qdata 输入序列虚部
% iout 输出序列实部
% qout 输出序列虚部
% filter 滤波器的系数
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
iout = conv2(idata,filter);
qout = conv2(qdata,filter);
%************************end of file**********************************
% ************************beginning of file*****************************
% sefade.m
%
% 此函数用于实现信道的频率选择性衰落设计
%
function[iout,qout,ramp,rcos,rsin]=sefade(idata,qdata,itau,dlvl,th,n0,itn,n1,nsamp,tstp,fd,flat)
%+++++++++++++++++++++++variables++++++++++++++++++++++++++++
% idata 输入的实部数据
% qdata 输入的虚部数据
% iout 输出的实部数据
% qout 输出的虚部数据
% ramp 幅度衰减
% rcos 正交分量衰减
% rsin 同相分量衰减
% itau 各径时延
% dlvl 各多径衰减量
% th 各多径初始相位
% n0 用于产生各径衰落的波数
% itn 各多径衰减计数
% n1 主径和各延迟波形总数
% nsamp 符号数
% tstp 最小时间分辨率
% fd 最大多普勒频移
% flat 是否是平坦衰落
% (1->flat (只有幅度衰落),0->nomal(相位和幅度都衰落)
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
iout = zeros(1,nsamp);
qout = zeros(1,nsamp);
total_attn = sum(10 .^( -1.0 .* dlvl ./ 10.0));
for k = 1 : n1
atts = 10.^( -0.05 .* dlvl(k));
if dlvl(k) >= 40.0
atts = 0.0;
end
theta = th(k) .* pi ./ 180.0;
[itmp,qtmp] = delay ( idata , qdata , nsamp , itau(k));
[itmp3,qtmp3,ramp,rcos,rsin] = fade (itmp,qtmp,nsamp,tstp,fd,n0(k),itn(k),flat);
iout = iout + atts .* itmp3 ./ sqrt(total_attn);
qout = qout + atts .* qtmp3 ./ sqrt(total_attn);
end
%************************end of file**********************************
% ************************beginning of file*****************************
% delay.m
% 此函数用以实现信号的延迟传输
%
function [iout,qout] = delay( idata, qdata , nsamp , idel )
%+++++++++++++++++++++++variables++++++++++++++++++++++++++++
% idata 输入序列实部
% qdata 输入序列虚部
% iout 输出序列实部
% qout 输出序列虚部
% nsamp 仿真的抽样值数量
% idel 延迟的抽样值数量
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
iout=zeros(1,nsamp);
qout=zeros(1,nsamp);
if idel ~= 0
iout(1:idel) = zeros(1,idel);
qout(1:idel) = zeros(1,idel);
end
iout(idel+1:nsamp) = idata(1:nsamp-idel);
qout(idel+1:nsamp) = qdata(1:nsamp-idel);
%************************end of file**********************************
% ************************beginning of file*****************************
% fade.m
%
% 此函数实现信道的瑞利衰减设计
%
function [iout,qout,ramp,rcos,rsin]=fade(idata,qdata,nsamp,tstp,fd,no,counter,flat)
%+++++++++++++++++++++++variables++++++++++++++++++++++++++++
% idata 输入序列实部
% qdata 输入序列虚部
% iout 输出序列实部
% qout 输出序列虚部
% ramp 幅度衰减
% rcos 正交分量衰减
% rsin 同相分量衰减
% nsamp 符号数
% tstp 最小时间分辨率
% fd 最大多普勒频移
% no 用于产生衰落的波数
% counter 衰落计数
% flat 是否是平坦衰落
% (1->flat (只有幅度衰落),0->nomal(相位和幅度都衰落)
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
if fd ~= 0.0
ac0 = sqrt(1.0 ./ (2.0.*(no + 1)));
as0 = sqrt(1.0 ./ (2.0.*no));
ic0 = counter;
pai = 3.14159265;
wm = 2.0.*pai.*fd;
n = 4.*no + 2;
ts = tstp;
wmts = wm.*ts;
paino = pai./no;
xc=zeros(1,nsamp);
xs=zeros(1,nsamp);
ic=[1:nsamp]+ic0;
for nn = 1: no
cwn = cos( cos(2.0.*pai.*nn./n).*ic.*wmts );
xc = xc + cos(paino.*nn).*cwn;
xs = xs + sin(paino.*nn).*cwn;
end
cwmt = sqrt(2.0).*cos(ic.*wmts);
xc = (2.0.*xc + cwmt).*ac0;
xs = 2.0.*xs.*as0;
ramp=sqrt(xc.^2+xs.^2);
rcos=xc./ramp;
rsin=xs./ramp;
if flat ==1
iout = sqrt(xc.^2+xs.^2).*idata(1:nsamp);
qout = sqrt(xc.^2+xs.^2).*qdata(1:nsamp);
else
iout = xc.*idata(1:nsamp) - xs.*qdata(1:nsamp);
qout = xs.*idata(1:nsamp) + xc.*qdata(1:nsamp);
end
else
iout=idata;
qout=qdata;
end
%************************end of file**********************************
% ************************beginning of file*****************************
% comb2.m
%
% 此函数实现信道的高斯白噪声
%
function [iout, qout] = comb2(idata, qdata, attn)
%%+++++++++++++++++++++++variables++++++++++++++++++++++++++++
% idata 输入序列实部
% qdata 输入序列虚部
% iout 输出序列实部
% qout 输出序列虚部
% attn 根据信噪比得到的信号衰减水平
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
v = length(idata);
h = length(attn);
iout = zeros(h,v);
qout = zeros(h,v);
for ii=1:h
iout(ii,:) = idata + randn(1,v) * attn(ii);
qout(ii,:) = qdata + randn(1,v) * attn(ii);
end
%************************end of file**********************************
% ************************beginning of file*****************************
% despread.m
%
% 此函数实现数据的解频
%
function [iout, qout] = despread(idata, qdata, code1)
%+++++++++++++++++++++++variables++++++++++++++++++++++++++++
% idata 输入序列实部
% qdata 输入序列虚部
% iout 输出序列实部
% qout 输出序列虚部
% code1 扩频码序列
% %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
switch nargin
case { 0 , 1 }
error('lack of input argument');
case 2
code1 = qdata;
qdata = idata;
end
[hn,vn] = size(idata);
[hc,vc] = size(code1);
vn = fix(vn/vc);
iout = zeros(hc,vn);
qout = zeros(hc,vn);
for ii=1:hc
iout(ii,:) = rot90(flipud(rot90(reshape(idata(ii,:),vc,vn)))*rot90(code1(ii,:),3));
qout(ii,:) = rot90(flipud(rot90(reshape(qdata(ii,:),vc,vn)))*rot90(code1(ii,:),3));
end
%************************end of file**********************************
% ************************beginning of file*****************************
% qpskdemod.m
%
% 此函数实现QPSK解调
%
function [demodata]=qpskdemod(idata,qdata,para,nd,ml)
%+++++++++++++++++++++++variables++++++++++++++++++++++++++++
% idata 输入数据的实部
% qdata 输入数据的虚部
% demodata 解调后的数据
% para 并行的信道数
% nd 输入数据个数
% ml 调制阶数
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
demodata=zeros(para,ml*nd);
demodata((1:para),(1:ml:ml*nd-1))=idata((1:para),(1:nd))>=0;
demodata((1:para),(2:ml:ml*nd))=qdata((1:para),(1:nd))>=0;
%************************end of file**********************************
% ************************beginning of file*****************************
% autocorr.m
%
% 此函数实现一个序列的自相关运算
%
function [out] = autocorr(indata, tn)
%%+++++++++++++++++++++++variables++++++++++++++++++++++++++++
% indata 输入序列
% tn 序列的周期长度
% out 自相关函数
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
if nargin < 2
tn = 1;
end
ln = length(indata);
out = zeros(1,ln*tn);
for ii=0:ln*tn-1
out(ii+1) = sum(indata.*shift(indata,ii,0));
end
%************************end of file**********************************
% ************************beginning of file*****************************
% crosscorr.m
%
% 此函数实现两个序列的互相关运算
%
function [out] = crosscorr(indata1, indata2, tn)
%+++++++++++++++++++++++variables++++++++++++++++++++++++++++
% indata1 第一个输入序列
% indata2 第二个输入序列
% tn 序列的周期长度
% out 互相关运算的结果
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
if nargin < 3
tn = 1;
end
ln = length(indata1);
out = zeros(1,ln*tn);
for ii=0:ln*tn-1
out(ii+1) = sum(indata1.*shift(indata2,ii,0));
end
%************************end of file**********************************
本文档为【DS_CDMA 仿真源程序】,请使用软件OFFICE或WPS软件打开。作品中的文字与图均可以修改和编辑,
图片更改请在作品中右键图片并更换,文字修改请直接点击文字进行修改,也可以新增和删除文档中的内容。
该文档来自用户分享,如有侵权行为请发邮件ishare@vip.sina.com联系网站客服,我们会及时删除。
[版权声明] 本站所有资料为用户分享产生,若发现您的权利被侵害,请联系客服邮件isharekefu@iask.cn,我们尽快处理。
本作品所展示的图片、画像、字体、音乐的版权可能需版权方额外授权,请谨慎使用。
网站提供的党政主题相关内容(国旗、国徽、党徽..)目的在于配合国家政策宣传,仅限个人学习分享使用,禁止用于任何广告和商用目的。