qpsk.m

qpsk的实现及其星座图希望大家有好的资源可以互相共享

T=100;        %原始信息序列长度
samp=100;     %采样点数
N=samp*T;
m=randint(1,T);     %生成序列
data1=2*m-1;                  %单/双极性变换

for i=1:T
    data1_orig(samp*(i-1)+1:samp*i)=data1(i);  %信息序列扩展
end

% 信源序列图及其频谱图
figure(1);
subplot(2,1,2);
pwelch(data1_orig,[],[],[],10);%data1_orig功率普密度
title('信源序列功率普密度');
subplot(2,1,1);
x1=0+0.01:0.01:100;
plot(x1,data1_orig);axis([0 100 -1.2 1.2]);
title('信源序列');

% 将信源分成两路
data11=zeros(1,N/2);            
for i1=1:T/2  
    data11(samp*(i1-1)+1:samp*i1)=data1(2*i1-1);
end
data12=zeros(1,N/2);              
for i2=1:T/2  
    data12(samp*(i2-1)+1:samp*i2)=data1(2*i2);
end

%  对载波抽样
for j1=1:N/2                      
    a1(j1)=cos((j1-1)*pi/20);
end
for j2=1:N/2                      
    a2(j2)=sin((j2-1)*pi/20);
end

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

data3_modul=data21+data22;                 %经过理想信道的QPSK信号
data3_gauss=awgn(data3_modul,10,'measured');  %经过AWGN信道后的QPSK信号

% 经过rayleigh信道
rly=raylrnd(0.6,1,T/2);                %rayleigh fading parameter

for i=1:T/2
   onl(i)=sqrt(rly(i)*rly(i)+rly(i)*rly(i));
   data21_ray(samp*(i-1)+1:samp*i)=data21(samp*(i-1)+1:samp*i)*onl(i);      %只考虑幅度瑞利衰落
   data21_ray_p(samp*(i-1)+1:samp*i)=data21(samp*(i-1)+1:samp*i)*rly(i)-data22(samp*(i-1)+1:samp*i)*rly(i);     %考虑幅度相位瑞利衰落
   data22_ray(samp*(i-1)+1:samp*i)=data22(i)*onl(i);                    %只考虑幅度瑞利衰落
   data22_ray_p(samp*(i-1)+1:samp*i)=data21(samp*(i-1)+1:samp*i)*rly(i)+data22(samp*(i-1)+1:samp*i)*rly(i);   %考虑幅度相位瑞利衰落
end

data3_ray=data21_ray+data22_ray;                   %生成只考虑幅度瑞利QPSK信号
data3_ray_p=data21_ray_p+data22_ray_p;        %生成考虑幅度相位瑞利QPSK信号

for i=1:2000
    data3_modul_t(i)=data3_modul(i);
    data3_gauss_t(i)=data3_gauss(i);
    data3_ray_t(i)=data3_ray(i);
    data3_ray_p_t(i)=data3_ray_p(i);
end

%调制后的信号波形和频谱图
figure(2);      
subplot(2,1,2);
pwelch(data3_modul,[],[],[],10);  title('调制后无衰落QPSK信号频谱图');                                    
subplot(2,1,1);
x2=0+0.01:0.01:20;
plot(x2,data3_modul_t);   title('调制后无衰落QPSK信号波形图');

%只经过高斯信道后的信号波形和频谱图
figure(3);                                      
subplot(2,1,2);
pwelch(data3_gauss,[],[],[],10);title('经AWGN信道后的信号频谱图(信噪比10dB)');
subplot(2,1,1);
x2=0+0.01:0.01:20;
plot(x2,data3_gauss_t);title('经AWGN信道后的信号波形图(信噪比10dB)'); 

%经过rayleigh衰落的信号波形和频谱图
figure(4);                                     
subplot(2,1,2);
pwelch(data3_ray_t,[],[],[],10);title('QPSK瑞利衰落幅度信号频谱图（参数0.6）');
subplot(2,1,1);
x2=0+0.01:0.01:20;
plot(x2,data3_ray_t);title('QPSK瑞利衰落幅度信号波形图（参数0.6）');
figure(5);                                     
subplot(2,1,2);
pwelch(data3_ray_p_t,[],[],[],10);title('QPSK瑞利衰落相位信号频谱图（参数0.6）');
subplot(2,1,1);
x2=0+0.01:0.01:20;
plot(x2,data3_ray_p_t);title('QPSK瑞利衰落相位信号波形图（参数0.6）');

%经过rayleigh信道的星座图
N=100;
s00=[1 0]; s01=[0 -1]; s11=[-1 0]; s10=[0 1];   % signal mapping
SNRindB1=30;
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:100             %在高斯信道下的星座图
    
    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(6);
for i=1:N
plot(rr(i),'b*');title('高斯干扰下的星座图（30dB）');
hold on;
end;
plot(s,'rd');
grid on;
hold off;

for i=1:N
    m=raylrnd(0.6);                                % 生成服从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(7);
for i=1:N
plot(rr(i),'b*');title('瑞利干扰下的星座图（参数0.6）');
hold on;
end;