qam_16.m

来自「简单的数字传输系统---主要是对16-QAM数字基带传输系统进行仿真」· M 代码 · 共 203 行

clc;close all;
clear;                                %清除变量

%定义变量
N_bit=10000;                          %二进制码流的比特数
Fs=3000;                              %采样率
Fc=500;                               %信道的截止频率(Hz)
SNR=20;                               %取值范围为：5-20dB
T=1/1000;                             %符号持续时间
Rd=1/T;                               %符号率
alfa=0.25;                            %升余弦滚降因子
R=Fs*T;                               %过采样因子

%%%%%%%%%QAM Coder%%%%%%%%%%
bit_source=randint(N_bit,1);          %产生10000bit的2进制随机序列
send=bi2de(flipud(reshape(bit_source,4,N_bit/4))');
qam=qammod(send,16);                  %产生16-QAM信号
Iqam=real(qam);
Qqam=imag(qam);
N=length(qam);

%16-QAM星座图映射图
figure(1);                                           
plot(Iqam,Qqam,'*');           
axis([-4 4 -4 4]);
title('16-QAM星座图');
xlabel('I');ylabel('Q');
grid on;

qam_tx=zeros(size(1:R*N));
qam_tx(1:R:R*N)=qam;                   %QAM信号补零
Iqam_tx=real(qam_tx);
Qqam_tx=imag(qam_tx);

%%%%%%%%%发送%%%%%%%%%%%
%发送滤波器
Pt=rcosine(Rd,Fs,'sqrt',alfa);         %升余弦成型滤波器
Delay_tx=(length(Pt)-1)/2;             %发送滤波器的延迟
Delay_tx=round(Delay_tx);              %取整
tx_zero=zeros(size(1:Delay_tx));
Iqam_tx=[Iqam_tx tx_zero];
Qqam_tx=[Qqam_tx tx_zero];

Iqam_tx=filter(Pt,1,Iqam_tx);          %信号经过发送滤波器
Qqam_tx=filter(Pt,1,Qqam_tx);          %信号经过发送滤波器

I_tx=Iqam_tx(1+Delay_tx:end);
Q_tx=Qqam_tx(1+Delay_tx:end);
tx=I_tx+j*Q_tx;

%16-QAM通过发送滤波器后的星座图
figure(2);                                             
plot(I_tx,Q_tx,'*');          
title('16-QAM通过发送滤波器后的星座');
xlabel('I');ylabel('Q');
grid on;

%信号通过发送滤波器后的眼图
figure(3);                                              
hold on;
for i=(5*R+1):2*R:R*N-2*R
    plot(I_tx(i:i+2*R));
end
hold off;
title('信号通过发送滤波器后的眼图');
xlabel('Time');ylabel('Amplitude');


%%原始信号和发送滤波器的频谱
figure(4);
f=-Fs/2:Fs/2;                             %要求正负半轴的频谱都有
Fa_qam=abs(freqz(qam,1,f,Fs));
Fa_filter=abs(freqz(Pt,1,f,Fs));
plot(f,20*log10(Fa_qam),f,20*log10(Fa_filter),'r');
title('原始信号和发送滤波器频谱');
xlabel('f');ylabel('A(dB)');
grid on;

%%信号通过发送滤波器后的频谱
fa_tx=abs(freqz(tx,1,f,Fs));
figure(5);
plot(f,20*log10(fa_tx));
title('信号通过滤波器后的频谱');
xlabel('f');ylabel('A(dB)');
grid on;

%%%%%%%%信道%%%%%%%%%%
[b,a]=butter(3,2*Fc/Fs);                   %截止频率为Fc的3阶巴特沃斯滤波器
%figure(6);                                 %信道的频谱图
%freqz(b,a,512);

%%信号经过信道
Delay_ch=(length(b)-1)/2;                  %计算延时
Delay_ch=round(Delay_ch);
ch_zero=zeros(size(1:Delay_ch));
I_tx=[I_tx ch_zero];
Q_tx=[Q_tx ch_zero];

I_ch=filter(b,a,I_tx);                     %信号经过信道滤波器
Q_ch=filter(b,a,Q_tx);                     %信号经过信道滤波器

I_ch=I_ch(1+Delay_ch:end);
Q_ch=Q_ch(1+Delay_ch:end);
%Signal_ch=I_ch+j*Q_ch;

I_ch=awgn(I_ch,SNR);                      %加高斯白噪声
Q_ch=awgn(Q_ch,SNR);                      %加高斯白噪声
ch=I_ch+j*Q_ch;

 %画信号经过信道加白噪声后的频谱
f=-Fs/4:Fs/4;                             %设置频率范围
fa_ch=abs(freqz(ch,1,f,Fs));                                
figure(7);
plot(f,20*log10(fa_ch));                  %采用对数坐标
xlabel('f');ylabel('A(dB)');
title('信号经过信道加白噪声后的频谱');
grid on; ; 

%%%%%%%%Receiver%%%%%%%%%
delay_rx=(length(Pt)-1)/2;
delay_rx=round(delay_rx);
delay_rx_zero=zeros(size(1:delay_rx));
I_ch=[I_ch delay_rx_zero];
Q_ch=[Q_ch delay_rx_zero];
I_rx=filter(Pt,1,I_ch);                   %信号通过接受滤波器
Q_rx=filter(Pt,1,Q_ch);                   %信号通过接受滤波器

I_rx=I_rx(1+delay_rx:end);
Q_rx=Q_rx(1+delay_rx:end);

%画信号通过接受滤波器后的眼图
figure(8);
hold on;
for i=(length(Pt)+1)/2:R*2:N-R*2,
    plot(I_rx(i:i+R*2));
end;
hold off;
title('信号通过接受滤波器后的眼图');
xlabel('Time');ylabel('Amplitude');

%%%%%%%%%采样%%%%%%%%%%%%%
Impulse = [1 zeros(1,100)];               % Create an impulse.
TxImpOut = filter(Pt,1,Impulse);          % Find out the combined impulse response of the Tx-filter, channel, 
                                          %and Rx-filter. Both Rx- and Tx-filters are
ChannelImpOut = filter(b,a,TxImpOut);     % square root raised cosine filters.
RxImpOut = filter(Pt,1,ChannelImpOut);
[Trash,Pos] = max(abs(RxImpOut));         % Find the maximum value and its position, now the position is delay. 
                                %Start sampling of the filtered received signal from
                                 %the position Pos. In the actual sampling, use, of course,
                                          %  your own received signal variable instead of RxImpOut.
RxSampOut_I = I_rx(1:R:end);
RxSampOut_Q = Q_rx(1:R:end);
RxSampOut=RxSampOut_I+j*RxSampOut_Q;

%采样后星座映射图
figure(9);
plot(RxSampOut_I,RxSampOut_Q,'*');          
axis([-4 4 -4 4]);
title('采样后信号星座映射');
xlabel('I');ylabel('Q');
grid on;

%%%%%%%%%%%均衡%%%%%%%%%%%%
beta=0.01;                                 % 步长
c=zeros(size(1:31))';                      %均衡系数的初始化
for i=Pos:960
 rk=flipud( RxSampOut_I(i:i+30)');         %Received signal vector
 Ek(i)=Iqam(i+15)-c'*rk;                   %Error signal
 c=c+beta*Ek(i)*rk;
 end 

delay_c=(length(c)-1)/2;                    %均衡器延迟
delay_c=round(delay_c);
temp_zero=zeros(size(1:delay_c));
RxSampOut_I=[RxSampOut_I temp_zero];
RxSampOut_Q=[RxSampOut_Q temp_zero];

RxSampOut_I=filter(c,1,RxSampOut_I);         %信号通过均衡器
RxSampOut_Q=filter(c,1,RxSampOut_Q);         %信号通过均衡器

RxSampOut_I=RxSampOut_I(1+delay_c:end);
RxSampOut_Q=RxSampOut_Q(1+delay_c:end);
RxSampOut=RxSampOut_I+j*RxSampOut_Q;

%%信号通过均衡器后的星座映射图
figure(10);                                            
plot(RxSampOut_I,RxSampOut_Q,'*');          
axis([-4 4 -4 4]);
title('信号通过均衡器后的星座映射');
xlabel('I');ylabel('Q');
grid on;

%%%%%%%%%%Decoder%%%%%%%%%%%
rcv = qamdemod(RxSampOut,16);

%%%%%误码率和无码数的计算%%%%%%
[Err_num,Err_ratio]=symerr(rcv,send')        %计算误码数和误码率

%close all;