- 1 -
目录
一、背景 .................................................................................................... 4 二、基本要求 ............................................................................................ 4 三、设计概述 .......................................................................................... 4 四、Matlab设计流程图 ........................................................................... 5 五、Matlab程序及仿真结果图 ............................................................... 6
1、生成m序列及m序列性质 ............................................................................ 6 2、生成50位随机待发送二进制比特序列，并进行扩频编码......................... 7 3、对扩频前后信号进行BPSK调制，观察其时域波形................................... 9 4、计算并观察扩频前后BPSK调制信号的频谱............................................. 10 5、仿真经awgn信道传输后，扩频前后信号时域及频域的变化 .................. 11 6、对比经信道前后两种信号的频谱变化......................................................... 12 7、接收机与本地恢复载波相乘，观察仿真时域波形..................................... 14 8、与恢复载波相乘后，观察其频谱变化......................................................... 15 9、仿真观察信号经凯萨尔窗低通滤波后的频谱............................................. 16 10、观察经过低通滤波器后无扩频与扩频系统的时域波形........................... 17 11、对扩频系统进行解扩，观察其时域频域 ................................................... 18 12、比较扩频系统解扩前后信号带宽............................................................... 19 13、比较解扩前后信号功率谱密度................................................................... 20 14、对解扩信号进行采样、判决....................................................................... 21 15、在信道中加入2040~2050Hz窄带强干扰并乘以恢复载波 ...................... 24 16、对加窄带干扰的信号进行低通滤波并解扩............................................... 25 17、比较解扩后信号与窄带强干扰的功率谱................................................... 27
六、误码率simulink仿真 ..................................................................... 28
1、直接扩频系统信道模型................................................................................. 28
2、加窄带干扰的直扩系统建模......................................................................... 29 3、用示波器观察发送码字及解扩后码字......................................................... 30 4、直接扩频系统与无扩频系统的误码率比较................................................. 31 5、不同扩频序列长度下的误码率比较............................................................. 32 6、扩频序列长度N=7时，不同强度窄带干扰下的误码率比较 ................... 33
七、利用Walsh码实现码分多址技术 ................................................. 34
1、产生改善的walsh码 ..................................................................................... 35
2、产生两路不同的信息序列............................................................................. 36 3、用两个沃尔什码分别调制两路信号............................................................. 38
- 2 -
4、两路信号相加，并进行BPSK调制............................................................. 39 5、观察调制信号频谱，并经awgn信道加高斯白噪和窄带强干扰 .............. 40 6、接收机信号乘以恢复载波，观察时域和频域............................................. 42 7、信号经凯萨尔窗低通滤波器......................................................................... 43 8、对滤波后信号分别用m1和m2进行解扩 .................................................. 44 9、对两路信号分别采样，判决......................................................................... 45
八、产生随机序列Gold码和正交Gold码 ......................................... 47
1、产生Gold码并仿真其自相关函数............................................................... 48
2、产生正交Gold码并仿真其互相关函数....................................................... 50
九、实验心得体会 .................................................................................. 51
- 3 -
直接序列扩频系统仿真
一、背景
直接序列扩频通信系统(DSSS)是目前应用最为广泛的系统。在发送端，直扩系统将发送序列用伪随机序列扩展到一个很宽的频带上去，在接受端又用相同的扩频序列进行解扩，回复出原有信息。由于干扰信息与伪随机序列不相关，扩频后能够使窄带干扰得到有效的抑制，提高输出信噪比。系统框图如下图所示：
二、基本要求：
1. 通过matlab建模，对直扩系统进行仿真，数据调制方式可以自由选择，可以
使用基带信号，但最好能使用频带信号，信道为高斯白噪信道。要仿真出扩频前的信号的频偏，扩频后的信号频谱，过信道之后的频谱以及解扩之后的频谱。
2. 研究并仿真产生m序列，写出生成m序列的算法。
3. 验证直扩系统对窄带干扰的抑制能力，在信道中加入一个窄带强干扰，仿真
出加了干扰后的频谱图和解扩后的频谱图，给出误码率等仿真图。
4. 在以上基础上仿真实现码分多址技术，使用Walsh码进行复用，实现多个信
号同时传输。(选做) 可选项：
1.在信道中加入多径，使用rake接收来抗多径效应。
2.产生除m序列之外的其他随机序列，如Gold码，正交Gold码等等。 3.对比无扩频的系统的误码率。
三、设计概述
本次课设完成基本要求，并选作了可选项码分多址，Gold码及误码率对比。通过matlab建模仿真了直扩系统BPSK调制的各点频偏及时域信号，并仿真了窄带强干扰对直扩系统的影响以及利用改善的WALSH码实现码分多址技术。另外，通过matlab的simulink工具盒bertool工具仿真对比了直扩系统和无扩频系统的误码率。
- 4 -
四、matlab设计流程图
基本扩频系统仿真流程图
100/7Hz 二进制比特信息 100Hz 7位双极性m序列 100Hz 扩频序列 2000Hz 载波cos4000πt BPSK调制信号 高斯白噪声 恢复载波cos4000πt 凯萨尔滤波器低通滤波 100Hz 7位双极性m序列 采样、判决 - 5 -
五、matlab程序及仿真结果图
1、生成m序列及m序列性质
实验产生7位m序列，频率100Hz，模拟线性反馈移位寄存器序列，原理图如下：
clear all; clc;
X1=0;X2=0;X3=1;
m=350; %重复50遍的7位单极性m序列 for i=1:m
Y3=X3; Y2=X2; Y1=X1; X3=Y2; X2=Y1; X1=xor(Y3,Y1); L(i)=Y1; end
for i=1:m
M(i)=1-2*L(i); %将单极性m序列变为双极性m序列 end
k=1:1:m; figure(1)
subplot(3,1,1) %做m序列图 stem(k-1,M); axis([0,7,-1,1]); xlabel('k');
ylabel('M序列');
title('移位寄存器产生的双极性7位M序列') ; subplot(3,1,2) ym=fft(M,4096); magm=abs(ym); %求双极性m序列频谱 fm=(1:2048)*200/2048;
plot(fm,magm(1:2048)*2/4096);
- 6 -
title('双极性7位M序列的频谱') axis([90,140,0,0.1]);
[a,b]=xcorr(M,'unbiased'); subplot(3,1,3) %求双极性m序列自相关函数 plot(b,a);
axis([-20,20,-0.5,1.2]);
title('双极性7位M序列的自相关函数');
由上图可以看出，7位m序列为1，-1，-1，-1，1，-1，1。另外，自相关函数的图形比较尖锐，最大值为1，最小值为-1/7，符合理论结果。
2、生成50位随机待发送二进制比特序列，并进行扩频编码
生成的信息码频率为100/7Hz，利用m序列编码后，频率变为100Hz。
N=50;a=0;
x_rand=rand(1,N); %产生50个0与1之间随机数 for i=1:N
if x_rand(i)>=0.5 %大于等于0.5的取1，小于0.5的取0 x(i)=1;a=a+1; else x(i)=0; end end
- 7 -
t=0:N-1; figure(2) %做信息码图 subplot(2,1,1) stem(t,x);
title('扩频前待发送二进制信息序列'); tt=0:349; subplot(2,1,2) l=1:7*N; y(l)=0; for i=1:N k=7*i-6;
y(k)=x(i);
k=k+1;y(k)=x(i);k=k+1;y(k)=x(i);k=k+1;y(k)=x(i);k=k+1;y(k)=x(i);k=k+1;y(k)=x(i);k=k+1;y(k)=x(i); end s(l)=0;
for i=1:350 %扩频后，码率变为100/7*7=100Hz s(i)=xor(L(i),y(i)); end
tt=0:7*N-1; stem(tt,s);
axis([0,350,0,1]);
title('扩频后的待发送序列码');
- 8 -
3、对扩频前后信号进行BPSK调制，观察其时域波形
BPSK调制采用2kHz信号cos(2*2000*t)作为载波
figure(3)
subplot(2,1,2) fs=2000;
ts=0:0.00001:3.5-0.00001;%为了使信号看起来更光滑，作图时采样频率为100kHz % ps=cos(2*pi*fs*ts);
s_b=rectpulse(s,1000); %将冲激信号补成矩形信号
s_bpsk=(1-2.*s_b).*cos(2*pi*fs*ts);%扩频后信号BPSK调制时域波形，(1-2.*s_b)
是1，-1序列
plot(ts,s_bpsk); xlabel('s');
axis([0.055,0.085,-1.2,1.2])
title('扩频后bpsk信号时域波形'); subplot(2,1,1)
s_bb=rectpulse(x,7000);
s_bpskb=(1-2.*s_bb).*cos(2*pi*fs*ts);%无扩频信号BPSK调制时域波形 plot(ts,s_bpskb); xlabel('s');
axis([0.055,0.085,-1.2,1.2]);
title('扩频前bpsk信号时域波形')
可以看出，100/7Hz的无扩频信号每0.07s时由于序列极性变换产生相位变
- 9 -
换，100Hz的扩频后调制信号每0.01s由于序列极性变换产生相位变换。
- 10 -
4、计算并观察扩频前后BPSK调制信号的频谱
对信号采用400000点fft计算，得到频谱 figure(4) N=400000;
ybb=fft(s_bpskb,N); %无扩频信号BPSK调制频谱 magb=abs(ybb);
fbb=(1:N/2)*100000/N; subplot(2,1,1)
plot(fbb,magb(1:N/2)*2/N); axis([1700,2300,0,0.8]);
title('扩频前调制信号频谱'); xlabel('Hz'); subplot(2,1,2)
yb=fft(s_bpsk,N); %扩频信号BPSK调制频谱 mag=abs(yb);
fb=(1:N/2)*100000/N; plot(fb,mag(1:N/2)*2/N); axis([1700,2300,0,0.8]);
title('扩频后调制信号频谱'); xlabel('Hz');
如图，扩频前信号主瓣宽度约为2*100/7=28Hz，扩频后，信号频谱展宽，主瓣1900~2100Hz约为200Hz，为无扩频信号频谱宽度的N=7倍，符合理论推算。
- 11 -
5、仿真经awgn信道传输后，扩频前后信号时域及频域的变化
awgn信道模拟了真实的信道，为传输信号增加了高斯白噪声。在本次仿真中，设定信道信噪比为3dB，即信噪比约为2。
figure(5)
subplot(2,2,1)
s_bpskba=awgn(s_bpskb,3,'measured');%经过信道加高斯白噪,信噪比为3dbw plot(ts,s_bpskb,ts,s_bpskba); axis([0,0.005,-1.2,1.2]); xlabel('t');
title('经过信道加噪后的信号与原信号时域波形对比');
subplot(2,2,3)
s_bpska=awgn(s_bpsk,3,'measured'); plot(ts,s_bpsk,ts,s_bpska);
title('扩频后经加噪过信道后的信号与原信号时域波形对比'); xlabel('t');
axis([0.0675,0.0725,-1.2,1.2]);
subplot(2,2,2)
ybba=fft(s_bpskba,N); %无扩频调制信号经信道后频谱分析 magba=abs(ybba);
plot(fbb,magba(1:N/2)*2/N);
title('扩频前经信道调制信号频谱'); axis([1700,2300,0,0.8]); xlabel('Hz');
subplot(2,2,4)
yba=fft(s_bpska,N); %扩频调制信号经信道后频谱分析 maga=abs(yba);
fb=(1:N/2)*100000/N; plot(fb,maga(1:N/2)*2/N); axis([1700,2300,0,0.8]); xlabel('Hz');
title('扩频后经信道调制信号频谱');
- 12 -
可以看出，蓝色为光滑余弦调制信号，绿色为加噪声后时域波形，出现较大锯齿。至于频谱变化，这张图并不明显，于是我在下一张图继续比较了加入高斯噪声后的频谱变化。
6、对比经信道前后两种信号的频谱变化
figure(6)
title('对比经信道前后的信号频谱'); subplot(2,2,1)
plot(fbb,magb(1:N/2)*2/N); axis([0,4000,0,0.04]);
title('扩频前调制信号频谱'); xlabel('Hz'); subplot(2,2,2)
plot(fbb,magba(1:N/2)*2/N); axis([0,4000,0,0.04]);
title('扩频前经信道调制信号频谱'); xlabel('Hz'); subplot(2,2,3)
plot(fb,mag(1:N/2)*2/N);
- 13 -
axis([0,4000,0,0.04]);
title('扩频后调制信号频谱'); xlabel('Hz'); subplot(2,2,4)
plot(fb,maga(1:N/2)*2/N); axis([0,4000,0,0.04]);
title('扩频后经信道调制信号频谱'); xlabel('Hz');
由上图可以清楚地对比经高斯白噪声信道前后的频谱对比。虽然整体的幅度趋势不变，但是能看出，经过信道加噪后，在所有的频率点上都产生了一定的振幅，符合高斯白噪声的原理。此处的信噪比为3dB。另外可以看出，BPSK调制将信号频谱搬移到了以2000Hz为中心频率的频段上。
下面，我将分别仿真解调解扩后译码输出及加入窄带强干扰后解调解扩译码输出的时域和频谱。
- 14 -
（1）不加窄带强干扰
7、接收机与本地恢复载波相乘，观察仿真时域波形
figure(7)
subplot(2,1,1)
reb=s_bpskba.*cos(2*pi*fs*ts); %无扩频系统接收信号乘以本地恢复载波信号 plot(ts,reb);
axis([0.055,0.085,-1.5,1.5]); xlabel('t');
title('扩频前接收信号乘以恢复载波'); subplot(2,1,2)
re=s_bpska.*cos(2*pi*fs*ts); %扩频系统接收信号乘以本地恢复载波信号 plot(ts,re);
axis([0.055,0.085,-1.5,1.5]); xlabel('t');
title('扩频后接收信号乘以恢复载波');
可以看出，接收信号乘以恢复载波后，已经能大致恢复出信号的变化。同时，无扩频系统符号速率仍然是100/7Hz即0.07s出现符号变化，扩频系统100Hz即0.01s出现符号变化。
- 15 -
8、与恢复载波相乘后，观察其频谱变化
figure(8)
subplot(2,1,1) yreb=fft(reb,N); magreb=abs(yreb);
freb=(1:N/2)*100000/N;
plot(freb,magreb(1:N/2)*2/N); axis([0,5000,0,0.5]);
title('扩频前乘以恢复载波后信号频谱'); subplot(2,1,2) yre=fft(re,N); magre=abs(yre);
plot(freb,magre(1:N/2)*2/N);
title('扩频后乘以恢复载波后信号频谱'); axis([0,5000,0,0.5]);
可以看出，信号乘以频率为2kHz的恢复载波后，在基带和4kHz处存在频谱分量，则下一步需要对信号进行低通滤波。
- 16 -
9、仿真观察信号经凯萨尔窗低通滤波后的频谱
figure(9)
subplot(2,1,1) fp=100; fc=200;
as=100;ap=1; %衰减100dB fsw=22000; wp=2*fp/fsw; wc=2*fc/fsw;
Nw=ceil((as-7.95)/(14.36*(wc-wp)/2))+1; %求凯萨尔窗低通滤波器阶数 beta=0.1102*(as-8.7);
window=kaiser(Nw+1,beta);
b=fir1(Nw,wc,window);%返回截止频率为wc的Nw阶的低通滤波器系数向量 bs=abs(freqz(b,1,400000,fsw))';%频率响应，400000点的fft变换 plot(bs)
magrebl=bs.*magreb; %频谱与低通滤波器相乘 plot(freb,magrebl(1:N/2)*2/N); axis([0,200/7,0,1]);
title('扩频前:信号经过凯萨尔窗函数低通滤波'); xlabel('Hz'); subplot(2,1,2)
magrel=bs.*magre;
plot(freb,magrel(1:N/2)*2/N);
title('扩频后:信号经过凯萨尔窗函数低通滤波'); axis([0,200,0,0.4]); xlabel('Hz');
- 17 -
10、观察经过低通滤波器后无扩频与扩频系统的时域波形
figure(10) subplot(2,1,1)
yrebl=real(ifft(bs.*yreb,400000)); %对无扩频系统频谱做ifft变换 tm=(1:N)/N*4; plot(tm,yrebl); xlabel('t');
title('扩频前经过凯萨尔窗函数滤波后时域波形'); subplot(2,1,2)
yrel=real(ifft(bs.*yre,400000)); %对扩频系统频谱做ifft变换 plot(tm,yrel); xlabel('t');
title('扩频后经过凯萨尔窗函数滤波后时域波形');
如图，经过低通滤波器后，高频分量基本消失，剩下的信号已经能够进行采样判决，时域波形与原信息基本吻合。
- 18 -
11、对扩频系统进行解扩，观察其时域频域
figure(11) subplot(2,1,1)
jj=rectpulse(M,1000); %扩频信号乘以解扩序列 yrej=jj.*yrel(1:350000); plot(ts(1:350000),yrej); xlabel('t');
axis([0,4,-0.5,0.5]);
title('解扩后信号波形'); subplot(2,1,2) yj=fft(yrej,N); magj=abs(yj);
plot(freb,magj(1:N/2)*2/N); axis([0,500,0,0.2]);
title('解扩后信号频谱'); xlabel('Hz');
由于扩频信号与m序列具有良好的相关性，故乘以m序列以后，能基本还原出原信号波形。同时可以看出，频谱已经由扩展带宽再次缩短，还原出原信号频谱。
- 19 -
12、比较扩频系统解扩前后信号带宽
figure(12)
title('解扩前后信号频偏对比'); subplot(2,1,1)
plot(freb,magrel(1:N/2)*2/N); axis([0,200,0,0.4]);
title('解扩前信号频偏'); subplot(2,1,2)
plot(freb,magj(1:N/2)*2/N); axis([0,200,0,0.4]);
title('解扩后信号频偏');
可以清楚看出，解扩前信号主瓣约为100Hz，解扩后恢复为100/7Hz，与发送信息吻合。
- 20 -
13、比较解扩前后信号功率谱密度
figure(13) subplot(2,1,1) yjb=fft(yrel,N);
prelb=yjb.*conj(yjb)/N; plot(freb,prelb(1:N/2)*2/N); axis([0,200,0,0.01]);
title('解扩前信号功率谱');xlabel('Hz'); subplot(2,1,2) yj=fft(yrej,N);
prel=yj.*conj(yj)/N;
plot(freb,prel(1:N/2)*2/N); axis([0,200,0,0.01]);
title('解扩后信号功率谱'); xlabel('Hz');
如图，解扩后信号的频谱被压缩，功率幅度增加，符合理论分析结果。
- 21 -
14、对解扩信号进行采样、判决
figure(14) subplot(2,1,1) for i=1:1:350
ij=i*1000-500; ss(i)=yrej(ij); end
stem(ss);
title('解扩信号采样'); subplot(2,1,2) for i=1:1:350 %判决信号算法 if ss(i)>0.2 ss(i)=1; elseif ss(i)<-0.2 ss(i)=-1; else ss(i)=0; end end
for i=1:1:50 ij=7*i-6; if ss(ij)==0
ss(ij)=ss(ij+4); end end
for i=1:1:348 if ss(i)==0
ss(i)=ss(i+2); end end
for i=1:1:50
S(i)=ss(i*7-3); if S(i)==0 S(i)=S(i)+1; end
S(i)=(1-S(i))/2; end stem(S);
title('判决后的最终信号');
- 22 -
至此，扩频系统仿真完成，比较最终判决信号与原发送信号完全吻合。
- 23 -
（2）信道中加入窄带强干扰
加入窄带干扰的matlab系统流程图
100/7Hz 二进制比特信息 100Hz 7位双极性m序列 100Hz 扩频序列 2000Hz 载波cos4000πt BPSK调制信号 高斯白噪声 2040~2050Hz 窄带强干扰 恢复载波cos4000πt 凯萨尔滤波器低通滤波 100Hz 7位双极性m序列 采样、判决 - 24 -
15、在信道中加入2040~2050Hz窄带强干扰并乘以恢复载波
产生方法为：生成单位冲激信号，其频谱为白色。之后，用2040~2050Hz的带通滤波器对其进行滤波。之后，与恢复载波相乘。
figure(15) subplot(2,1,1) fd=200000;
Wp1=2*2040/fd; 40~2050Hz带通滤波器 Wp2=2*2050/fd; Wc1=2*2030/fd; Wc2=2*2060/fd; Ap=1; As=100;
W1=(Wp1+Wc1)/2; W2=(Wp2+Wc2)/2;
wdth=min((Wp1-Wc1),(Wc2-Wp2)); Nd=ceil(11*pi/wdth)+1; bd=fir1(Nd,[W1 W2]); zd(1)=1; %产生单位冲激信号 for i=2:1:350000 zd(i)=0; end
ds=abs(freqz(bd,1,400000,fd))';
ybz=fft(zd,N)*100000; magz=abs(ybz);
dz=ds.*magz; %冲激信号频谱乘以带通滤波器 dsz=maga+dz; %magrelz=magrel;
plot(freb,dz(1:N/2)*2/N,freb,maga(1:N/2)*2/N); xlabel('Hz');
axis([1700,2300,0,0.6]);
title('经信道后，加上一窄带强干扰的频谱'); subplot(2,1,2)
rez=real(ifft(dz,N)); ts=(1:N)/N*4;
yzz=rez.*cos(2*pi*2000*ts); %信号乘以恢复载波 yz=fft(yzz,N); magyz=abs(yz);
renz=real(ifft(maga,N));
ynzz=renz.*cos(2*pi*2000*ts);
- 25 -
ynz=fft(ynzz,N); magynz=abs(ynz);
plot(freb,magyz(1:N/2)*2/N,freb,magynz(1:N/2)*2/N); axis([0,5000,0,0.2]);
title('加窄带干扰的信号乘以恢复载波后频谱');
如图，绿色部分为原扩频信号，蓝色部分为新加上窄带干扰频谱。从图上可以看出，窄带干扰非常强。
16、对加窄带干扰的信号进行低通滤波并解扩
figure(16) subplot(3,1,1)
magyzl=bs.*magyz; %用凯赛尔窗低通滤波器滤波 magynzl=bs.*magynz;
plot(freb,magyzl(1:N/2)*2/N,freb,magynzl(1:N/2)*2/N); axis([0,200,0,0.2]); xlabel('Hz');
title('信号与窄带干扰经过凯萨尔窗低通滤波后频谱');
- 26 -
subplot(3,1,2)
yrnzl=real(ifft(bs.*yre,400000)); %对滤波后带干扰信号进行ifft变换 yrzl=real(ifft(magynzl,400000)); tm=(1:N)/N*4;
yrnzlj=jj.*yrnzl(1:350000); yrzlj=jj.*yrzl(1:350000);
plot(ts(1:350000),yrnzlj+yrzlj); xlabel('t');
axis([0,4,-1.5,1.5]);
title('加窄带干扰后解扩信号波形');
subplot(3,1,3) %分别画出信号与干扰经低通滤波后频谱 yzj=fft(yrzlj,N); magzj=abs(yzj); ynzj=fft(yrnzlj,N); magnzj=abs(ynzj);
plot(freb,magzj(1:N/2)*2/N,freb,magnzj(1:N/2)*2/N); axis([0,500,0,0.2]);
title('信号与窄带干扰解扩后频谱'); xlabel('Hz');
能够看出，由于窄带干扰与扩频序列不相关，经解扩后频谱幅度明显降低，对信号无太大影响，信号的时域波形只是增加了幅度很小的干扰，对采样判决没有影响。说明扩频系统能够很好抑制窄带强干扰。
- 27 -
17、比较解扩后信号与窄带强干扰的功率谱
figure(17)
prelnz=ynzj.*conj(ynzj)/N; prelz=yzj.*conj(yzj)/N;
plot(freb,prelnz(1:N/2)*2/N,freb,prelz(1:N/2)*2/N); axis([0,100,0,0.007]); xlabel('Hz');
title('信号与窄带干扰经解扩后的功率谱');
比较功率谱。蓝色部分为信号功率谱，绿色部分为窄带强干扰功率谱，可以发现窄带强干扰已经完全淹没在信号功率里。由理论计算可知，信噪比与N=1无扩频相比能够提高N2倍，抗干扰能力明显增强。
- 28 -
六、误码率simulink仿真
误码率仿真采用了MATLAB里面的SIMULINK及BERTOOL工具。SIMULINK是MATLAB最重要的组件之一，它提供一个动态系统建模、仿真和综合分析的集成环境。SIMULINK具有适应面广、结构和流程清晰及仿真精细、贴近实际、效率高、灵活等优点，并基于以上优点SIMULINK已被广泛应用于控制理论和数字信号处理的复杂仿真和设计。同时有大量的第三方软件和硬件可应用于或被要求应用于SIMULINK。
1、直接扩频系统信道模型
- 29 -
2、加窄带干扰的直扩系统建模
仿真过程中，仍然采用了100Hz的扩频序列，100/7Hz的数据码流，2kHz的BPSK调制并加入了与中心频点频偏20Hz的窄带强干扰。误码率仿真采用了SIMULIMK自带的Error Rate Calculation模块，来对比发送与接收的码流。建模的过程中，考虑了信号的频率，极性，窄带干扰和awgn信道的高斯白噪声，完全符合一般的通信系统的组成结构。仿真时间为100s。
- 30 -
3、用示波器观察发送码字及解扩后码字
上图为仿真过程中截取的部分发送与接收码字，上半部分为发送码字序列，下半部分为最终解扩后的码字序列。可以看出，大多数时间信号吻合，在第25.86s左后出现两个误码。
- 31 -
4、直接扩频系统与无扩频系统的误码率比较
采用SIMULINK下的BERTOOL工具，可以轻松地仿真出信号的误码率。在误码率计算中，我分别仿真了不同m序列长度和不同窄带干扰强度下，误码率与awgn高斯信道信噪比的关系图。
下图为无窄带干扰时，无扩频系统与N=7的m序列直接扩频BPSK系统的误码率比较，横坐标为信道信噪比。
可以看出，扩频序列误码率与无扩频系统的理论误码率基本相同，说明在没有窄带干扰的情况下扩频与否对于误码率影响不大。这说明，在干扰为高斯白噪声的情况下，扩频系统与无扩频系统的抗干扰能力相同。
- 32 -
5、不同扩频序列长度下的误码率比较
如图，扩频系统的误码率与扩频序列的长度有关，可以看出，在信噪比很小的情况下，不同扩频序列长度的误码率基本相同。当信噪比逐渐增大时，系统的误码情况有所好转。扩频码越长，误码率越低，抗干扰能力越强。
- 33 -
6、扩频序列长度N=7时，不同强度窄带干扰下的误码率比较
如图，当窄带强干扰的振幅为信号幅度的10倍时，扩频系统的误码率不变，说明扩频系统对窄带强干扰有很强的抑制能力，符合之前的matlab仿真结果。当窄带强干扰的振幅增至信号幅度的250倍和500倍时，系统误码率随干扰振幅增大而增大。窄带干扰幅度为信号幅度250倍时，误码率小于0.01，可以接受。当窄带干扰幅度为信号幅度的500倍时，误码率较大，不能忍受，但这种情况基本不会出现。
- 34 -
七、利用Walsh码实现码分多址技术
沃尔什码的自相关特性和互相关特性都不理想，这意味着CDM信号经过多径信道传输时，每个用户的不同径之间会产生严重干扰，同时不同用户不同径之间也会产生严重干扰。为此，可以用相关性较好的伪随机序列，即m序列与沃尔什码模二加，得到改善的沃尔什码。这样，既保留了沃尔什码的正交特性，又大大改善了其相关特性。
12.5Hz二进制比特信息1 12.5Hz二进制比特信息2 100Hz8位改善walsh码m1 100Hz8位改善walsh码m2 100Hz 扩频序列 2000Hz 载波cos4000πt BPSK调制信号 高斯白噪声 2040~2050Hz 窄带强干扰 恢复载波cos4000πt 凯萨尔滤波器低通滤波 100Hz8位改善walsh码1 采样、判决信息1 - 35 -
100Hz8位改善walsh码m2 采样、判决信息2
1、产生改善的walsh码
在本实验中，我采用了前面的7位m序列改善的沃尔什码。产生方法为在双极性m序列后加“1”，并与正交沃尔什码相乘，并用改善后的沃尔什码对两路信息序列分别进行扩频，相加后共同传输，实现码分多址技术。
由于m序列为7位，补位后0和1出现概率相等，成为了真正的随机序列。两路沃尔什函数采用了八阶沃尔什函数的的第一个和第五个，分别为1,1,1,1,1,1,1,1和1,1,1,1,0,0,0,0，为正交码。
clear all; clc;
X1=0;X2=0;X3=1; m=7; for i=1:m
Y3=X3; Y2=X2; Y1=X1; X3=Y2; X2=Y1; X1=xor(Y3,Y1); L(i)=Y1; end
L(8)=0; %单极性7位m序列末尾补0 m=m+1; for i=1:m
M(i)=1-2*L(i); end
k=1:1:m; figure(1)
subplot(3,1,1) stem(k-1,M); axis([0,9,-1,1]); xlabel('k');
title('移位寄存器产生的双极性7位M序列,末尾加一个1') ; subplot(3,1,2)
for i=1:m %生成改善的沃尔什码m1 m1(i)=M(i); end
stem(k-1,m1); axis([0,8,-1,1]); xlabel('k');
title('改善的沃尔什码m1') ; subplot(3,1,3)
for i=1:4 %生成改善的沃尔什码m2 m2(i)=M(i); end
for i=5:8
- 36 -
m2(i)=-M(i); end
stem(k-1,m2); axis([0,8,-1,1]); xlabel('k');
title('改善的沃尔什码m2') ;
如图，改善的沃尔什码m1与m2正交。下面将用这两路沃尔什码分别对两路信息进行扩频调制。
2、产生两路不同的信息序列
figure(2) N=50;a=0;
x_rand=rand(1,N); for i=1:N
if x_rand(i)>=0.5 x1(i)=1;a=a+1;
- 37 -
else x1(i)=0; end end
t=0:N-1;
subplot(2,1,1) stem(t*0.08,x1); xlabel('t/s');
title('扩频前待发送二进制信息序列1'); x_rand=rand(1,N); for i=1:N
if x_rand(i)>=0.5 x2(i)=1;a=a+1; else x2(i)=0; end end
subplot(2,1,2) stem(t*0.08,x2);
title('扩频前待发送二进制信息序列2'); xlabel('t/s');
- 38 -
3、用两个沃尔什码分别调制两路信号
figure(3)
y1=rectpulse(x1,8);y2=rectpulse(x2,8); for i=1:1:50
M1(8*i-7)=m1(1);M2(8*i-7)=m2(1); M1(8*i-6)=m1(2);M2(8*i-6)=m2(2); M1(8*i-5)=m1(3);M2(8*i-5)=m2(3); M1(8*i-4)=m1(4);M2(8*i-4)=m2(4); M1(8*i-3)=m1(5);M2(8*i-3)=m2(5); M1(8*i-2)=m1(6);M2(8*i-2)=m2(6); M1(8*i-1)=m1(7);M2(8*i-1)=m2(7); M1(8*i)=m1(8);M2(8*i)=m2(8); end
subplot(2,1,1) tt=0:8*N-1; for i=1:400;
x11(i)=1-2*y1(i);x12(i)=1-2*y2(i);s1(i)=M1(i)*x11(i);s2(i)=M2(i)*x12(i); end
subplot(2,1,1) stem(tt/100,s1);
title('扩频后的待发送序列码1'); xlabel('t/s'); subplot(2,1,2) stem(tt/100,s2);
title('扩频后的待发送序列码2'); xlabel('t/s');
- 39 -
4、两路信号相加，并进行BPSK调制 figure(4)
subplot(2,1,1) s=s1+s2;
stem(tt/100,s);
title('两路信息相加,即进行walsh码复用'); xlabel('t/s'); subplot(2,1,2) fs=2000;
ts=0:0.00001:4-0.00001; s_b=rectpulse(s,1000);
s_bpsk=s_b.*cos(2*pi*fs*ts); plot(ts,s_bpsk); xlabel('t/s');
axis([0.065,0.095,-2.2,2.2])
title('walsh码分复用做BPSK变换');
对比walsh码进行复用后的信号与无码分多址系统的信号，可以看出，无码分多址系统的信息码只有两个取值-1和1，但是复用后存在三个值-1,0,1，所以BPSK调制信号存在为零的时刻。
- 40 -
5、观察调制信号频谱，并经awgn信道加高斯白噪和窄带强干扰
figure(5)
subplot(3,1,1) N=400000;
yb=fft(s_bpsk,N); mag=abs(yb);
fb=(1:N/2)*100000/N; plot(fb,mag(1:N/2)*2/N); axis([1000,3000,0,0.25]);
title('码分复用扩频后调制信号频谱'); xlabel('f/Hz'); subplot(3,1,2)
s_bpska=awgn(s_bpsk,3,'measured');%经过awgn信号，信噪比3dB plot(ts,s_bpsk,ts,s_bpska);
title('码分复用扩频后经加噪过信道后的信号与原信号时域波形对比'); xlabel('t/s');
axis([0.0775,0.0825,-2.4,2.4]); subplot(3,1,3)
yba=fft(s_bpska,N); maga=abs(yba);
fb=(1:N/2)*100000/N; fd=200000;
Wp1=2*2040/fd;%%%%%%%%%%带通滤波器 Wp2=2*2050/fd; Wc1=2*2030/fd; Wc2=2*2060/fd; Ap=1; As=100;
W1=(Wp1+Wc1)/2; W2=(Wp2+Wc2)/2;
wdth=min((Wp1-Wc1),(Wc2-Wp2)); Nd=ceil(11*pi/wdth)+1; bd=fir1(Nd,[W1 W2]); zd(1)=1;%单位冲激 for i=2:1:350000 zd(i)=0; end
ds=abs(freqz(bd,1,400000,fd))'; ybz=fft(zd,N)*40000; magz=abs(ybz); dz=ds.*magz;
dsz=maga+dz;%将窄带干扰叠加到信号上
- 41 -
plot(fb,dsz(1:N/2)*2/N); axis([1000,3000,0,0.4]); xlabel('f/Hz');
title('码分复用扩频后经信道调制加窄带强干扰信号频谱');
- 42 -
6、接收机信号乘以恢复载波，观察时域和频域
figure(6)
subplot(2,1,1)
rez=real(ifft(dz,N));
re=(s_bpska+rez).*cos(2*pi*fs*ts); plot(ts,re);
axis([0.065,0.095,-2.8,2.8]); xlabel('t');
title('扩频后接收信号乘以恢复载波'); xlabel('t/s'); subplot(2,1,2) yre=fft(re,N); magre=abs(yre);
plot(fb,magre(1:N/2)*2/N);
title('码分复用扩频后乘以恢复载波后信号频谱'); axis([0,5000,0,0.16]); xlabel('f/Hz');
- 43 -
7、信号经凯萨尔窗低通滤波器
figure(7)
subplot(2,1,1)
magrel=bs.*magre;
plot(fb,magrel(1:N/2)*2/N);
title('码分复用信号经过凯萨尔窗函数低通滤波'); axis([0,200,0,0.25]); xlabel('f/Hz'); subplot(2,1,2) tm=(1:N)/N*4;
yrel=real(ifft(bs.*yre,400000)); plot(tm,yrel); xlabel('t/s');
title('经过凯萨尔窗函数滤波后时域波形');
可以看出，信号经过滤波后的时域波形前后部分畸变较大，说明用改善的walsh码进行码分复用时，不同用户间仍然后一定的码间干扰。
- 44 -
8、对滤波后信号分别用m1和m2进行解扩
figure(8)
subplot(2,1,1)
j1=rectpulse(M1,1000); j2=rectpulse(M2,1000); yrej1=j1.*yrel; plot(ts,yrej1);
title('时域波形与解扩码m1相乘波形'); xlabel('t/s'); subplot(2,1,2) yrej2=j2.*yrel; plot(ts,yrej2); xlabel('t/s');
title('时域波形与解扩码m2相乘波形');
对比码分复用系统和无码分复用系统，可以看出，解扩后波形出现了0项，每个有效样值的持续时间比无码分复用的系统减少了一半。而且解扩信号时域波形前后部分的干扰比解扩前还要严重，说明码分复用引入了码间干扰。
- 45 -
9、对两路信号分别采样，判决
figure(9) for i=1:1:50
S1(i)=0;S2(i)=0; end
for i=1:1:800
ij=i*500-250;
ss1(i)=yrej1(ij);ss2(i)=yrej2(ij); end
for i=1:1:800 if ss1(i)>0.3 ss1(i)=1; elseif ss1(i)<-0.3 ss1(i)=-1; else ss1(i)=0; end end
for i=1:1:800 if ss2(i)>0.3 ss2(i)=1; elseif ss2(i)<-0.3 ss2(i)=-1; else ss2(i)=0; end end
n1=1;n2=1; for ii=1:1:100 i=ii*8-7;
k1=ss1(i)+ss1(i+1)+ss1(i+2)+ss1(i+3)+ss1(i+4)+ss1(i+5)+ss1(i+6)+ss1(i+7); k2=ss2(i)+ss2(i+1)+ss2(i+2)+ss2(i+3)+ss2(i+4)+ss2(i+5)+ss2(i+6)+ss2(i+7); if k1>=5
S1(n1)=0;n1=n1+1; elseif k1<=-5
S1(n1)=1;n1=n1+1; end if k2>=5
S2(n2)=0;n2=n2+1; elseif k2<=-5
S2(n2)=1;n2=n2+1; end end
subplot(2,1,1)
- 46 -
stem(S1);
title('序列1判决后的最终信号'); subplot(2,1,2) stem(S2);
title('序列2判决后的最终信号');
比较原信息序列，可以看出最终判决信号与原信号相同，但是有时会出现误码。所以，为了防止多用户带来的码间干扰，CDMA除了进行码分复用以外还应采用其他措施抗干扰。至此，walsh正交码的码分复用仿真完成。
- 47 -
八、产生随机序列Gold码和正交Gold码
Gold码是由m序列派生出的一种伪随机码，它具有类似于m序列具有的伪随机性质，但其同长度不同序列的数目比m序列的多得多。
Gold码发生器框图：
m1 Gold码
Gold码是由m序列的优选对移位模二加构成，m1和m2为同长度的两个不同m序列并为优选对。Gold码的自相关特性：Gold证明过Gold码序列的自相关函数的所有非最高峰的取值为三值，如下式所示，其中p=2n-1，为Gold码序列的周期。
?1??p??tR???
?p?t?2??p当n为奇数时，t?2n?12m2 +1，当n为偶数且不是4的整数倍时，t?2n?12?1。
实验中，我采用了长度n=7的一对优选对构成Gold码。采用matlab里面的gfprimfd(7,'all')命令得到全部的18个本原多项式表达式，如下图：
之后选取的本原多项式的八进制数表示为m1（211），m2（217）。
- 48 -
1、产生Gold码并仿真其自相关函数
clear all; clc;
X1=1;X2=1;X3=1;X4=1;X5=1;X6=1;X7=-1; m=2^7-1; for i=1:1:m
Y7=X7;Y6=X6;Y5=X5;Y4=X4;Y3=X3;Y2=X2;Y1=X1; X7=Y6;X6=Y5;X5=Y4;X4=Y3;X3=Y2;X2=Y1; X1=Y3*Y7; L1(i)=Y1; end
for i=1:1:m
Y7=X7;Y6=X6;Y5=X5;Y4=X4;Y3=X3;Y2=X2;Y1=X1; X7=Y6;X6=Y5;X5=Y4;X4=Y3;X3=Y2;X2=Y1; X1=Y1*Y2*Y3*Y7; L2(i)=Y1; end
for i=1:1:m
L(i)=L1(i)*L2(i); L(i)=1-L(i); L(i)=0.5*L(i); end
figure(1)
subplot(2,1,1) i=1:127; stairs(i,L);
axis([0,130,-1.2,1.2]); title('n=7Gold码'); xlabel('k'); subplot(2,1,2)
a=L;b=a;Nb=127;N=2*m; for k=1:N
c=xor(a,b); D=sum(c); A=Nb-D;
R(k)=(A-D)/(A+D);%计算相关系数 b=[b(Nb),b(1:(Nb-1))]; end
k=0:(N-1); plot(k,R,'--r.');
axis([64,191,-0.5,1.2]);
title('n=7Gold码自相关函数');
- 49 -
xlabel('k');
?1??p?t计算出，除了在k=127*n即序列周期的整数倍时自相关可以根据R?????p?t?2??p函数为1以外，其余自相关数值应为三值函数，取值分别为-0.0787,-0.1339和0.1181。而从图上看出，Gold自相关函数满足以上条件，具有尖锐的自相关特性，说明仿真正确。
- 50 -