16QAM调制解调和误码率

课程设计作业

16QAM调制解调

clear all;
close all;
%%
%16QAM信号调制解调
%通过瑞利/高斯信道
%理论和仿真误码率

%%
%系统参数设计
T_start=0;%开始时间
T_stop=40;%截止时间
T=T_stop-T_start;%仿真持续时间

rs=10;%传输速率
NumBits=rs*T;%传输bit数

fc=20;%载波频率

L=100;%码元长度
N_sample=NumBits*L;% 采样点数
T_sample=T/N_sample;%采样间隔
f_sample=1/T_sample;%采样频率
f_res=f_sample/N_sample;%频率分辨率

%载波
n=1:N_sample/4;
c1=cos(2*pi*fc*n*T_sample);
c2=sin(2*pi*fc*n*T_sample);

%%
%16QAM信号的调制

%随机产生原始01信息比特
a=rand(1,NumBits)>0.5;

%串并变换
A=reshape(a,2,NumBits/2);
a_I=A(1,:);
a_Q=A(2,:);

%2-4电平变换，00对应-3，01对应-1，10对应1，11对应3
a_4I=zeros(1,NumBits/4);
for i=1:NumBits/4
    g=a_I(2*i-1)*2+a_I(2*i);
    if g==0
        a_4I(i)=-3;
    elseif g==1
        a_4I(i)=-1;
    elseif g==2
        a_4I(i)=1;
    else
        a_4I(i)=3;
    end
end
a_4Q=zeros(1,NumBits/4);
for i=1:NumBits/4
    g=a_Q(2*i-1)*2+a_Q(2*i);
    if g==0
        a_4Q(i)=-3;
    elseif g==1
        a_4Q(i)=-1;
    elseif g==2
        a_4Q(i)=1;
    else
        a_4Q(i)=3;
    end
end
%
%将信息比特扩展成方波，得到两路基带信号
e_I=[];
for i=1:NumBits/4
    B(1:L)=a_4I(i);
    e_I=[e_I,B];
end
e_Q=[];
for i=1:NumBits/4
    B(1:L)=a_4Q(i);
    e_Q=[e_Q,B];
end


%基带信号与载波相乘
sI=e_I.*c1;%I支路调制信号
sQ=e_Q.*c2;%Q支路调制信号
s=sI-sQ;%16QAM调制信号


%%
%瑞利信道
SNR=15;%输入信噪比
snr=10^(SNR/10);
N0=1/2/snr;%计算噪声功率
N0_dB=10*log10(N0);%将噪声功率转换为dBW
h=raylrnd(1/sqrt(2),1,N_sample/4);%产生瑞利信号
n_i=wgn(1,N_sample/4,N0_dB);%产生高斯噪声
y_n=s.*h+n_i;%16QAM调制信号通过瑞利信道
y_n0=s+n_i;%16QAM信号通过高斯信道

%%
%16QAM信号的解调

%通过带通滤波器，去除噪声
[b1,a1]=butter(4,[0.02,0.06]);%计算带通滤波器的H(z)系数
y=filtfilt(b1,a1,y_n);%对信号y_i进行滤波，得到信号y;使用filtfilt进行滤波几乎没有延时

%与恢复载波相乘
x1_I=y.*c1;
x1_Q=-y.*c2;

%通过低通滤波器，分离出两路信号
[b2,a2]=butter(4,0.02,'low');%计算低通滤波器H(z)系数
x2_I=filtfilt(b2,a2,x1_I);%对信号x1进行滤波，得到信号x2_I
x2_Q=filtfilt(b2,a2,x1_Q);%对信号x1进行滤波，得到信号x2_Q

%判决
x_I=fun_panjue16QAM(x2_I);%调用函数，对I支路进行判决
x_Q=fun_panjue16QAM(x2_Q);

%提取信息比特，，得到-3,-1,1,3信息
d_I0=fun_tiqu16QAM(x_I,L);
d_Q0=fun_tiqu16QAM(x_Q,L);

%4-2进制转换,将-3,-1,1,3转换成01
d_I=zeros(1,NumBits/2);
for i=1:NumBits/4
    if d_I0(i)==3
        d_I(2*i-1)=1;
        d_I(2*i)=1;
    elseif d_I0(i)==1
        d_I(2*i-1)=1;
        d_I(2*i)=0;
    elseif d_I0(i)==-1
        d_I(2*i-1)=0;
        d_I(2*i)=1;
    else
        d_I(2*i-1)=0;
        d_I(2*i)=0;
    end
end
d_Q=zeros(1,NumBits/2);
for i=1:NumBits/4
    if d_Q0(i)==3
        d_Q(2*i-1)=1;
        d_Q(2*i)=1;
    elseif d_Q0(i)==1
        d_Q(2*i-1)=1;
        d_Q(2*i)=0;
    elseif d_Q0(i)==-1
        d_Q(2*i-1)=0;
        d_Q(2*i)=1;
    else
        d_Q(2*i-1)=0;
        d_Q(2*i)=0;
    end
end
        
%并串变换
d(1:2:NumBits)=d_I;
d(2:2:NumBits)=d_Q;

%%
%误码率计算   
err=length(find(d~=a));%计算解调信号中错误码元个数
BER_test=err/NumBits;%计算误码率

%%
%图像
figure(1);
subplot(4,2,1:2);stem(a);title('原始01信息序列a');axis([1,20,0,1]);
subplot(4,2,3);stem(a_I);title('I支路01信息序列a_I');axis([1,10,0,1]);
subplot(4,2,4);stem(a_Q);title('Q支路01信息序列a_Q');axis([1,10,0,1]);
subplot(4,2,5);stem(a_4I);grid on;
title('四进制I支路信息序列');axis([1,10,-4,4]);
subplot(4,2,6);stem(a_4Q);grid on;
title('四进制Q支路信息序列');axis([1,10,-4,4]);
subplot(4,2,7);plot(e_I);hold on;plot(sI);grid on;
legend('I支路基带信号e_I','I支路调制信号sI');axis([0,1000,-4,4]);
subplot(4,2,8);plot(e_Q);hold on;plot(sQ);grid on;
legend('Q支路基带信号e_Q','Q支路调制信号sQ');axis([0,1000,-4,4]);

figure(2);
subplot(3,1,1);plot(s);title('16QAM调制信号s');axis([0,2000,-6,6]);
subplot(3,1,2);plot(y_n);title('16QAM调制信号s通过瑞利信道后的信号y_n');axis([0,2000,-6,6]);
subplot(3,1,3);plot(y_n0);title('16QAM调制信号通过高斯信道后的信号y_n0');axis([0,2000,-6,6]);


figure(3);
subplot(3,2,1:2);plot(y);title('y_i通过带通滤波器后的信号y');axis([0,4000,-6,6]);
subplot(3,2,3);plot(x1_I);title('y与恢复载波c1相乘后的信号x1_I');axis([0,2000,-6,6]);
subplot(3,2,4);plot(x1_Q);title('y与恢复载波c2相乘后的信号x1_Q');axis([0,2000,-6,6]);
subplot(3,2,5);plot(x2_I);hold on;plot(e_I);
grid on;legend('x1_I通过低通滤波器后的信号x2_I','原始I支路基带信号e_I');
axis([0,2000,-6,6]);
subplot(3,2,6);plot(x2_Q);hold on;plot(e_Q);
grid on;legend('x1_Q通过低通滤波器后的信号x2_Q','原始Q支路基带信号e_Q');
axis([0,2000,-6,6]);

figure(4);
subplot(4,2,1);plot(x_I);hold on;plot(e_I);
grid on;legend('x2_I经过判决后得到的信号x_I','原始I支路基带信号e_I');
axis([0,2000,-4,4]);
subplot(4,2,2);plot(x_Q);hold on;plot(e_Q);
grid on;legend('x2_Q经过判决后得到的信号x_Q','原始Q支路基带信号e_Q');
axis([0,2000,-4,4]);
subplot(4,2,3);stem(d_I0);grid on;
title('x2_I经过信息比特提取后得到的信息序列d_I0');axis([1,5,-4,4]);
subplot(4,2,4);stem(d_Q0);grid on;
title('x2_Q经过信息比特提取后得到的信息序列d_Q0');axis([1,5,-4,4]);
subplot(4,2,5);stem(d_I);title('d_I0经过4-2进制转换后得到的信息序列d_I');axis([1,10,0,1]);
subplot(4,2,6);stem(d_Q);title('d_I0经过4-2进制转换后得到的信息序列d_Q');axis([1,10,0,1]);
subplot(4,2,7:8);stem(d);hold on;stem(a,'*');
legend('d_I和d_Q并串变换后得到解调01信息序列','原始01信息序列');axis([1,20,0,1.5]);

子程序1

function w=fun_panjue16QAM(w)
%以1,0,-1为基准进行判决
N=length(w);
for i=1:N
    if w(i)>=1
        w(i)=3;
    elseif w(i)<1&&w(i)>=0;
        w(i)=1;
    elseif w(i)<0&&w(i)>-1
        w(i)=-1;
    else
        w(i)=-3;
    end
end

子程序2

function h=fun_tiqu16QAM(x,L)
%提取出信息比特
%x为输入信号，L为信号的码元长度
for i=1:length(x)/L;
    k=x((i-1)*L+1:i*L);
    v1=length(find(k==3));
    v2=length(find(k==1));
    v3=length(find(k==-1));
    v4=length(find(k==-3));    
    if v1>v2&&v1>v3&&v1>v4
        h(i)=3;
    elseif v2>v1&&v2>v3&&v2>v4
        h(i)=1;
    elseif v3>v1&&v3>v2&&v3>v4
        h(i)=-1;
    elseif v4>v1&&v4>v2&&v4>v3
        h(i)=-3;
    else
        h(i)=0;
    end          
end

16QAM误码率

clear all;
close all;
%误码率曲线的绘制
SNR=1:1:20;%信噪比变化范围
snr=10.^(SNR/10);%将信噪比转化成直角坐标
N=1000000;%仿真点数
M=16;%进制数
x=randi([0,1],1,N);%产生随机信号
y=qammod(x,M,pi/M);%调用matlab自带的qam调制函数
for i=1:length(SNR);
    N0=1/2/snr(i);%计算噪声功率
    N0_dB=10*log10(N0);%将噪声功率转换为dBW
    ni=wgn(1,N,N0_dB);%产生高斯噪声
    h=raylrnd(1/sqrt(2),1,N);%产生瑞利信号

    yAn=y+ni;%通过高斯信道
    yA=qamdemod(yAn,M,pi/M);%调用matlab自带的qam解调函数
    bit_A=length(find(x~=yA));%统计错误比特数
    BPSK_s_AWGN(i)=bit_A/N;%计算误码率
    
    yRn=y.*h+ni;%通过瑞利信道
    yR=qamdemod(yRn,M,pi/M);%调用matlab自带的qam解调函数
    bit_R=length(find(x~=yR));%统计错误比特数
    BPSK_s_Ray(i)=bit_R/N;%计算误码率 
end
p1=6*log2(M)/(M-1)*snr;
p2=(1-1/sqrt(M))*erfc(sqrt(p1));
BPSK_t_AWGN=2*p2.*(1-p2/2);%AWGN信道下16QAM理论误码率
BPSK_t_Ray=1/2*(1-sqrt((1/2*p1)./(1+1/2*p1)));%Rayleigh信道下16QAM理论误码率

%绘制图形
figure;
semilogy(SNR,BPSK_s_AWGN,'-k*');hold on;
semilogy(SNR,BPSK_t_AWGN,'-go');hold on;
semilogy(SNR,BPSK_s_Ray,'-b*');hold on
semilogy(SNR,BPSK_t_Ray,'-ro');grid on;
axis([-1,20,10^-4,1]);
legend('AWGN仿真','AWGN理论','瑞利仿真','瑞利理论');
title('16QAM误码性能分析');
xlabel('信噪比（dB）');ylabel('BER');