qpsk.m

基于MATLAB的QPSK信道仿真。

%  the simulation of QPSK modulation by simulink
% 1.基带信号波形及功率谱:figure(1)
% 2.已调制信号的波形及功率谱:figure(2)
% 3.经过加性白高斯噪声（AWGN）信道后的信号波形及功率谱:figure(3)
% 4.经过Rayleigh信道前信号的星座图:figure(4)
% 5.经过Rayleigh信道后信号的星座图:figure(5)
% 6.经过Rayleigh衰落信道下的误码率BER图:figure(6)
%     2007.11.08----2007.12.23修改原来图象有问题的现象
clear all;
num=100;                                  %设置原始信息序列的长度
sample=100;                                 %采样点数
num_T=1/sample;                           %采样间隔
N=sample*num;                               %总的采样点数
msg_orig = randint(1,num);                  %生成num个随机序列

data1=2*msg_orig-1;                       %单/双极性变换

ray_I=raylrnd(0.8,1,num/2);                %以0.8为参数的瑞利衰落
ray_Q=raylrnd(0.8,1,num/2);                %以0.8为参数的瑞利衰落

for i=1:num
    data1_orig(sample*(i-1)+1:sample*i)=data1(i); %对data1进行采样
end
for i=1:2000
    data(i)=data1_orig(i);
end
%--------------------------------------------画基带信号及频谱图-----------------------

figure(1);
subplot(2,1,2);                              %在一个窗口里画两个图
pwelch(data1_orig,[],[],[],10);              %画基带信号的频谱图
subplot(2,1,1);
x1=0+0.01:0.01:100;
plot(x1,data1_orig);                          %画基带信号的图型
axis([0 100 -1.2 1.2]);                       %画基带信号图时规定两坐标轴的范围
title('基带信号波形及功率谱');

%----------------------------------将基带信号分成两路 ，并对第I路信号抽样-------------------------
data11=zeros(1,sample*num/2);              
for i1=1:num/2  
    data11(sample*(i1-1)+1:sample*i1)=data1(2*i1-1);%I路信号是基带信号的第奇数个码元，并对它抽样
end

%----------------------------------将基带信号分成两路，并对第Q路信号抽样-------------------------
data12=zeros(1,sample*num/2);              
for i2=1:num/2  
    data12(sample*(i2-1)+1:sample*i2)=data1(2*i2);%Q路信号是基带信号的第偶数个码元，并对它抽样
end

%------------------------------------------对余弦载波抽样------------------------------
for j1=1:N/2                      
    a1(j1)=cos((j1-1)*2*pi/100);                 %每个余弦波每个周期采50个点
end

%------------------------------------------对正弦载波抽样------------------------------
for j2=1:N/2                      
    a2(j2)=sin((j2-1)*2*pi/100);
end

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

%--------------------------------------给两路输出信号叠加上瑞利衰落--------------------------
for i=1:num/2
    data21_ray(i)=data21(i)*ray_I(i);
    data22_ray(i)=data22(i)*ray_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
%-----------------------------------------画调制后的信号波形和频谱图---------------
figure(2); 
subplot(2,1,2);
pwelch(data2_modul,[],[],[],10);       %画调制后的信号波形的功率谱图
subplot(2,1,1);
x2=0+0.01:0.01:10;
plot(x2,data2_modul_t);                 %画调制后的信号波形
title('已调信号波形及功率谱');

data3=awgn(data2_modul,1,'measured');           %Add white Gaussian noise to a signal 经过AWGN信道
for i=1:1000
    data3_t(i)=data3(i);
end
%-----------------------------------画经过高斯信道后的信号波形和功率谱图----------------------
figure(3);  
subplot(2,1,2);                        
pwelch(data3,[],[],[],10);                        %画经过高斯信道后的信号波形的功率谱图
subplot(2,1,1);
plot(x2,data3_t);                                 %画经过高斯信道后的信号波形
title('经过加性白高斯噪声（AWGN）信道后的信号波形及功率谱');

%-----------------------------------画经过Raleigh信道后的信号波形和功率谱图----------------------
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; 

% 信源产生 :生成取值在0和1之间的服从均匀分布的随机变量;信源输出00，10，11，01的概率均为1/4 
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                                       %在Raleigh信道下的星座图
    
    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),'gp');
hold on;
end;
plot(s,'ks');
title('经过Rayleigh信道前信号的星座图');
grid on;
hold off;



for i=1:N
    m=raylrnd(0.7);                                % 生成服从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;
figure(5);
for i=1:1000
plot(rr(i),'gp');
hold on;
end;
plot(s,'ks');
title('经过Rayleigh信道后信号的星座图');
grid on;
hold off;


%------------------------------Rayleigh信道，BER------------------------------
clear;
SNRindB1=0:1:8;                                     % 设置SNR范围
SNRindB2=0:0.1:8;
for ii=1:length(SNRindB1),
counter=0;                                          % 计数器置0
numofsymbolerror=0;
numofbiterror=0;

E=1;                                                % 码元能量
snr=10^(SNRindB1(ii)/10);                           % 信噪比
sgma=sqrt(E/snr)/2;                                 % 噪声功率
s00=[1 0]; s01=[0 -1]; s11=[-1 0]; s10=[0 1];       % 信号映射
while (numofbiterror<100)
% -----------------信源产生----------------------------------------------
    temp=rand;                                      % 生成取值在0和1之间的服从均匀分布的随机变量
    if (temp<0.25),                                 % 信源输出00，10，11，01的概率均为1/4
        dsource1=0; dsource2=0;
    elseif (temp<0.5),                          
        dsource1=0; dsource2=1;
    elseif (temp<0.75),                         
        dsource1=1; dsource2=0;
    else                                        
        dsource1=1; dsource2=1;
    end;
% ------------------------差错检测及误码率计算----------------------
    mm=raylrnd(0.7);                                % 生成服从Rayleigh分布的随机变量
    n=sgma*randn(1,2);                              % 生成两个服从N(0,sgma)的随机变量
    if ((dsource1==0) & (dsource2==0)), 
        r=mm*s00+n;                                 % 加入噪声
    elseif ((dsource1==0) & (dsource2==1)), 
        r=mm*s01+n;
    elseif ((dsource1==1) & (dsource2==0)), 
        r=mm*s10+n;
    else 
        r=mm*s11+n;
    end;
    
    %----------------------- Rayleigh判决-----------------------
    c00=dot(r,s00); 
    c01=dot(r,s01); 
    c10=dot(r,s10); 
    c11=dot(r,s11);
    c_max=max([c00,c01,c10,c11]);
    if (c00==c_max), decis1=0; decis2=0;
    elseif (c01==c_max), decis1=0; decis2=1;
    elseif (c10==c_max), decis1=1; decis2=0;
    else decis1=1; decis2=1;
    end;
    % ---------------计算误码率及误符号率----------------------------
    symbolerror=0;
    if (decis1~=dsource1), numofbiterror=numofbiterror+1; symbolerror=1; 
    end;
    if (decis2~=dsource2), numofbiterror=numofbiterror+1; symbolerror=1; 
    end;
    if (symbolerror==1), numofsymbolerror=numofsymbolerror+1; 
    end;
counter=counter+1;
end;
ps=numofsymbolerror/(counter); 
pb=numofbiterror/(2*counter); 
    smld_bit_err_prb(ii)=pb;
    smld_symbol_err_prb(ii)=ps;
    disp([ps,pb]);
end;

%---------------------------------- 理论值-----------------------------
for i=1:length(SNRindB2),   
    SNR=10^(SNRindB2(i)/10);                        % 信噪比
    theo_err_prb(i)=0.5*erfc(sqrt(SNR));            % AWGN信道下的误码率理论值
    theo_err_ray(i)=0.5*(1-1/sqrt(1+1/SNR));        % Rayleigh衰落信道下的误码率理论值
end;
% ---画出高斯和锐利信道下的BER，误码率和理论值的对比图------------
figure(6)                                        
semilogy(SNRindB1,smld_bit_err_prb,'b-s');
hold on;
semilogy(SNRindB1,smld_symbol_err_prb,'b-o');
hold on;
semilogy(SNRindB2,theo_err_prb,'b');
hold on;
semilogy(SNRindB2,theo_err_ray,'r');
grid on;
title('经过Rayleigh衰落信道下的误码率BER图');
xlabel('Eb/No');
ylabel('Bit erro Rate');
legend('Rayleigh误比特率','Rayleigh误码率','AWGN理论误比特率','RAY理论误比特率');
hold off;