qammodu.m

16QAM调制解调原码

% 16QAM 调制

clear
clc
close all

% data_in = [ 0 0 0 0 0 0 0 1 0 0 1 0 0 0 1 1 0 1 0 0 0 1 0 1 0 1 1 0 0 1 1 1 1 0 0 0 1 0 0 1 1 0 1 0 1 0 1 1 1 1 0 0 1 1 0 1 1 1 1 0 1 1 1 1];
data_in = round(rand(1,5000));
L_data_in=length(data_in);
if rem(L_data_in,4)~=0
    msgbox('输入数据长度应该是4的倍数','提示','error');
    return;
end

QAM_table1=[0 1; -1 2;1 2 ; 0 3;1 0; 2 -1;2 1 ;3 0; -1 0; -2 1; -2 -1 ;-3 0; 0 -1; 1 -2;-1 -2;0 -3];            % 16QAM 映射表
QPSK_table_1=[0 1; 1 0;-1 0; 0 -1];                                                                             % QPSK

for i=1:16
    QAM_table(i)=(QAM_table1(i,1)+QAM_table1(i,2)*j)/sqrt(5);
end

for i=1:4
    QPSK_table(i)=(QPSK_table_1(i,1)+QPSK_table_1(i,2)*j);
end
inp = reshape(data_in,4,L_data_in/4);                               %将data_in矩阵内容转换为4bit并行数据
mod_symbols = QAM_table([8 4 2 1] * inp + 1);                       % 计算出映射后的值
dsource = [8 4 2 1] * inp + 1;                                      % a number between 1 and 16
L_dsource=length(dsource);
QAM_table_real=real(QAM_table);
QAM_table_imag=imag(QAM_table);
%==========================================================================

SNR=0:2:20;                                                       % 信噪比 dB
% SNR=10

for ii=1:length(SNR)
%     Symbols_chaNNL = mod_symbols;
	Symbols_chaNNL = awgn(mod_symbols,SNR(ii),'measured');          % 通过AWGN信道 加高斯白噪声
    L_Symbols_chaNNL=length(Symbols_chaNNL);
    Symbols_chaNNL_real=real(Symbols_chaNNL);                       % 实部    
    Symbols_chaNNL_imag=imag(Symbols_chaNNL);                       % 虚部
end;   
    %===============================
    % 第一种解调方法：按照距离来解
    %===============================
	code_errA=0;
    bit_errA=0;
%     tic
	for k=1:L_Symbols_chaNNL
        % 随后为量度计算
        for jj=1:16
            metrics(jj)=(Symbols_chaNNL_real(k)-QAM_table_real(jj))^2+(Symbols_chaNNL_imag(k)-QAM_table_imag(jj))^2;
        end
        [min_metric decisA]=min(metrics);                            % 找出和16个QAM点距离最近的一个 这个点就作为判决点
		if (decisA~=dsource(k)),
			code_errA=code_errA+1;
		end
        Charge_out(k,:)=de2bi((decisA-1),4,'left-msb');              % 解码比特
        de_modulation_word(k) = decisA;                              % 保存判决符号到 de_modulation_word
	end
%     toc
    p_A(ii)=code_errA/L_Symbols_chaNNL;                                % 误符号率
        
	Charge_out_1=reshape(Charge_out',1,L_data_in);
	for i = 1:L_data_in
        if Charge_out_1(i)~= data_in(i)
            bit_errA=bit_errA+1;
		end
	end
	BER_A(ii)=bit_errA/L_data_in;                                       % 误比特率
    
    %===============================
    % 第二种解调方法：对数似然解
    %===============================
    bit_errL = 0;
    code_errL=0;
%     tic
    in_vit=qampsk_demodulate_x_y(Symbols_chaNNL,16,QAM_table);          % 对数似然解 解码QAM调制
%     toc
    
    for i = 1:L_data_in
        if in_vit(i) >= 0
            in_vit(i) = 1;                                
        else 
            in_vit(i) = 0;
        end
        if in_vit(i) ~= data_in(i)
            bit_errL = bit_errL+1;
        end
	end
    Judge_out=reshape(in_vit,4,L_data_in/4)';
    Judge_out_word = (bi2de(Judge_out,'left-msb')+1)'; 
	BER_L(ii)=bit_errL/L_data_in;                                   % 误比特率
	
	for k=1:L_Symbols_chaNNL
        if (Judge_out_word(k)~=dsource(k)),
			code_errL=code_errL+1;
		end
	end
    p_L(ii)=code_errL/L_Symbols_chaNNL;                             % 误符号率
    
    %===============================
    % 第三种解调方法：改进距离
    %===============================
    code_errB=0;
    bit_errB=0;
    Table=reshape(QAM_table,4,4)';
	Table_real=real(Table);
	Table_imag=imag(Table).*(-1);
%     tic
    for k=1:L_Symbols_chaNNL 
        decisB(k) = demoualtion_16QAM(Symbols_chaNNL(k),Table_real,Table_imag);     % 16QAM 解调
        if (decisB(k)~=dsource(k))
			code_errB=code_errB+1;
		end
        Charge_outB(k,:)=de2bi((decisB(k)-1),4,'left-msb');              % 解码比特
    end
%     toc
    p_B(ii)=code_errB/L_Symbols_chaNNL;
    
	Charge_outB_1=reshape(Charge_outB',1,L_data_in);
	for i = 1:L_data_in
        if Charge_outB_1(i)~= data_in(i)
            bit_errB=bit_errB+1;
		end
	end
	BER_B(ii)=bit_errB/L_data_in;                                       % 误比特率
    
end
dsource(1:16);
de_modulation_word(1:16);
decisB(1:16);

% diary Command_out.txt                   % 保存命令窗口的内容到文本 Command out.txt
%     Symbols_chaNNL;
%     mod_symbols;
%     QAM_table;
%     Judge_out_word(end-10:end);
%     de_modulation_word(end-10:end);
%     bit_errL;
% 	code_errA;
%     BER_L;
%     BER_A;
%     BER_B;
% 	p_L;
%     p_A;
%     p_B;
%     in_vit;
% 	inputdata=reshape(data_in,4,length(data_in)/4)';
% 	inputdata(1:20,:)'
% 	Charge_out(1:20,:)'
%     Charge_outB(1:20,:)'
% diary off

%==========================================图形输出========================


figure(1)
plot(QAM_table,'k*');
hold on
plot(QPSK_table,'rs');
xlabel('实部');
ylabel('虚部');
legend('16QAM','QPSK');

% figure(2)
% x=-1.5:0.5:1.5;
% plot(QAM_table,'k*');
% hold on
% plot(x,x,'r')
% % hold on
% % plot(x,-x,'r')
% text(0.8,-1,'第2比特为 1');
% text(-1.2,1,'第2比特为 0');
% xlabel('实部');
% ylabel('虚部');

figure(3)
plot(Symbols_chaNNL,'bx');
title('经过AWGN信道后的符号');
xlabel('Real');
ylabel('Imag');
grid on;
axis([-2,2,-2,2]);

figure(4)
semilogy(SNR,BER_L,'r');
hold on
semilogy(SNR,BER_A);
hold on
semilogy(SNR,BER_B,'k*');
ylabel('BER');
xlabel('SNR(dB)');
legend('似然解','距离判决','改进算法',3);
grid on

figure(5)
semilogy(SNR,p_L,'r');
hold on
semilogy(SNR,p_A);
hold on
semilogy(SNR,p_B,'k*');
ylabel('误码率');
xlabel('SNR(dB)');
legend('似然解','距离判决','改进算法',3);
grid on
%==========================================================================
figure(6)
plot(SNR,BER_A,'go');

% tt=[-3+3*j,-1+3*j,3+3*j,1+3*j,-3+j,-1+j,3+j,1+j,-3-3*j,-1-3*j,3-3*j,1-3*j,-3-j,-1-j,3-j,1-j];
% a=real(tt)
% b=imag(tt)
% figure
% plot(a,b,'k*')
% axis([-5,5,-5,5]);
% xlabel('实部');
% ylabel('虚部');
% title('标准星座图')