bpsk1.m

利用matlab对BPSK系统的误码率进行仿真

%本程序中时间单位是微秒
%频率单位为MHz
%码元速率单位是Mb/s

global dt df t f N
close all
clear Eb_N0 Pe



N=2^13;   %采样点数
L=64;    %每码元的采样点数
M=N/L     %码元数
Rb=2;
Ts=0.5		%码元宽度是0.5us
dt=Ts/L;        
df=1/(N*dt)     %MHz
T=N*dt          %截短时间
Bs=N*df/2       %系统带宽
Na=4;           %示波器扫描宽度为4个码元
alpha=1;
fc=20;

t=[-T/2+dt/2:dt:T/2]; %时域横坐标
f=[-Bs+df/2:df:Bs];   %频域横坐标


g1=sin(pi*t/Ts)./(pi*t/Ts);
g2=cos(alpha*pi*t/Ts)./(1-(2*alpha*t/Ts).^2);
g=g1.*g2;		 %升余弦脉冲波形
G=abs(t2f(g));
GT=sqrt(G);
GR=GT;

h1=cos(2*pi*fc*t)
h2=h1;



figure(1)
 set(1,'Position',[10,50,300,200])
 
                %设定窗口位置及大小 
figure(2)
 set(2,'Position',[400,50,300,200])
                %设定窗口位置及大小 
                
 hold on
 grid 
 xlabel('t in us')
 ylabel('s(t) in V')

figure(3)
 set(3,'Position',[10,50,500,400])
                %设定窗口位置及大小 
figure(4)
 set(4,'Position',[400,50,500,400])
 
                %设定窗口位置及大小 
figure(5)
 set(5,'Position',[10,50,100,200])
                %设定窗口位置及大小 
                
                                                                                                                                                                                                                                                                                                             
EP=zeros(size(f));

for loop1=1:20
 Eb_N0(loop1)=(loop1-1)  %Eb/N0 in dB
 eb_n0(loop1)=10^(Eb_N0(loop1)/10);
 Eb=1;	
 n0=Eb/eb_n0(loop1); %信道的噪声谱密度   
 sita=n0*Bs; %信道中噪声功率

 n_err=0;     %误码计数

 for loop2=1:5 
 a=sign(randn(1,M));
 imp=zeros(1,N);    %产生冲激序列
 imp(L/2:L:N)=a/dt;
 IMP=t2f(imp)
 
S=t2f(imp).*G       %升余弦信号的傅氏变换
ss=f2t(t2f(imp).*G);  %升余弦信号的时域波形
s=real(ss.*h1);
ss=real(ss);
s1=s.*h2
S=t2f(s);
 
P=S.*conj(S)/T;      %升余弦信号的功率谱
EP=(EP*(loop2-1)+P+eps)/loop2;
 
n_ch=sqrt(sita)*randn(size(t)); %信道噪声
nr=real(f2t(t2f(n_ch).*GR));    %输出噪声
sr=real(s1)+nr;       %接收信号                        

y=sr(L/2:L:N);		  %取样
aa=sign(y);       %判决                 
n_err=n_err+length(find(aa~=a))
 
 
c=zeros(L,M);
for i=1:L,c(i,:)=aa;
    
end
 
c=reshape(c,1,N);
 
 
end
Pe(loop1)=n_err/(M*loop2)
 
if loop1==10
                                                                                 
figure(1)                               %功率谱
plot(f,30+10*log10(EP),'g');
grid
axis([-40,+40,-50,50])
xlabel('f (MHz)')
ylabel('Ps(f)  (dBm/MHz)')
 
 
figure(3)                               %调制信号
plot(t,s,'b')
grid
axis([-2,+2,-1.5,1.5])
xlabel('t (us)')
ylabel('s(t) (V)')
 
 
end


figure(4)                               %解调信号                           
plot(t,c,'b')  
grid
axis([-2,+2,-1.5,1.5])
xlabel('t (us)')
ylabel('s(t) (V)')


figure(5)                                %误码率
semilogy(Eb_N0,Pe,'r');
eb_n0=10.^(Eb_N0/10);
hold on
semilogy(Eb_N0,0.5*erfc(sqrt(eb_n0)));
axis([0,9,1e-4,1])
xlabel('Eb/N0')
ylabel('Pe')

figure(2)                               %眼图
tt=[0:dt:Na*L*dt];
for jj=1:Na*L:N-Na*L
plot(tt,ss(jj:jj+Na*L));
end

 
end