基于MATLAB的数字通信系统仿真

clear all;                  
close all;                  
clc;                        
%% 基本参数
M=10;                       % 产生码元数    
L=100;                      % 每码元复制L次,每个码元采样次数
Ts=0.001;  % 每个码元的宽度,即码元的持续时间
Rb=1/Ts;   % 码元速率1Kbps
dt=Ts/L;   % 采样间隔
TotalT=M*Ts;             % 总时间
t=0:dt:TotalT-dt;        % 时间
Fs=1/dt;    % 采样间隔的倒数即采样频率

%% 产生单极性波形
wave=[0 0 1 1 0 0 1 1 0 0];
% wave=randi([0,1],1,M);      % 产生二进制随机码,M为码元个数
fz=ones(1,L);               % 定义复制的次数L,L为每码元的采样点数
x1=wave(fz,:);              % 将原来wave的第一行复制L次，称为L*M的矩阵
jidai=reshape(x1,1,L*M);    % 产生单极性不归零矩形脉冲波形，将刚得到的L*M矩阵，按列重新排列形成1*(L*M)的矩阵

%% 2ASK调制
fc=10000;                   % 载波频率10kHz       
zb=cos(2*pi*fc*t);   % 载波
ask2=jidai.*zb;      % 2ASK的模拟调制 
figure(1);                  
subplot(3,2,1);               
plot(t,jidai,'LineWidth',2);% 绘制基带码元波形
title('基带信号波形');     
xlabel('时间/s');           
ylabel('幅度');             
axis([0,TotalT,-0.1,1.1])   % 坐标范围限制

subplot(3,2,3)              
plot(t,ask2,'LineWidth',1); % 绘制2ASK的波形 
title('2ASK信号波形');
axis([0,TotalT,-1.1,1.1]);
xlabel('时间/s');
ylabel('幅度');
%% 信号经过高斯白噪声信道
tz=awgn(ask2,20);           % 信号ask2中加入白噪声，信噪比为SNR=20dB
                 
subplot(3,2,5);               
plot(t,tz,'LineWidth',1);   % 绘制2ASK信号加入白噪声的波形
axis([0,TotalT,-1.5,1.5]);  
title('通过高斯白噪声信道后的信号');
xlabel('时间/s');          
ylabel('幅度');            
%% 带通滤波器
%计算滤波器指标数
%B=2Rb=2000HZ,t=1000，DFT=1000
fp1 = 9000; fp2 = 11000;
fs1 = 1000; fs2 = 30000;
wp1 = 2*fp1/Fs; wp2 = 2*fp2/Fs; %通带截止频率，数字归一化频率
ws1 = 2*fs1/Fs; ws2 = 2*fs2/Fs; %阻带截止频率，数字归一化频率
Ap = 1; As = 30;
[N, wc] = buttord([wp1, wp2], [ws1, ws2], Ap, As);
[b, a] = butter(N, wc, 'bandpass');
[h1,w1]=freqz(b, a);
tz=filter(b, a, tz);

%% 解调部分
tz=tz.*zb*2;   % 相干解调，乘以相干载波
subplot(3,2,2)               
plot(t,tz,'LineWidth',1)    % 绘制乘以相干载波后的信号
axis([0,TotalT,-0.5,1.5]);  % 设置坐标范围
title("乘以相干载波后的信号")
xlabel('时间/s');
ylabel('幅度');
%% 加噪信号经过滤波器

fp = 6000;     
fs = 10000;    
wp = 2 * fp / Fs;   
ws = 2 * fs / Fs;   
Ap = 1; As = 30;
% 设计低通滤波器
[N, wc] = buttord(wp, ws, Ap, As);
[b, a] = butter(N, wc, 'low');
[h,w]=freqz(b,a);
% 使用滤波器对信号进行滤波
lvbo = filter(b, a, tz);
 
subplot(3,2,4); % 窗口分割成3*1的，当前是第2个子图 
plot(t,lvbo,'LineWidth',2); % 绘制经过低通滤波器后的信号
axis([0,TotalT,-0.1,1.1]);  % 设置坐标范围
title("经过低通滤波器后的信号");
xlabel('时间/s');
ylabel('幅度');            
%% 抽样判决
k=0.5;              % 设置抽样限值
pdst=1*(lvbo>0.5); % 滤波后的向量的每个元素和0.25进行比较，大于0.25为1，否则为0
subplot(3,2,6)                % 窗口分割成2*1的，当前是第3个子图 
plot(t,pdst,'LineWidth',2)  % 画出经过抽样判决后的信号
axis([0,TotalT,-0.1,1.1]); %设置坐标范用
title("经过抽样判决后的信号");
xlabel('时间/s'); 
ylabel('幅度');
%% 绘制频谱
%% 2ASK信号频谱
T=t(end);                % 时间
df=1/T;                  % 频谱分辨率
N=length(ask2);         % 采样长度
f=(-N/2:N/2-1)*df;       % 频率范围
sf=fftshift(abs(fft(ask2)));% 对2ASK信号采用快速傅里叶变换并将0-fs频谱移动到-fs/2-fs/2
figure(2)      
subplot(222);
plot(f,sf,'LineWidth',2);
title("2ASK信号频谱")
xlabel('频率/Hz');
ylabel('幅度');
%% 信源频谱
mf=fftshift(abs(fft(jidai)));%对信源信号采用快速傅里叶变换并移到矩阵中心
subplot(221);              
plot(f,mf,'LineWidth',2);   % 绘制信源频谱波形
title("基带信号频谱");
xlabel('频率/Hz');
ylabel('幅度');
%% 乘以相干载波后的频谱
mmf=fftshift(abs(fft(tz))); % 对相干载波信号采用快速傅里叶变换并移到矩阵中心             
subplot(223);              
plot(f,mmf,'LineWidth',2)  
title("乘以相干载波后的频谱")
xlabel('频率/Hz');
ylabel('幅度');
%% 经过低通滤波后的频谱
dmf=fftshift(abs(fft(lvbo)));%对低通滤波信号采用快速傅里叶变换并移到矩阵中心
subplot(224);             
plot(f,dmf,'LineWidth',2)  
title("经过低通滤波后的频谱");
xlabel('频率/Hz');
ylabel('幅度');     
figure(3)           
 subplot(211);             
plot(w1/pi*Fs/2,abs(h1),'LineWidth',2);% 绘制滤波器的幅频响应
 xlabel('频率/Hz');      
 ylabel('幅度/dB');         
 title('带通滤波器频谱');  
subplot(212)
plot(w*Fs/(2*pi),abs(h),'LineWidth',2); % 绘制滤波器的幅频响应
title('低通滤波器的频谱');  % 标题
xlabel('频率/Hz');
ylabel('幅度/dB');

3.1.1 时域波形分析

首先生成随机的二进制码元wave=randi([0,1],1,M);这里为了方便观察产生固定的二进制码元wave=[0 0 1 1 0 0 1 1 0 0];一个码元的持续时间位0.001s，在一个码元上采了100个点。

wave=[0 0 1 1 0 0 1 1 0 0];

然后设计了一个10000hz的载波，进行2ASK调制。调制信号2ASK信号波形如图3-2所示。

fc=10000;                   % 载波频率10kHz       

zb=cos(2*pi*fc*t);          % 载波

ask2=jidai.*zb;

    

之后为模仿信道噪声，我们人工设计给信号加上高斯白噪声，信噪比为20dB。加噪信号的时域波形如图3-3所示。tz=awgn(ask2,20);% 信号ask2中加入白噪声，信噪比为SNR=20dB.

 

进入信道后的信号受高斯白噪声影响，我们需要设计一个带通滤波器滤除带外噪声，带通滤波器的频谱图如3-4所示。

fp1 = 9000; fp2 = 11000;

fs1 = 1000; fs2 = 30000;

wp1 = 2*fp1/Fs; wp2 = 2*fp2/Fs;   % 通带截止频率，数字归一化频率

ws1 = 2*fs1/Fs; ws2 = 2*fs2/Fs;   % 阻带截止频率，数字归一化频率

Ap = 1; As = 30;

[N, wc] = buttord([wp1, wp2], [ws1, ws2], Ap, As);

[b, a] = butter(N, wc, 'bandpass');

[h1,w1]=freqz(b, a);

tz=filter(b, a, tz);

  

将通过带通滤波器得到的信号进行解调即乘以同步载波，tz=tz.*zb*2;乘以相干载波后的信号时域波形图3-5所示。这里乘以2倍的载波目的是将其幅度与基带信号幅度相同。具体会在频域分析中说明。

     

解调之后需要一步低通滤波，进一步将除所需信号之外的信号滤除。低通滤波器的频响特性图如图3-6所示，滤除后的波形如图3-7所示。

% 设计低通滤波器

fp = 6000;     

fs = 10000;    

wp = 2 * fp / Fs;   

ws = 2 * fs / Fs;   

Ap = 1; As = 30;

[N, wc] = buttord(wp, ws, Ap, As);

[b, a] = butter(N, wc, 'low');

[h,w]=freqz(b,a);

% 使用滤波器对信号进行滤波

lvbo = filter(b, a, tz)

       

最后是抽样判决。因为2ASK为单极性信号，其判决门限定为0.5，滤波后的向量的每个元素和0.5进行比较，大于0.5为1，否则为0。抽样判决后的信号如图3-8所示。

k=0.5;              % 设置抽样限值

pdst=1*(lvbo>0.5); %滤波后的向量的每个元素和0.5进行比较，大于0.5为1，否则为0

  

通过对比初始基带信号与抽样判决后的信号可以发现，二者几乎完全相同，但存在一定时间的延时。其误码率分析将在后面与2FSK等一起进行分析。

3.1.2 频域分析

初始基带信号的频域波形，2ASK信号频谱图如图3-9所示。

由图可以看出，初始信号的带宽为1000Hz，初始信号频谱图在经过2ASK调制后被搬移到了10000Hz处，且幅度变成了原来的1/2。之后进行解调时，信号的频谱再次进行搬移，由于乘的是同步载波，信号的频谱一部分会搬移到原来的位置，一部分会搬移到20000Hz的位置，乘以相干载波后的信号频谱图如图3-10所示。这里乘的是2倍的载波，目的是使其频谱幅度跟初始频谱幅度相同，如此判决门限需要设为0.5。

由乘以相干载波后的信号频谱图可以看到，我们想要恢复出原始信号，需要将多余频率滤除即可，即通过一个低通滤波器，我们便可以得到原始信号频谱。经过低通滤波器的信号频谱如图3-11所示。

通过理论分析以及实际对比发现，频谱图几乎完全相同，可以有效的恢复出初始基带信号波形。

2FSK，2PSK，2DPSK仿真可私聊，这里不过多分析。下面是4种二进制数字通信系统仿真的实际和理论误码率对比图。