project1.m

16－QAM基带传输系统的Matlab仿真 实现了主要模块二进制信息源输入

% Project1

% 日期:2006.12.28

close all;clear;clc
%Variable declare
N=10000;                                                                    % 二进制码流比特数
T=0.001;                                                                    % 16QAM符号长度(s)
Fd=1/T;                                                                     % 符号率
Fs=4000;                                                                    % 过采样率(Hz)
r=Fs*T;                                                                     % 过采样因子
alpha=0.25;                                                                 % 滚降因子
Fc=500;                                                                     % 信道截止频率(Hz)
SNR=15;                                                                     % 信道的信噪比

%Bit Source
source=randint(N,1);                                                        % 产生二进制码流

%QAM Coder
send=bi2de(flipud(reshape(source,4,N/4))');
QAM=qammod(send,16);
h1=scatterplot(QAM);
set(h1,'Name','Constellation Diagram');
title('Constellation of 16QAM in source');

%Transmitter
p=rcosine(Fd,Fs,'sqrt',alpha);                                              % 均方根升余弦滤波器

len=length(QAM);
QAMS=zeros(size(1:r*len));
QAMS(1:r:r*len)=QAM;
QAM_x=real(QAMS);
QAM_y=imag(QAMS);
delay=(length(p)-1)/2;
temp1=zeros(size(1:delay));
T_x_a=[QAM_x temp1];
T_y_a=[QAM_y temp1];
T_x_a=filter(p,1,T_x_a);
T_y_a=filter(p,1,T_y_a);
T_x=T_x_a(1+delay:10000+delay);
T_y=T_y_a(1+delay:10000+delay);
figure(2);
hold on;
for i=(length(p)+1)/2:r*2:len-r*2,
    plot(T_x(i:i+r*2));
end;
hold off;
figure(3);
freqz(T_x,1,256);

%Channal
%[B A]=butter(2,2*Fc/Fs);
B=fir1(8,2*Fc/Fs);A=1;
delay=(length(B)-1)/2;
temp1=zeros(size(1:delay));
C_x_a=[T_x temp1];
C_y_a=[T_y temp1];
C_x_a=filter(B,A,C_x_a);
C_y_a=filter(B,A,C_y_a);
C_x=C_x_a(1+delay:10000+delay);
C_y=C_y_a(1+delay:10000+delay);
C_x=awgn(C_x,SNR);
C_y=awgn(C_y,SNR);
figure(4);
freqz(C_x,1,256);

%Equalizer
load c;
freqz(B,A,256);
hold on;
freqz(c,1,256);
hold off;
delay=(length(c)-1)/2;
temp1=zeros(size(1:delay));
E_x_a=[C_x temp1];
E_y_a=[C_y temp1];
E_x_a=filter(c,1,E_x_a);
E_y_a=filter(c,1,E_y_a);
E_x=E_x_a(1+delay:10000+delay);
E_y=E_y_a(1+delay:10000+delay);

%E_x=C_x;
%E_y=C_y;

%Receiver
delay=(length(p)-1)/2;
temp1=zeros(size(1:delay));
R_x_a=[E_x temp1];
R_y_a=[E_y temp1];
R_x_a=filter(p,1,R_x_a);
R_y_a=filter(p,1,R_y_a);
R_x=R_x_a(1+delay:10000+delay);
R_y=R_y_a(1+delay:10000+delay);
figure(5);
hold on;
for i=(length(p)+1)/2:r*2:len-r*2,
    plot(R_x(i:i+r*2));
end;
hold off;

%Sample
DQAM_x=R_x(1:4:10000);
DQAM_y=R_y(1:4:10000);

Error_x=DQAM_x'-real(QAM);
Error_y=DQAM_y'-imag(QAM);

DQAM=DQAM_x+j*DQAM_y;
h2=scatterplot(DQAM);

%Decoder
rcv = qamdemod(DQAM,16);

%Error
[Err_num,Err_ratio]=symerr(rcv,send')