bpsk_correlation.m

There are two files in the zip folder. bpsk_spread.m and jakesmodel.m Steps for simulation: 1] Run

clc
clear all;

%-----------data generation--------------------------
fs=2400;
ts=1/fs;
data_L=24000;          %产生2400个数据
databits=randsrc(data_L,1);
%------------------差分编码---------------------------
ai=(databits+1)/2;
bi=zeros(1,data_L);
bb=0;     %假设初始参考相位为0;
for k=1:data_L
    if k==1
       bi(1)=xor(ai(1),bb);       
    else
      bi(k)=xor(ai(k),bi(k-1));     
    end
end
bi1=2*bi-1;
dataI=bi1;
%-----------------upsample---------------------------
dataI=repmat(dataI,40,1);   %将数据上采样到射频caiyang频率
dataI=reshape(dataI,1,40*length(dataI));
%----------------波形成形----risecosine---------------------
td=ts/40;
Fd=1/td;
Fs=2*Fd;
sampletime=Fs/Fd*40;
Delay=3;
alpha=0.35;
I0=rcosflt(dataI,Fd,Fs,'fir',alpha,Delay);   %对数据比特流进行波形成形，目的是压缩带宽。
I0(1:Delay*Fs/Fd-fix(Fs/Fd/2))=[];          %经过滤波器以后波形会有一段延迟，这里截去前边延迟的无用数据，
                                            %虽然实际中并不截去，但是为了处理方便，这里截去无用数据
% ----- --- --zaipin generation-----------------------
fc=2/ts;    %载频是4800hz，说明一个码元内有2个完整的载波波形  
w=2*pi*fc; %--仿真时发现，载波频率越高（此时采样点还是1:length(I0)个，并未提高），误码率越低，为什么？
n=(1:length(I0))*(td/2);  %将载波离散化，间隔和数目应与数据一致
zaipin1=cos(w*n+pi/4);
%---------Qpsk modulation--------------------------
dataImodulate=I0'.*zaipin1;  %I路BPSK调制
finalsenddata=dataImodulate;
%-------------figure------------------------------
figure
subplot(2,1,1)
plot(dataI)
title('原始数据I差分编码后')
axis([0 300 -1.5 1.5])
subplot(2,1,2)
plot(dataImodulate)
title('I路bpsk调制信号')
axis([0 600 -1.5 1.5])
grid on

%------------- Design filter.-----------------------
%因为要滤掉2w频段，所以截止频率要<2*w/(2*pi)=9600hz;这里设计成4000hz和4800hz
fpass=2400;       %通带截止频率  
fstop=3000;       %阻带截止频率
bandwith=fc;
wp=fpass/fc/2;  %归一化
ws=fstop/fc/2;  
Rp=3;           %通带允许最大衰减
Rs=50;          %阻带允许最小衰减
[stepnum,wn]=buttord(wp,ws,Rp,Rs); %设计巴特沃思滤波器
[fenzi,fenmu]=butter(stepnum,wn);
% H=freqz(fenzi,fenmu,100,100);
% figure
% plot(H)
%---------------加噪声---------------------
snr_in_dB=[6 8 10 12 15];
L=3;   %对每个信噪比要计算L遍pe
pe=zeros(1,length(snr_in_dB));%初始化
Ipower=sum(dataImodulate.^2)/length(dataImodulate); %I路能量

signalpower=Ipower;   %数据能量

for q=1:length(snr_in_dB)
   for p=1:L
snr=10^(snr_in_dB(q)/10);
N_hope=signalpower/snr;       
nnnn=randn(2222,1);        
noise=randn(1,length(finalsenddata)); %产生噪声数据
noiseshape=filter(fenzi,fenmu,noise);  %限制噪声带宽与数据一致
N_power=sum(noiseshape.^2)/length(noiseshape);
bizhi=sqrt(N_hope/N_power);
noise=bizhi*noiseshape;
receivedata=finalsenddata+noise;    %最终发送数据
%----------------------------------------------------
%-------相干解调----------demodulation-----------------------
%----------------------------------------------------
% --------------Upsample and demodulate--------------
m=(1:length(receivedata))*td/2;
I1=receivedata.*cos(w*m+pi/4);         %相干解调
%-------------- filter----------------------
I2=filter(fenzi,fenmu,I1);           %低通滤波，滤掉2*w分量
I3=I2;
I3(1:fix(stepnum/3))=[];             %截去延迟的无用数据
% figure
% subplot(3,1,1)
% stem(I1)
% title('chengyi jietiao zaibo')
% axis([0 50 -1.5 1.5])
% subplot(3,1,2)
% stem(I2)
% title('lowpass')
% axis([0 100 -1.5 1.5])
% subplot(3,1,3)
% stem(I3)
% title('lowpass and cut')
% axis([0 100 -1.5 1.5])
%-----------------downsample and judge---------------
var1=fix(length(I3)/sampletime);
var2=fix(sampletime/2)+1;                 %取中间点
% receivedataI=zeros(var1);
receivedataI=I3((0:var1-1)*sampletime+var2);  %取中间点进行下采样
receivedataI1=sign(receivedataI);
%---------------I-Q-差分码反变换------------------------------
bii=(receivedataI1+1)/2;
bbb=0;
aii=zeros(1,length(receivedataI1));
for k=1:length(receivedataI1)
    if k==1
        aii(1)=xor(bii(1),bbb);    
    else
        aii(k)=xor(bii(k),bii(k-1));    
    end
end
aii1=2*aii-1;
receivedataI2=aii1;
 %----------------判决-------------------------------
 finaldataout=sign(receivedataI2);
 %-------------计算误码率-----------------
  errornum=0;    %初始化，以防万一；
  sumpe=[];
  pe1=0;
 errornum=nnz(finaldataout(1:data_L)-databits');  %计算误码个数
 pe1=errornum/data_L;
 sumpe(p)=pe1;
   end
   pe(q)=sum(sumpe)/L;
end
 disp('pe=')
 disp(pe)
 %---------------画图------------------------
 figure
 subplot(2,1,1)
stem(databits)
title('原始数据')
axis([0 64 -1.5 1.5])
subplot(2,1,2)
stem(finaldataout)
 title('解调数据')
 axis([0 64 -1.5 1.5])
 figure
 semilogy(snr_in_dB,pe,'-x')
 title('误码率曲线')
 grid