clear all
clc
N=100000;
snr=1:2:24;
NUM_error=0;
h1=zeros(2,N/2);
h2=zeros(2,N/2);
data_1=zeros(2,N/2); %定义data_1矩阵
QPSK_re=zeros(1,N/2); %定义调制的实部
QPSK_im=zeros(1,N/2); %定义调制的虚部
f_frame=rand(1,N)>0.5;%产生随机的二进制比特流
count=0;
for SNR=snr
count=count+1;
QPSK_temp=[-sqrt(2)/2 sqrt(2)/2];
QPSK_re=QPSK_temp(f_frame(:,1:2:end)+1);%二进制比特流的奇数列是QPSK的实部
QPSK_im=QPSK_temp(f_frame(:,2:2:end)+1);%二进制比特流的偶数列是QPSK的虚部
QPSK=QPSK_re+j*QPSK_im; %进行调制
data=reshape(QPSK,2,N/4);
%进行分组空时编码
data_1(:,1:2:end)=data; %产生data_1奇数列
data_1(1,2:2:end)=-conj(data(2,:));%产生偶数列的第一行
data_1(2,2:2:end)=conj(data(1,:)); %产生偶数列的第二行
h=randn(2,N/2)/sqrt(2)+j*randn(2,N/2)/sqrt(2); %产生随机矩阵,均值为0,方差为0.5,randn产生均值为0,标准(方)差为1的正态分布
h1=h(:,1:2:end);
h2=h(:,2:2:end);
H1=kron(h1,ones(1,2)); %构造H矩阵
H2=kron(h2,ones(1,2));
H_1=H1;
H_2=H2;
H1(:,2:2:end)=flipud(H1(:,2:2:end)).*kron(ones(1,N/4),[1;-1]);%构造合并时需要的H矩阵
H1(1,:)=conj(H1(1,:));
H2(:,2:2:end)=flipud(H2(:,2:2:end)).*kron(ones(1,N/4),[1;-1]);
H2(1,:)=conj(H2(1,:));
Hnew=[H1;H2];
n=sqrt(0.5/(10^(SNR/10))) * (randn(2,N/2) + j*randn(2,N/2)); %产生噪声,均值为0,标准差为sqrt(0.5/(10^(SNR/10)))的正态分布
r1=sum(H_1.*data_1,1)+n(1,:);% 接收信号
r2=sum(H_2.*data_1,1)+n(2,:);
r_1=kron(reshape(r1,2,N/4),[1 1]);%构造合并时需要的接收信号矩阵
r_2=kron(reshape(r2,2,N/4),[1 1]);
r_1(2,:)=conj(r_1(2,:));
r_2(2,:)=conj(r_2(2,:));
rnew=[r_1;r_2];
s=sum(Hnew.*rnew,1);%合并信号
dh=sqrt(2)/2*[1+j;1-j;-1+j;-1-j];
s0=repmat(s,size(dh));%复制S矩阵成4行
s1=repmat(dh,1,N/2);%复制dh矩阵成N列
d0=s0-s1;
distance=d0.*conj(d0)+((abs(dh)).^2)*(sum((abs(Hnew)).^2)-1);%最大似然判决
[d,p]=min(distance); %取最小的距离,并记下其位置
decoded=dh(p); %输出该位置的数
%解调
QPSK_output_re = real(decoded.'); %取解码之后数据的实部
QPSK_output_im = imag(decoded.'); %取解码之后数据的虚部
output_frame = zeros(1,N);
id = find(QPSK_output_re>0); %找出实部大于0的,idx记下其位置
QPSK_output_re(id) = 1; %则该位置的实部为1
id = find(QPSK_output_re<0);%找出实部小于0的,idx记下其位置
QPSK_output_re(id) = 0; %则该位置的实部为0
id = find(QPSK_output_im>0);
QPSK_output_im(id) = 1;
id = find(QPSK_output_im<0);
QPSK_output_im(id) = 0;
output_frame(1:2:end) = QPSK_output_re;
output_frame(2:2:end) = QPSK_output_im;
NUM_error=sum(f_frame~=output_frame); %计算解码后的错误数据个数
ber(count)=NUM_error/N; %计算误比特率
end
figure
semilogy(snr,ber,'bp-');
grid on;
xlabel('snr');
ylabel('ber');
title('ber for STBC 2t2r with QPSK modulation and hard decision(Rayleigh channel )');
axis([0 25 10^-5 0.5])