clc
clear all;

%-----------data generation--------------------------
fs=2400;
ts=1/fs;
data_L=24000;          %產(chǎn)生2400個數(shù)據(jù)
databits=randsrc(data_L,1);
%------------------差分編碼---------------------------
ai=(databits+1)/2;
bi=zeros(1,data_L);
bb=0;     %假設(shè)初始參考相位為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);   %將數(shù)據(jù)上采樣到射頻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);   %對數(shù)據(jù)比特流進(jìn)行波形成形，目的是壓縮帶寬。
I0(1:Delay*Fs/Fd-fix(Fs/Fd/2))=[];          %經(jīng)過濾波器以后波形會有一段延遲，這里截去前邊延遲的無用數(shù)據(jù)，
                                            %雖然實際中并不截去，但是為了處理方便，這里截去無用數(shù)據(jù)
% ----- --- --zaipin generation-----------------------
fc=2/ts;    %載頻是4800hz，說明一個碼元內(nèi)有2個完整的載波波形  
w=2*pi*fc; %--仿真時發(fā)現(xiàn)，載波頻率越高（此時采樣點(diǎn)還是1:length(I0)個，并未提高），誤碼率越低，為什么？
n=(1:length(I0))*(td/2);  %將載波離散化，間隔和數(shù)目應(yīng)與數(shù)據(jù)一致
zaipin1=cos(w*n+pi/4);
%---------Qpsk modulation--------------------------
dataImodulate=I0'.*zaipin1;  %I路BPSK調(diào)制
finalsenddata=dataImodulate;
%-------------figure------------------------------
figure
subplot(2,1,1)
plot(dataI)
title('原始數(shù)據(jù)I差分編碼后')
axis([0 300 -1.5 1.5])
subplot(2,1,2)
plot(dataImodulate)
title('I路bpsk調(diào)制信號')
axis([0 600 -1.5 1.5])
grid on

%------------- Design filter.-----------------------
%因為要濾掉2w頻段，所以截止頻率要<2*w/(2*pi)=9600hz;這里設(shè)計成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); %設(shè)計巴特沃思濾波器
[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;   %數(shù)據(jù)能量

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)); %產(chǎn)生噪聲數(shù)據(jù)
noiseshape=filter(fenzi,fenmu,noise);  %限制噪聲帶寬與數(shù)據(jù)一致
N_power=sum(noiseshape.^2)/length(noiseshape);
bizhi=sqrt(N_hope/N_power);
noise=bizhi*noiseshape;
receivedata=finalsenddata+noise;    %最終發(fā)送數(shù)據(jù)
%----------------------------------------------------
%-------相干解調(diào)----------demodulation-----------------------
%----------------------------------------------------
% --------------Upsample and demodulate--------------
m=(1:length(receivedata))*td/2;
I1=receivedata.*cos(w*m+pi/4);         %相干解調(diào)
%-------------- filter----------------------
I2=filter(fenzi,fenmu,I1);           %低通濾波，濾掉2*w分量
I3=I2;
I3(1:fix(stepnum/3))=[];             %截去延遲的無用數(shù)據(jù)
% 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;                 %取中間點(diǎn)
% receivedataI=zeros(var1);
receivedataI=I3((0:var1-1)*sampletime+var2);  %取中間點(diǎn)進(jìn)行下采樣
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');  %計算誤碼個數(shù)
 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('原始數(shù)據(jù)')
axis([0 64 -1.5 1.5])
subplot(2,1,2)
stem(finaldataout)
 title('解調(diào)數(shù)據(jù)')
 axis([0 64 -1.5 1.5])
 figure
 semilogy(snr_in_dB,pe,'-x')
 title('誤碼率曲線')
 grid