Communications Toolbox 学习记录一:仿真框架
说明：这个系列主要用来记录对matlab的communication toolbox（通信工具箱）的学习，并非教程！
简介
作为一个工科生、特别是通信专业的学生如果不知道matlab 那就太让人匪夷所思了
matlab是美国mathworks公司旗下的一款产品，主要启用于工业仿真（据说最早是从美国一位大学教授为了授课而研发出来的，tql），其产品的核心竞争力就在于软件已经集成了大量的工具箱，就正如现在很多编程语言其强大之处不在于他的语言本身，而在于其强大的工具库，如python在数据处理、机器学习方面的、java在后端开发当中的…… 这个系列博文主要针对matlab在通信仿真当中常用到的Communications Toolbox工具箱 进行学习与记录。
案例一：从QPSK&OFDM仿真开始
通过一个案例讲解一般的仿真基本框架
教程地址
这个教程,展示了怎样构建一个最基本的通信系统，
这个系统采用了QPSK&OFDM调制方式
简单来说，整体上是OFDM，但是每一个OFDM的子载波采用QPSK进行“预先”调制，具体理解的话可以试试看看北邮通信原理第四版P368-371,上面讲了两个OFDM的例子，子载波分别采用BPSK与MQAM调制
参数配置
主要参数
一个通信仿真程序除了最基本的参数之外，就应该开门见山的写配置信息，
% QPSK参数
M = 4; % Modulation alphabet ：调制阶数
k = log2(M); % Bits/symbol： 每符号对应的比特数
% OFDM参数
numSC = 128; % Number of OFDM subcarriers ：OFDM子载波数量
cpLen = 32; % OFDM cyclic prefix length：循环前缀长度
%
maxBitErrors = 100; % Maximum number of bit errors：最大比特错误
maxNumBits = 1e7; % Maximum number of bits transmitted ：最大传输比特数目仿真”对象”
matlab自带了的很多类似的对象函数，目的是减少你重新建模的麻烦，这里更改一些对象的参数就能使用matlab自带的仿真对象
下面是一些初始化对象的操作，对象的具体参数要在下一节调节
QPSK 调制解调器对象
设成直接输入输出二进制 数据的 调制模式
%
qpskMod = comm.QPSKModulator('BitInput',true);
qpskDemod = comm.QPSKDemodulator('BitOutput',true);OFDM 调制解调器
ofdmMod = comm.OFDMModulator('FFTLength',numSC,'CyclicPrefixLength',cpLen);
ofdmDemod = comm.OFDMDemodulator('FFTLength',numSC,'CyclicPrefixLength',cpLen);信道
将NoiseMethodAWGN 通道对象的属性设置为Variance并定义该VarianceSource属性，以便可以从输入端口设置噪声功率。
channel = comm.AWGNChannel('NoiseMethod','Variance', ...
'VarianceSource','Input port');错误率计算器
将ResetInputPort属性设置为true以使错误率计算器能够在模拟过程中被重置。
errorRate = comm.ErrorRate( 'ResetInputPort' ,true);仿真对象的详细参数解释 或者配置
OFDM 调制解调器对象
用info 函数可以查看ofdmMod 对象的输入输出维度
ofdmDims = info(ofdmMod)显示
ofdmDims = struct with fields:
DataInputSize: [117 1]
OutputSize: [160 1]
表明输入一个117*1的二进制数据，得到一组160*1的数据
换句话说，128个子载波中有117个子载波是数据子载波，11个是其他的子载波如空子载波、训练用的?
总共得到160个输出，是因为还有32个循环前缀 ：$128+32=160$
%数据子载波数量 number of data subcarriers
numDC = ofdmDims.DataInputSize(1)显示
numDC = 117
那么一帧OFDM数据（一个OFDM符号）所需要的比特数为：
frameSize = [k*numDC 1];现在理一下帧结构，OFDM帧包括数据帧与循环前缀部分（时域上的划分），在频域上，分为各种子载波，每一个数据子载波“承载”了一个QPSK映射的符号。
信道（高斯白噪声）对象
设置信噪比向量（仿真多个信噪比）
EbNoVec = (0:10) ';
snrVec = EbNoVec + 10 * log10 (k) + 10 * log10 (numDC / numSC);这里不是特别理解？
按照最常用基本的理论
这个实验应该默认sps=1？ k是每个符号对应的比特数 一般$k=R·log_2M$, $R$是编码码率
将numDC / numSC视为编码码率可以好理解一些，就是sps=1 不是很理解，因为sps原本是每个符号的采样点数，OFDM 应该是整体一帧作为一个符号，不可能只采样一个点。如果把sps 理解为：采样增加倍数，这样似乎好理解一些。
后面我再看看相关的书，如果有满意的解释我再修改
错误率计算器对象
初始化BER向量（每一个信噪比一个BER）
初始化错误数据矩阵（暂时不知道，后面看看）
berVec = zeros(length(EbNoVec),3);
errorStats = zeros(1,3);仿真主体
主体循环采用”蒙特卡洛仿真”：再简单BER场景下就是达到一定误码比特数量maxBitErrors，一般设置为100，或者，传输帧数（传输比特数）大于一定数量maxNumBits。
for m = 1:length(EbNoVec) % 每个 信噪比跑一遍
snr = snrVec(m);
while errorStats(2) <= maxBitErrors && errorStats(3) <= maxNumBits
dataIn = randi([0,1],frameSize); % Generate binary data
qpskTx = qpskMod(dataIn); % Apply QPSK modulation
% 由于设置了直接能接受二进制数据所以直接用
txSig = ofdmMod(qpskTx); % Apply OFDM modulation
powerDB = 10*log10(var(txSig)); % Calculate Tx signal power
noiseVar = 10.^(0.1*(powerDB-snr)); % Calculate the noise variance
% 根据信噪比推算噪声能量
rxSig = channel(txSig,noiseVar); % Pass the signal through a noisy channel
% 通过一个信道对象
qpskRx = ofdmDemod(rxSig); % Apply OFDM demodulation
dataOut = qpskDemod(qpskRx); % Apply QPSK demodulation
errorStats = errorRate(dataIn,dataOut,0); % Collect error statistics
end
berVec(m,:) = errorStats; % Save BER data
errorStats = errorRate(dataIn,dataOut,1); % Reset the error rate calculator
% 重置错误率
end用awgn 函数计算 QPSK 的理论误码率（作为参照）
berTheory = berawgn(EbNoVec,'psk',M,'nondiff');画图
figure
semilogy(EbNoVec,berVec(:,1),'*')
hold on
semilogy(EbNoVec,berTheory)
legend('Simulation','Theory','Location','Best')
xlabel('Eb/No (dB)')
ylabel('Bit Error Rate')
grid on
hold off案例二：Parallel Computing ToolBox
案例讲述了如何用Parallel Computing ToolBox 并行计算来加速一般的仿真实验。
教程地址
这个案例没啥太大价值，
其实就是 告诉你 matlab 的for 循环非常慢，可以用parfor代替。
简单来说，原来的for循环采用CPU的一个线程（具体是一个核，还是一个线程我也不知道，就不深入讨论了）来跑，现在，可以用parfor将计算任务分配到多个线程中，这些线程并行的计算来缩减仿真时间。
当然这种方式是有一些限制条件的
简单来说就是
- 1、循环次数多一些的
- 2、循环之间没有迭代关系的
- 3、如果循环用到同一个数组、数组最好循环次数一一对应
具体看看官方文档
案例三：Matlab 仿真AWGN信道下的16QAM 的BER曲线
step 1：使用 MATLAB 检查 16-QAM
基本参数
M = 16; % Modulation order (alphabet size or number of points in signal constellation)
k = log2(M); % Number of bits per symbol
n = 30000; % Number of bits to process
sps = 1; % Number of samples per symbol (oversampling factor)信源
将rng函数设置为其默认状态或任何静态种子值，以便示例产生可重复的结果。然后使用该randi函数生成随机二进制数据。
rng default;
dataIn = randi([0 1],n,1); % Generate vector of binary data检查生成的前40个二进制数据
stem(dataIn(1:40),'filled');
title('Random Bits');
xlabel('Bit Index');
ylabel('Binary Value');
将二进制数据转化为多进制（16进制）数据 （因为我们用到的matlab调制器 对象/函数 需要输入多进制数据），并画图检查
dataSymbolsIn = bit2int (dataIn, k);
figure; % Create new figure window.
stem(dataSymbolsIn(1:10));
title('Random Symbols');
xlabel('Symbol Index');
ylabel('Integer Value');
利用qammod 函数将16进制数据转换为星座点，映射规则默认为格雷映射。
dataMod = qammod(dataSymbolsIn,M,'bin'); % Natural-encoded binary
dataModG = qammod(dataSymbolsIn,M); % Gray-encoded binary噪声
awgn函数加噪声，scatterplot函数查看 星座图（加噪声后的散点图）
EbNo = 10;
snr = EbNo+10*log10(k)-10*log10(sps);
receivedSignal = awgn(dataMod,snr,'measured');
receivedSignalG = awgn(dataModG,snr,'measured');
sPlotFig = scatterplot(receivedSignal,1,0,'g.');
hold on
scatterplot(dataMod,1,0,'k*',sPlotFig)
解调
dataSymbolsOut = qamdemod(receivedSignal,M,'bin'); % Natural-encoded binary data symbols
dataSymbolsOutG = qamdemod(receivedSignalG,M); % Gray-coded binary data symbols变回二进制数据
dataOut = int2bit(dataSymbolsOut,k);
dataOutG = int2bit(dataSymbolsOutG,k);计算误比特
% 自然映射
[numErrors,ber] = biterr(dataIn,dataOut);
fprintf('\nThe binary coding bit error rate is %5.2e, based on %d errors.\n', ...
ber,numErrors)
% 格雷映射
[numErrorsG,berG] = biterr(dataIn,dataOutG);
fprintf('\nThe Gray coding bit error rate is %5.2e, based on %d errors.\n', ...
berG,numErrorsG)
The binary coding bit error rate is 2.27e-03, based on 68 errors.
The Gray coding bit error rate is 1.63e-03, based on 49 errors.
小技巧，显示 映射关系 的星座图
M = 16; % Modulation order
x = (0:15); % Integer input
symbin = qammod(x,M,'bin'); % 16-QAM output (natural-coded binary)
symgray = qammod(x,M,'gray'); % 16-QAM output (Gray-coded)scatterplot(symgray,1,0,'b*'); % 先把点 画出来
for k = 1:M % 对16个点
% 格雷映射
% PS：text 函数text（x轴位置、y轴位置，文本，其他参数（颜色））
text(real(symgray(k)) - 0.0,imag(symgray(k)) + 0.3,...
dec2base(x(k),2,4),'Color',[1 1 0]);
text(real(symgray(k)) - 0.5,imag(symgray(k)) + 0.3,...
num2str(x(k)),'Color',[1 1 0]);
% 自然映射
text(real(symbin(k)) - 0.0,imag(symbin(k)) - 0.3,...
dec2base(x(k),2,4),'Color',[1 0 0]);
text(real(symbin(k)) - 0.5,imag(symbin(k)) - 0.3,...
num2str(x(k)),'Color',[1 0 0]);
end
title('16-QAM Symbol Mapping')
axis([-4 4 -4 4])step2：对 16-QAM 信号使用脉冲整形
上一步以为采样倍率是1 所以一个符号只采样一次 一个符号只用一个数据表示、
为了让仿真会更加贴合实际的效果，还可以加上脉冲成型。
总体来说这个实验就是在 step1 的基础上加上脉冲成型
基本参数
M = 16; % Modulation order
k = log2(M); % Bits per symbol
numBits = k*7.5e4; % Bits to process ；
sps = 4; % Samples per symbol (oversampling factor)根升余弦滤波器 RRC滤波器
filtlen = 10; % Filter length in symbols 滤波器冲激响应持续时间
rolloff = 0.25; % Filter rolloff factor 滚降因子 记得没错的话（1+α）fs 等于带宽（基带）
rrcFilter = rcosdesign(rolloff,filtlen,sps);
fvtool(rrcFilter,'Analysis','Impulse') % 查看冲击相应
其他操作基本上与step1 基本相同
信源
rng default; % Use default random number generator
dataIn = randi([0 1],numBits,1); % Generate vector of binary data
%dataSymbolsIn = bit2int(dataIn,k);
% 由于 bit2int函数r2021b才有
% 我手写实现了一个
temp=reshape(dataIn,k,numBits/k);
table=(2.^((k:-1:1)-1));
dataSymbolsIn=table*temp;
dataMod = qammod(dataSymbolsIn,M);上采样 过成型滤波器
(使用该upfirdn函数以sps倍率的过采样倍数对信号进行上采样，并将信号通过 RRC 滤波器。该upfirdn函数在最后用零填充上采样信号以刷新滤波器。) %不是很懂这个什么叫做 用零填
txFiltSignal = upfirdn(dataMod,rrcFilter,sps,1);
%我运行了一下 dataMod size 1×75000， 按理说sps=4 上采样后 应该为 300000
% 通过滤波器（滤波器10 * 4 +1 =41） 如果是卷积操作，应该大小为 300000+41-1=300040
% 但是结果是一个 300037 不知道为啥？噪声
EbNo = 10；
snr = EbNo + 10*log10(k) - 10*log10(sps)；
rxSignal = awgn (txFiltSignal, snr, 'measured' );接收端
对信号 通过RRC滤波器， 下采样
rxFiltSignal_temp = ...
upfirdn(rxSignal,rrcFilter,1,sps); % Downsample and filter
rxFiltSignal = ...
rxFiltSignal_temp(filtlen + 1:end - filtlen); % Account for delaydataSymbolsOut = qamdemod(rxFiltSignal,M);
% dataOut = int2bit(dataSymbolsOut,k);
temp= de2bi(dataSymbolsOut,k,'left-msb');
dataOut=reshape(temp',[],1);
[numErrors,ber] = biterr(dataIn,dataOut);
fprintf('\nFor an EbNo setting of %3.1f dB, the bit error rate is %5.2e, based on %d errors.\n', ...
EbNo,ber,numErrors)
EbNo = 20;
snr = EbNo + 10*log10(k) - 10*log10(sps);
rxSignal = awgn(txFiltSignal,snr,'measured');
rxFiltSignal_temp = ...
upfirdn(rxSignal,rrcFilter,1,sps); % Downsample and filter
rxFiltSignal = ...
rxFiltSignal(rxFiltSignal_temp + 1:end - filtlen); % Account for delayeyediagram(txFiltSignal(1:2000),sps*2);
eyediagram(rxSignal(1:2000),sps*2);
eyediagram(rxFiltSignal(1:2000),2);scatplot = scatterplot(sqrt(sps)*...
rxSignal(1:sps*5e3),...
sps,0);
hold on;
scatterplot(rxFiltSignal(1:5e3),1,0,'bx',scatplot);
title('Received Signal, Before and After Filtering');
legend('Before Filtering','After Filtering');
axis([-5 5 -5 5]); % Set axis ranges
hold off;M = 16; % Modulation order
k = log2(M); % Bits per symbol
numBits = k*2.5e5; % Total bits to process
sps = 4; % Samples per symbol (oversampling factor)
filtlen = 10; % Filter length in symbols
rolloff = 0.25; % Filter rolloff factorrng default; % Use default random number generator
dataIn = randi([0 1],numBits,1); % Generate vector of binary data
% 卷积码编码
constrlen = [5 4]; % Code constraint length 约束长度
genpoly = [23 35 0; 0 5 13] ; % Generator polynomials
tPoly = poly2trellis(constrlen,genpoly);
codeRate = 2/3;
dataEnc = convenc(dataIn,tPoly);
% 进制转换
%dataSymbolsIn = bit2int(dataEnc,k);
temp=reshape(dataEnc,k,numel(dataEnc)/k);
table=(2.^((k:-1:1)-1));
dataSymbolsIn=table*temp;
% 调制
dataMod = qammod(dataSymbolsIn,M);
txSignal = upfirdn(dataMod,rrcFilter,sps,1);
EbNo = 10;
snr = EbNo+10*log10(k*codeRate)-10*log10(sps);
rxSignal = awgn(txSignal,snr,'measured'); % 下采样
rxFiltSignal = ...
upfirdn(rxSignal,rrcFilter,1,sps); % Downsample and filter
rxFiltSignal = ...
rxFiltSignal(filtlen + 1:end - filtlen); % Account for delay
% 解调
dataSymbolsOut = qamdemod(rxFiltSignal,M);
% 进制转换
% dataOut = int2bit(dataSymbolsOut,k);
temp= de2bi(dataSymbolsOut,k,'left-msb');
codedDataOut =reshape(temp',[],1);
%因为是硬判决所以只用输入 比特数据就好，不用LLR数据
traceBack = 16; % Decoding traceback length
% wtf 什么意思？？？？
numCodeWords = ...
floor(length(codedDataOut)*2/3); % Number of complete codewords
dataOut = ...
vitdec(codedDataOut(1:numCodeWords*3/2), ...
tPoly,traceBack,'cont','hard'); % Decode datadecDelay = 2*traceBack; % Decoder delay, in bits 按比特来算的解码器时延
[numErrors,ber] = ...
biterr(dataIn(1:end - decDelay),dataOut(decDelay + 1:end));
fprintf('\nThe bit error rate is %5.2e, based on %d errors.\n', ...
ber,numErrors)总结
哦吼吼 写（抄了）好多了，这一节就先到这里，后面的再继续写
到这里，已经对一个通信仿真系统基本流程已经了解差不多了
其实里面还有很多不懂得地方，后面搞懂了以后再写写，修改修改（尽量不咕）
评论
