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;