%**************************m序列1************************
clear;
%a=[1 0 0 0 0 1 0 0 1 0 1 1 0 0 1 1 1 1 1 0 0 0 1 1 0 1 1 1 0 1 0];
fbconnection1=[ 0 0 1 0 1 ];
register1 = [0 0 0 0 1];
n1 = length(fbconnection1);
N1 = 2^n1-1;
mseq1(1)= register1(n1);
for i = 2:N1
newregister1(1)= mod(sum(fbconnection1.*register1),2); %寄存器与反馈的模2加
for j = 2:n1
newregister1(j)= register1(j-1);
end;
register1 = newregister1;
mseq1(i) = register1(n1);
end
figure(1)
subplot(5,1,1);
stem(mseq1)
axis([0 35 -0.2 1.2]);
title('m序列1');
grid on;
%*****************m1自相关********************
L=length(mseq1);
m1=mseq1;
N1=120;
for i=1:1:N1
output1(i)=sum((1-2*m1).*(1-2*mseq1));
m1=circshift(m1',1)'; %加一撇表示转置,一个参数默认对行进行移位,等价于circshift(m1,[0,1])
subplot(5,1,2);
end
plot(1:N1,output1(1:N1),'g');
title('m1序列自相关');
%**************m序列2******************
%b=[1 1 1 1 1 0 1 1 1 0 0 0 1 0 1 0 1 1 0 1 0 0 0 0 1 1 0 0 1 0 0];
fbconnection2=[ 1 1 1 0 1 ];
register2 = [1 1 1 1 1];
n2 = length(fbconnection2);
N2 = 2^n2-1;
mseq2(1)= register2(n2);
for i = 2:N2
newregister2(1)= mod(sum(fbconnection2.*register2),2);
for j = 2:n2
newregister2(j)= register2(j-1);
end;
register2 = newregister2;
mseq2(i) = register2(n2);
end
subplot(5,1,3);
stem(mseq2) %绘制离散函数
axis([0 35 -0.2 1.2]);%设置坐标轴
title('m序列2');
grid on;
%*********************m2自相关*********************
L=length(mseq2);
m2=mseq1;
N2=120;
for i=1:1:N2
output1(i)=sum((1-2*m1).*(1-2*mseq1));
m=output1(i);
m1=circshift(m1',1)';
subplot(5,1,4);
end
plot(1:N2,output1(1:N2),'g');
title('m2序列自相关');
%***************m1,m2的互相关函数******************
L=length(mseq1);
m=mseq1;
N3=120;
for k=1:1:N3
output2(k)=sum((1-2*m).*(1-2*mseq1));
m=circshift(m',1)';
subplot(5,1,5);
end
plot(1:N3,output2(1:N3),'g');
title('m1,m2序列互相关');
%******************gold序列*******************
goldseq1=mod(mseq1+mseq2,2);
mseq1=circshift(mseq1',1)';
goldseq2=mod(mseq1+mseq2,2);
figure(2)
subplot(5,1,1);
stem(goldseq1);
axis([0 35 -0.2 1.2]);
title('gold序列1');
subplot(5,1,2);
stem(goldseq2);
axis([0 35 -0.2 1.2]);
title('gold序列2');
grid on;
%******************gold序列自相关******************
L=length(goldseq1);
m3=goldseq1;
N4=160;
for i=1:1:N4
output3(i)=sum((1-2*m3).*(1-2*goldseq1));
m3=circshift(m3',1)';
subplot(5,1,3);
end
%g=output3(1:160);
plot(1:N4,output3(1:N4),'g');
output3(1:N4)
title('gold1序列自相关');
L=length(goldseq2);
m4=goldseq2;
N4=160;
for i=1:1:N4
output4(i)=sum((1-2*m4).*(1-2*goldseq2));
m4=circshift(m4',1)';
subplot(5,1,4);
end
plot(1:N4,output4(1:N4),'g');
title('gold2序列自相关');
%****************gold1,gold2的互相关函数**************
L2=length(goldseq2);
m2=goldseq2;
N5=160;
for k=1:1:N5
output5(k)=sum((1-2*m2).*(1-2*goldseq1));
m2=circshift(m2',1)';
subplot(5,1,5);
end
plot(1:N5,output5(1:N5),'g');
title('gold1,gold2序列互相关');
%*******************功率谱*********************
Fs=1000; %采样频率
n=0:1/Fs:1; %产生含有噪声的序列 xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
nfft=1024; %采样点数
w=blackman(100);
figure(3)
CXk=fft(output3,nfft); %gold序列自相关函数FFT得功率谱
Pxx=abs(CXk);
index=0:round(nfft/2-1); %round函数用于四舍五入取整
k=index*Fs/nfft;
plot_Pxx=(Pxx(index+1));
plot(k,plot_Pxx);
axis([200,250,0,600])
%***********************扩频前序列*******************
N6=100;
i=0;
x_rand=rand(1,N6);
for i=1:1:N6
if x_rand(i)>=0.5
x1(i)=1;
else x1(i)=0;
end
end
t1=0:N6-1;
figure(4)
subplot(3,1,1);
stem(t1,x1);
title('扩频前序列');
%*******************扩频前bpsk信号波形******************
ts=0:0.00001:3.5-0.00001;
fs=2000;
subplot(3,1,2)
s_bb=rectpulse(x1,3500);%将冲击信号补成矩形信号
s_bpskb=(1-2.*s_bb).*cos(2*pi*fs*ts);
plot(ts,s_bpskb);
xlabel('s');
axis([0.026,0.045,-1.2,1.2])
title('扩频前bpsk信号调制信号波形');
%********************扩频后bpsk时域波形*******************
g1=goldseq2;
m2=goldseq1;
m3=[m2 m2 m2 m2 m2 m2 m2 m2 m2 m2 m2 m2];
tt=0:349;
l=1:7*N6;
y(l)=0;
for i=1:N6
k=7*i-6;
y(k)=x1(i);k=k+1;
y(k)=x1(i);k=k+1;
y(k)=x1(i);k=k+1;
y(k)=x1(i);k=k+1;
y(k)=x1(i);k=k+1;
y(k)=x1(i);k=k+1;
y(k)=x1(i);
end
s(l)=0;
for i=1:350
s(i)=xor(m3(i),y(i));
end
subplot(3,1,3)
fs=2000;
ts=0:0.00001:3.5-0.00001;
s_b=rectpulse(s,500);
s_bpsk=(1-2.*s_b).*cos(2*pi*fs*ts);
plot(ts,s_bpsk);
xlabel('t');
axis([0.055,0.085,-1.2,1.2])
title('扩频后bpsk信号调制后的时域波形');
%*************************频谱***********************
figure(5)
N1=400000;
subplot(2,1,1)
yb=fft(s_bpsk,N1);
mag=abs(yb);
fb=(1:N1/2)*100000/N1;
plot(fb,mag(1:N1/2)*2/N1);
axis([1200,2800,0,0.8]);
title('扩频前调制信号频谱');
xlabel('kHz');
N1=400000;
ybb=fft(s_bpskb,N1);
magb=abs(ybb);
fbb=(1:N1/2)*100000/N1;
subplot(2,1,2)
plot(fbb,magb(1:N1/2)*2/N1);
axis([1200,2800,0,0.8]);
title('扩频后调制信号频谱');
xlabel('kHz');
%**********************误码率*********************
frames=2000;%frames为发送的循环次数,次数越多曲线越平滑
EbN0=0:1:10;
SNR = 10.^(EbN0/10);%%信噪比转化为线性值即有dB值转换为数值
slength=100;
% t=0:slength-1;
errtotal=0;
lengthtotal=frames*slength;
pdata=zeros(1,lengthtotal);
for loop= 1:length(SNR)
data=randi([0,1],1,lengthtotal);
bpsk=1-2*data;%bpsk 调制映射,逻辑0映射为pi相位1,逻辑1映射为0相位-1
spow = norm(bpsk).^2/lengthtotal;%求每个符号的平均值,其中norm是求向量2范数函数
sigma = sqrt(spow/(2*SNR(loop)));%根据符号功率求噪声功率
s_receive = bpsk+sigma*(randn(1,length(bpsk)));%添加高斯白噪声;bpsk可以只考虑实部
%解调 硬判决
for i=1:length(s_receive)
if s_receive(i)>=0
s_receive(i)=0;
else
s_receive(i)=1;
end
end
[err,ser(loop)] = symerr(data,s_receive);%误码率
end
figure(6)
ser_theory = qfunc(sqrt(2*SNR));%理论误码率,注意Q函数和误差函数的对应关系
semilogy(EbN0,ser,'-r*',EbN0,ser_theory,'-bo');
title('BPSK信号在AWGN信道下的误码率');
xlabel('EbN0/dB');ylabel('误码率');
legend('未扩频bpsk误码率','bpsk理论误码率');