DS_CDMA 仿真源程序

% 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**********************************