BPSK系统仿真—MATLAB

程序代码：
主程序：
%基带信号s(t)
clear all
close all
clc
num=10; %码元个数
tnum=200;%码元长度
N=num*tnum;
a=randint(1,num,2); %产生1行num列的矩阵，矩阵内0和1随机出现
fc=0.1; %载波频率
t=0:0.05:9.99;
s=[];c=[];
for i=1:num
if(a(i)==0)
A=zeros(1,tnum);
else
A=ones(1,tnum);
end
s=[s A];
cs=sin(2*pi*fc*t);
c=[c cs]; %产生载波信号
end
figure(1);subplot(4,1,1);plot(s);grid on;
xlabel('基带信号s(t)');
axis([0 N -2 2]);

%BPSK已调信号。采用模拟调制方法
s_NRZ=[];
for i=1:num
if(a(i)==0)
A=ones(1,tnum);
else
A=-1*ones(1,tnum);
end
s_NRZ=[s_NRZ,A];
end
e=s_NRZ.*c;
figure(1);subplot(4,1,2);plot(e);grid on;
xlabel('BPSK调制信号');
axis([0 N -2 2]);

%加高斯噪声
am=0.7; %输入信号经信道后振幅由1衰减为0.7
SNR=5;
snr=10^(SNR/10);
N0=(am*am)/2/snr;%计算噪声功率
N0_db=10*log10(N0);%将噪声功率转换为dBW
ni=wgn(1,N,N0_db);% 产生1行N列的高斯噪声
yi=e+ni;%BSK已调信号中加入白噪声，输入信噪比为SNR
figure(1);subplot(4,1,3);plot(yi);grid on;
xlabel('加入高斯白噪声的已调信号yi(t)');

% %带通滤波器
[b1,a1] = BUTTER(3,[2*pi*0.0001,2*pi*0.01]);%计算带通滤波器的H(z)系数
y=filter(b1,a1,yi);%对信号yi进行滤波，得到信号y
figure(1);subplot(4,1,4);plot(y);grid on;
xlabel('经带通滤波器后信号');

%与恢复载波相乘
x1=2*c.*y;
figure(2);subplot(4,1,1);plot(x1);grid on
xlabel('与恢复载波相乘后的信号x1(t)');

%低通滤波器
[b2,a2]=butter(2,0.005); %计算H(z)系数 ，频率为(1/200)
x=filter(b2,a2,x1); %对信号x1滤波，得到信号x
figure(2);subplot(4,1,2);plot(x);grid on
axis([0 N -2 2]);xlabel('经低通滤波器后信号波形')
x=fun_panjue(x);%调用函数，进行抽样判决
figure(2);subplot(4,1,3);plot(x);grid on;
xlabel('加噪后解调信号x(t)');
axis([0 N -2 2]);
x=fun_yanc(x); %调用函数，进行消除延迟
figure(2);subplot(4,1,4);plot(x);grid on;
xlabel('加噪后去掉延迟的解调信号x(t)');
axis([0 N -2 2]);

%误码率计算
Err1=length(find(x~=s))%计算解调信号中错误码元个数
Pe_test1=Err1/N%计算实际误码率
Pe1=(1/2)*erfc(sqrt(snr))%计算系统理论误码率

%% 理论误码率曲线
Pe=[];
for SNR=1:10
am=0.7; %输入信号经信道后振幅由1衰减为0.7
E=am*am/2;
snr=10^(SNR/10);
N0=(am*am)/(2*snr);
no=N0/(2*200); %计算噪声功率
N0_db=10*log10(N0);%将噪声功率转换为dBW
ni=wgn(1,N,N0_db);% 产生1行N列的高斯噪声
yi=e+ni; %BSK已调信号中加入白噪声，输入信噪比为SNR
y=filter(b1,a1,yi);%对yi进行滤波(带通滤波器)，得到信号y
x1=2*c.*y; %与恢复载波相乘
xx=filter(b2,a2,x1); %经低通滤波器滤波
xx=fun_panjue(xx);%抽样判决
xx=fun_yanc(xx); %消除延迟
snr=10^(SNR/10);
Pe=[Pe,(1/2)*erfc(sqrt(snr))];%计算理论误码率
end
Pe;
figure;
SNR=1:10;
semilogy(SNR,Pe,'b--');hold on %以log10(Pe)为纵坐标画图
grid on

子程序一：抽样判决
（由于存在“倒pi”现象，故以0为基准进行判决，大于0，判为0；小于0，判为1）
function w=fun_panjue(w)
N=length(w);
if w(100)>0
w(1:100)=0;
else
w(1:100)=1;
end
for i=101:N
if w(i)>0;
w(i)=0;
else
w(i)=1;
end
end

子函数二：消除延迟
（由于抽样判决后会产生延迟，故认为对解调信号进行时移，以减小甚至消除延迟，便于与
基带信号进行比较，计算实际误码率）
function m=fun_yanc(m)
N=length(m);
leng=0;
if m(1)==0
for i=1:N
if m(i)==1
leng=i;
break;
end
end
else
for i=1:N
if m(i)==0
leng=i;
break;
end
end
end
leng1=leng-(floor(leng/200))*200;
for i=1:(N-leng1)
m(i)=m(i+leng1);
end
for i=(N-leng1):N
m(i)=m(N-200+10);
end
运行结果：
输入信噪比为5 dB时，计算的错误码元，实际误码率，理论误码率
Err1 =

15

Pe_test1 =
0.0075

Pe1 =
0.0060
各波形图如下所示：