%QPSK_Rayleigh fading.m
%Simulation program to realize QPSK transmission system
%(Both the AWGN channel and  Rayleigh fading channel)
%徐亞沖 編寫整理
clear all;
clc;
%******************** 準備階段 ************************************
nd=100;                                  %原始信息序列長度
samp=100;                                 %采樣點數
delta_T=1/samp;                           %采樣間隔
N=samp*nd;                                %總的采樣點數
data_orig=randint(1,nd);                  %序列生成
data1=2*data_orig-1;                       %單/雙極性變換

ry_I=raylrnd(0.7,1,nd/2);                %rayleigh fading parameter I支路
ry_Q=raylrnd(0.7,1,nd/2);                %rayleigh fading parameter Q支路

for i=1:nd
    data2(samp*(i-1)+1:samp*i)=data1(i);
end
for i=1:2000
    data(i)=data2(i);
end

%******************************將信息源分成兩路并對第I路信號抽樣*************************
data11=zeros(1,samp*nd/2);              
for i1=1:nd/2  
    data11(samp*(i1-1)+1:samp*i1)=data1(2*i1-1);
end

%**********************************將信息源分成兩路并對第Q路信號抽樣**********************
data12=zeros(1,samp*nd/2);              
for i2=1:nd/2  
    data12(samp*(i2-1)+1:samp*i2)=data1(2*i2);
end

%*****************************************對余弦載波抽樣***********************************

for j1=1:samp*nd/2                      
    a1(j1)=cos((j1-1)*pi/25);
end

%*********************************************對正弦載波抽樣*******************************

for j2=1:samp*nd/2                      
    a2(j2)=sin((j2-1)*pi/25);
end

%******************************************QPSK調制****************************************
data21=data11.*a1;  %I路用余弦調制                         
data22=data12.*a2;  %Q路用正弦調制

%*****************************************為兩路信號加上瑞利衰落****************************

for i=1:nd/2
    data21_ray(i)=data21(i)*ry_I(i);
    data22_ray(i)=data22(i)*ry_Q(i);
end

data2_modul=data21+data22;
data2_ray=data21_ray+data22_ray;

for i=1:1000
    data2_modul_t(i)=data2_modul(i);
end
%****************************************經過AWGN信道******************************************
data3=awgn(data2_modul,1,'measured');           
for i=1:1000
    data3_t(i)=data3(i);
end
%***************************************畫信源及頻譜圖*****************************************

figure(1);
subplot(2,1,2);
pwelch(data2,[],[],[],10);
title('基帶信號的功率譜');
subplot(2,1,1);
x1=0+0.01:0.01:100;
plot(x1,data2);
title('基帶信號波形');
axis([0 100 -1.2 1.2]);

%************************************QPSK調制后的信號波形和頻譜圖********************************
figure(2);                                      
subplot(2,1,1);
x2=0+0.01:0.01:10;
plot(x2,data2_modul_t);
title('QPSK調制后波形');
subplot(2,1,2);
pwelch(data2_modul,[],[],[],10);
title('QPSK調制后功率譜');

 %**********************************經過高斯信道后的信號波形和頻譜圖******************************                                 
figure(3);                                     
subplot(2,1,1);
plot(x2,data3_t);
title('經過AWGN信道后的輸出波形');
subplot(2,1,2);
pwelch(data3,[],[],[],10);
title('經過AWGN信道后的信號功率譜');

%************************************************星座圖********************************************


clear;                                          
N=1000;
s00=[1 0]; s01=[0 -1]; s11=[-1 0]; s10=[0 1];   % signal mapping
SNRindB1=8;
E=1;                                            % energy per symbol
snr=10^(SNRindB1/10); 
sgma=sqrt(E/snr)/2; 
for i=1:N
    temp=rand;                                  % a uniform random variable between 0 and 1
    if (temp<0.25),                             % with probability 1/4, source output is "00"
        dsource1(i)=0; dsource2(i)=0;
        s(i)=complex(s00(1),s00(2));
    elseif (temp<0.5),                          % with probability 1/4, source output is "01"
        dsource1(i)=0; dsource2(i)=1;
        s(i)=complex(s01(1),s01(2));
    elseif (temp<0.75),                         % with probability 1/4, source output is "10"
        dsource1(i)=1; dsource2(i)=0;
        s(i)=complex(s10(1),s10(2));
    else                                        % with probability 1/4, source output is "11"
        dsource1(i)=1; dsource2(i)=1;
        s(i)=complex(s11(1),s11(2));
    end;
end;


for i=1:N                                      
    n=sgma*randn(1,2);                          
    if ((dsource1(i)==0) & (dsource2(i)==0)), 
        r=s00+n;
    elseif ((dsource1(i)==0) & (dsource2(i)==1)), 
        r=s01+n;
    elseif ((dsource1(i)==1) & (dsource2(i)==0)), 
        r=s10+n;
    else 
        r=s11+n;
    end;
    rr(i)=complex(r(1),r(2));
end;
%在高斯信道下的星座圖
figure(4);
for i=1:1000
plot(rr(i),'b*');
hold on;
end;
plot(s,'gd');    
grid on;
hold off;
title('經過高斯信道后的星座圖');
%--------------------------------------------------------------------------

for i=1:N
    m=raylrnd(0.8);                                 % 生成服從Rayleigh分布的隨機變量
    n=sgma*randn(1,2);                              % 生成兩個服從N(0,sgma)的隨機變量
    if ((dsource1(i)==0) & (dsource2(i)==0)), 
        r=m*s00+n;
    elseif ((dsource1(i)==0) & (dsource2(i)==1)), 
        r=m*s01+n;
    elseif ((dsource1(i)==1) & (dsource2(i)==0)), 
        r=m*s10+n;
    else 
        r=m*s11+n;
    end;
    rr(i)=complex(r(1),r(2));
end;
%在Raleigh信道下的星座圖
figure(5);
for i=1:1000
plot(rr(i),'b*');
hold on;
end;
plot(s,'gd');
grid on;
hold off;    
title('經過瑞利信道后的星座圖');
%************************************ 程序結束 ********************************************