1、% 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,
2、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:10switch seq case 1 % m 序列code = mseq(stage,pt
3、ap1,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); %+%+信道衰减初值设定
4、+ 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 = 10
5、00; % 仿真循环次数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
6、 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
7、 * 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 = qpskd
8、emod(ich9,qch9,user,nd,ml); % QPSK 解调%+%+误码率分析+ noe2 = sum(sum(abs(data-demodata); nod2 = user * nd * ml; noe = noe + noe2; nod = nod + nod2; % fprintf(%dt%en,ii,noe2/nod2); end %+%+数据文件+ ber = noe / nod; fprintf(%dt%dt%dt%en,ebn0,noe,nod,noe/nod); fid = fopen(BER.dat,a); fprintf(fid,%dt%et%ft%ftn,e
9、bn0,noe/nod,noe,nod); fclose(fid); err_rate_final(ebn0+6)=ber;end%+%+性能仿真图+ figuresemilogy(SNR,err_rate_final,b-*);xlabel(信噪比/dB)ylabel(误码率)axis(-5,10,0,1)grid ondisp(-end-); %+% *beginning of file*% hrollfcoef.m % % 此函数用于产生 Nyquist 滤波器的系数% function xh = hrollfcoef(irfn,ipoint,sr,alfs,ncc) %+variabl
10、es+ % 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 =
11、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
12、)+(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 %
13、 function mout = mseq(stg, taps, inidata, n) %+variables+ % stg m 序列阶数% taps 线性移位寄存器的系数% inidata 序列的初始化% n 输出序列的数目 % mout 输出的 m 序列%+ if nargin 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(in
14、regi,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 shift
15、r 0 outregi(1 :h-shiftu,:) = inregi(1+shiftu:h, :); outregi(h-shiftu+1:h, :) = inregi(1 :shiftu,:); elseif shiftu 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,:) = resh
16、ape(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*sam
17、ple); 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 输出序列虚部
18、% filter 滤波器的系数%+ iout = conv2(idata,filter); qout = conv2(qdata,filter); %*end of file*% *beginning of file*% sefade.m % % 此函数用于实现信道的频率选择性衰落设计% functioniout,qout,ramp,rcos,rsin=sefade(idata,qdata,itau,dlvl,th,n0,itn,n1,nsamp,tstp,fd,flat) %+variables+ % idata 输入的实部数据 % qdata 输入的虚部数据 % iout 输出的实部数据%
19、 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 .(
20、-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(
21、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
22、,nsamp);qout=zeros(1,nsamp);if idel = 0 iout(1:idel) = zeros(1,idel);qout(1:idel) = zeros(1,idel);endiout(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,ns
23、amp,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.*
24、(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 = x
25、s + 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:nsa
26、mp); 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 = l
27、ength(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
28、+ % 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,:) =
29、 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 输入数据的
30、实部 % 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
31、(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*
