【Pluto SDR】基于OFDM的数字通信系统

引言

我是一名通信工程的一名大三学生。这是，第一次写小文章分享呢。以前包括现在也一直钦佩CSDN博客上一些大佬们的分享，所以这次拿着自己的期末报告来稍稍分享一下。如果以下分享有什么问题，大家可以在评论区指点哈~

在之前的实验中，我们基于16QAM数字调制，搭建了数字带通传输通信系统，实现了数字调制、消除码间串扰的根升余弦滤波器、帧同步、相位同步、信道估计等功能。

但是数字带通传输通信系统所占频段较大，且频带利用率低。针对这些问题，当下现实系统中，有新型数字带通调制可对其性能改善，例如正交2FSK系统、高斯最小频移键控系统、OFDM 通信系统。

接下来，我们将介绍OFDM通信系统的基本理论以及相应的 Pluto Radio SDR 实现。

OFDM 基本概念

OFDM (Orthogonal Frequency Division Multiplexing)即正交频分复用技术，实际上OFDM是MCM(Multi Carrier Modulation)，多载波调制的一种。通过频分复用实现高速串行数据的并行传输, 它具有较好的抗多径衰弱的能力，能够支持多用户接入。

OFDM 技术由MCM（Multi-Carrier Modulation，多载波调制）发展而来。OFDM 技术是多载波传输方案的实现方式之一，它的调制和解调是分别基于IFFT和FFT来实现的，是实现复杂度最低、应用最广的一种多载波传输方案。

OFDM 是一种多载波调制方案，它的基本原理是将传输信道分解为若干个正交的子信道，将需要传输的高速数据信号转换为并行的低速数据流，并调制到每个正交信道的子载波上进行传输，叠加传输的若干正交信号在接收端通过一定方式进行解调分离。

为了使N路子信道信号在接收时能够完全分离，需要满足在每个码元持续时间 内，任意两个载波相互正交。三角函数系$1,sint,cost,...,sinnt,cosnt,...$满足任意两个不同函数之间相互正交（在$[-\pi ,\pi ]$内两个不同函数乘积积分为0），利用这一性质，可以得到彼此相互正交的载波。

从时域上来看，发送端信号在空间中的叠加可以写作：
$$f(t)=\Sigma A_{i}sin(2\pi f_{i}t)+\Sigma B_{i}cos(2\pi f_{i}t)$$
写做复数形式：
$$f(t)=\Sigma B_i{e}^{j2\pi f_{i}t}$$
从复数形式上可以看出，当对时间t进行离散化，OFDM 其实就是对各路子载波上的传输信号 进行了求傅里叶反变换的操作，OFDM 在进行调制时也是利用 IFFT 模块实现的。

从频域上来看OFDM。调制时，时域上是一个正弦载波和一个码元波形(门函数)的乘积，频域上体现的是对门函数频谱的搬移(门函数的频谱为sinc函数), OFDM 为了保证载波的正交性，最小子载频间隔 $\Delta f=1/T_s$

图示频谱为理想状态，即码元波形为边缘垂直的门函数，假设有N路子信道，每路子载波采用M进制调制，则OFDM的带宽：
$$B=\frac{N+1}{2}\frac{2}{T_s}=\frac{N+1}{T_s}$$
信息传输速率：
$$R_b=R_{B}lo{g}_{2}M=\frac{N}{T_s}lo{g}_{2}M$$
故频带利用率：
$$B=\frac{2N}{T_s}$$
分析单载波，为了得到相同的码元传输速率$N/T_s$，所以每个码元的持续时间变为$T_s/N$，带宽变为：

故频带利用率：
$${\eta }_b=\frac{R_b}{B}=\frac{T_s}{2N}\frac{N}{T_s}lo{g}_{2}M=\frac{1}{2}lo{g}_{2}M$$
当N较大时，OFDM 与单载波相比，频带利用率大约变为两倍。

所以，我们应该如何实现 OFDM 呢？

从时域上来看，OFDM 实在发射端发送N路子载波，这些信号在空间中叠加，在接收端接收到叠加的信号，将叠加信号分别对各个子载波做相干解调，就可以得到各路子载波信号，实现起来较为复杂。各路子载波的叠加公式的物理意义可以看作对不同载波上的信号 做傅里叶反变换
$$f(t)=\Sigma B_i{e}^{j2\pi f_{i}t}$$
对$f(t)$用$f_s=N/T_s$频率进行采样：
$$f(k{T}_s/N)={\Sigma }_{i=0}^{N-1}{B}_i{e}^{j\frac{2\pi it}{N}},(0\le k\le N-1)$$
归一化后即为DFT运算：
$$f(k{T}_s/N)=\frac{1}{\sqrt{N}}{\Sigma }_{i=0}^{N-1}{B}_i{e}^{j\frac{2\pi it}{N}},(0\le k\le N-1)$$
因次OFDM可以利用IFFT实现，发射端的IFFT模块对各个载波叠加后的波形进行了计算，从而无需发送大量载波；而接收端的FFT模块可以直接计算得到各路信号，而无需对每一个载波进行相干解调。

OFDM 通信系统理论建模

首先，将我们要传输的信号（例如0~95）进行预处理，将信号二进制化。

将比特流进行QPSK映射，将每两位比特为一个字符，分别映射实部和虚部。其中，“1”表示为数值的1，“0”表示为数值的-1。

将复数信号进行串并变换，并加载子载波。其中，子载波分为空子载波、数据子载波和导频子载波。

其中，这边的导频子载波波段前三个为1，后面一个为-1，主要的目的是为了信道估计，已达到信道均衡。

为了便于接收端帧同步和相位补偿，在信号前加入同步字，增强其抗干扰能力。

将装载好的波段信号分别搬移到相应的波段中。根据前面的理论，可以将信号搬移搬移到不同的子载波理解成傅里叶逆变换，即将频域数据转换为时域波形。由于发送数据为数字信号，即进行IFFT。

为了避免多径干扰和实现时间预估计与频率同步，添加循环前缀（Cyclic Prefix, CP)。

将信号并串变换回序列，定义Pluto Radio传输的中心频率，将其搬移到高频信道进行传输。由此发送端任务完成。

与此类似，接收端可以理解成对发送端进来了一系列的逆过程。不过，由于传输过程之间有延时、信道有干扰、接收器有增益控制等的影响，所以在恢复过程也并非完全的逆过程。

首先，可以利用同步字和循环前缀进行帧同步。利用同步字中的发送前后的相关度不变形，可以设计帧同步的方案，以实现目的。

利用获取点位的同步字，探索频率偏移情况，进行频偏估计与频偏补偿。

利用转载的导频子载波，可进行信道估计，实现信道补偿。

最后，将信号去除CP，FFT 恢复频域信号，获取数据子载波，QPSK 逆映射，二进制转十进制，即获得恢复的信号。

基于SDR的OFDM通信系统实现

我们将利用计算机上的MATLAB软件，对Pluto Radio的收发；通过MATLAB对信号发送前信号数字化处理，调制处理，对接收信号进行时间同步、相位同步，根据接收信号的导频字段进行信道估计与信道均衡等操作。

首先，将待传输信号转换为二进制序列，

将二进制序列装载到子载波框架中，并注入相应的空白子载波与导频子载波，

从上图可以发现，这里的子载波转载后，星座图中出现了0、-1、1，符合转载情况。

在序列头部加上两个同步字，其中第一个同步字在时域有一定特性：实数部分偶对称，虚数部分奇对称——这一特性可以为后续的接收端时间同步、帧同步提供信息。

根据上图，64点的同步字段，实部呈现偶对称特性，而虚部呈现奇对称特性；同时实虚部序列分布，在仅考虑相对幅度变化情况下，虚部是实部的相移90度，形成正交关系。

由于该序列的特性，当其前32点对后32点进行相关函数模长运算，则可获得相关度较高的点。

将序列进行FFT运算后，添加循环前缀CP，本次实验取扩展CP，即前缀长度为符号长度的1/4。

我们一般将FFT变化出来的序列（这里是64点FFT），再加上循环前缀CP，称为一个OFDM符号。让我们看一下待发送的一个OFDM符号的情况以及整个待发送实信号波形。

接下来，我们利用Pluto Radio接收到的信号，对原信号进行恢复。

首先，利用同步字和循环前缀CP，对信号每两个相邻的32点序列进行相关性分析，可绘制接收信号的相关性波形图。
根据下图，发现相关性曲线出现间隔性高峰。将其中一个相关性尖峰的时间轴放大，发现尖峰并非真的如此，而是一个“高原”。

导致这种现象发生的主要原因就是循环前缀CP与同步字的特性。当我们增加16点循环CP，原理上会多出现16个连续的高相关度的采样时刻（加上原先同步字的特性所形成的一个点，即有17点），又因为同步字中第一个点的模极小，所以当偏移一个点时，对相关度影响较小。这也导致出现先前分析的连续18点“高原现象”[1]。

由于，以上操作只是为了帧同步，所以“高原”向后偏移现象是可以容忍的。所以我们只需要保证同步字的第32或第64的模要尽可能大即可。

利用以上相关度“高原”分析法，我们获得一帧数据。

除此之外，同步字还可以用于基于数据辅助的频偏估计，利用相邻32点之间的相角偏移，可以对信号进行频偏补偿——主要的方法就是利用已经估计出来的相角偏移，匹配对应子载波的载波序列，可以根据频率偏移知识对信号进行频偏补偿。

由此我们对信号进行了帧同步、频偏补偿等操作。

将该信号进行FFT，获得频域信号，绘制其星座图。

此时的信号非常凌乱（排除直流分量区域的散点分布），存在大量的错乱现象。

根据我们使用的梳状导频子载波，可进行信道估计，实现信道均衡。将均衡后的星座图展示如图

此时，大量的数据点相对聚集在了四个区域，除了直流分量处与(1,0)、(-1,0)。

（正常情况下，在帧同步前会进行时间同步，由于本次实验未使用过采样，所以就没有进行帧同步；通常在频偏补偿、信道估计后也需要进行相位补偿，不过在实验过程中频偏补偿、信道估计后星座图有较好的恢复，所以就没有再进行）

不过，在虚轴正半轴上依然存在一些零星散点分布，这种情况可能会造成后续的判决器的误判，所以也是后续代码改进的方向。

但是，在这里，我们继续进行实验，暂时保留这种不完美。去除空白子载波和导频子载波，获得所谓的恢复二进制信号。绘制该二进制信号的星座图如图。

简单的观察，我们很难发现有任何误码情况。只能简单分析，有零星散点未能完全的被恢复出来，可能存在误码情况。

为了进一步探究恢复情况，我们对以上二进制序列与原发送序列的二进制流进行相减取模比较。

如图，发现在768个二进制序列中粗出现了一个点的码间串扰现象。对应上图去除载波的二进制信号星座图的涂红圆点。

由此可计算该基于SDR的OFDM通信系统的传输误码率大约为0.13%。

总结与分析

本次实验我们实现了一个基于SDR的OFDM通信系统。该系统利用载波频率正交性，极大的提高了频带利用率；同时，该通信系统，包含了接收端的帧同步、频偏估计及其频偏补偿、信道估计及其信道补偿等技术。

不过，本通信系统尚存在一定的问题，未能在通信过程中实现无误码传输。

本次实验，我对OFDM通信从理论到实践有了更加深刻的理解，也为未来的数字通信的进一步学习奠定了基础。

代码

以下是发送端的MATLAB代码

clc;clear;close all;
data =  0:95; % 待传输数据（96个十进制数）
% 生成比特流
bits_stream = dec2bin(data,8)-48;  % 每行对应一个字节
bits_stream = bits_stream'; % 每列对应一个字节
% QPSK映射
bits_stream = bits_stream*2-1;%将二进制中0变为-1,1变为+1
qpsk_symbol = (1i*bits_stream(1:2:end)+bits_stream(2:2:end))/sqrt(2);% MSB
%绘制星座图
figure(1)
hold on,grid on
scatter(real(qpsk_symbol),imag(qpsk_symbol));
xlabel("Real");ylabel("Imag");title("Binary Sequences");
% 串并变换
qpsk_symbol_s2p = reshape(qpsk_symbol,48,8);
% OFDM符号数据子载波分配（共8个OFDM数据符号/帧）
ofdm_DateCarrier_allocator = zeros(64,8);
ofdm_DateCarrier_allocator(33+[-21 -7 7],:) = 1; % pilot
ofdm_DateCarrier_allocator(33+21,:) = -1; % pilot
ofdm_DateCarrier_allocator(33+(-26:-22),:) = qpsk_symbol_s2p(1:5,:);
ofdm_DateCarrier_allocator(33+(-20:-8),:)   = qpsk_symbol_s2p(6:18,:);
ofdm_DateCarrier_allocator(33+(-6:-1),:)     = qpsk_symbol_s2p(19:24,:);
ofdm_DateCarrier_allocator(33+(1:6),:)       = qpsk_symbol_s2p(25:30,:);
ofdm_DateCarrier_allocator(33+(8:20),:)     = qpsk_symbol_s2p(31:43,:);
ofdm_DateCarrier_allocator(33+(22:26),:)   = qpsk_symbol_s2p(44:48,:);
%绘制星座图
figure(2);hold on,grid on
scatter(real(ofdm_DateCarrier_allocator),imag(ofdm_DateCarrier_allocator));
xlabel("Real");ylabel("Imag");title("Binary Sequences with SubCarrier ");

load sync_word_data.mat
sync_word = [sync_word1;sync_word2];
ofdm_freq = [sync_word',ofdm_DateCarrier_allocator]; % 加同步字-频域
ofdm_time = ifft(ifftshift(ofdm_freq,1),64,1); % IFFT
ofdm_time_addCP = [ofdm_time(49:64,:);ofdm_time]; % 加循环前缀
ofdm_time_addCP = ofdm_time_addCP(:); % 并串变换

% 比较第k+1个OFDM符号
k = 3; 
figure(3);
a1 = 64*fftshift(fft(ofdm_time_addCP(80*k+17:80*k+80)));
hold on;stem(real(a1),'r*-.'); stem(imag(a1),'bo-.')
xlabel('Time')
ylabel('Amplitude')
title("One of OFDM Symbol of Transmition Signal")
% 比较整个时域信号
figure(4);hold on
plot(64*real(ofdm_time_addCP),'r-*')
xlabel('Time')
ylabel('Amplitude')
title("Whole Transmition Signal")


%%
%Pluto 发送
Tdata = ofdm_time_addCP;
tx = sdrtx('Pluto','CenterFrequency',2.5e+09,...
    'BasebandSampleRate',1e6,...
    'Gain',0);
rx = sdrrx('Pluto','CenterFrequency',2.5e+09,...
    'OutputDataType','double',...
    'BasebandSampleRate',1e6,...
    'SamplesPerFrame',1e5,...
    'GainSource','Manual',...
    'Gain',40);                                
%Pluto反复发送数据
tx.transmitRepeat(Tdata);
%接收数据
for k=1:3
    rx();
end
Rdata = rx();
release(tx);
% save OFDM_RX_Pluto.mat Rdata

以下是接收端的MATLAB代码

close all;clear;clc;
load sync_word_data.mat;
load OFDM_RX_Pluto.mat;
Rx_Sig = Rdata(end-2000:end);
N = 0;%画图编号

%%
%----------------------接收信号---------------------
N=N+1;figure(N);hold  on
plot(real(Rx_Sig),'r');
plot(imag(Rx_Sig),'b');
legend("实部","虚部");
title("接收信号")
%%
% -------------------帧同步----------------------------
%同步码转为时域
sync_word1_time = ifft(ifftshift(sync_word1)); 
sync_word2_time= ifft(ifftshift(sync_word2,1));  
%时域同步字+CP
sync_word1_time_cp = [sync_word1_time(end-15:end),sync_word1_time];
sync_word2_time_cp = [sync_word2_time(end-15:end),sync_word2_time];
N=N+1;figure(N);subplot(2,1,1);hold  on;grid on
plot(real(sync_word1_time_cp(17:end)),'b');
set(gca, 'XTickMode', 'manual', 'XTick', [1,17,33,49,64]);
title("real of sync word1 time with circuit prefix")
subplot(2,1,2);hold  on;grid on
plot(imag(sync_word1_time_cp(17:end)),'r');
set(gca, 'XTickMode', 'manual', 'XTick', [1,17,33,49,64]);
title("imag of sync word1 time with circuit prefix")
% N=N+1;figure(N);hold  on
% plot(imag(sync_word2_time_cp),'b');
% title("real of sync word2 time with circuit prefix")
%%
%时域同步字1+CP 的对称特性——出现plateau,即17个连续的最大值
%发现第18个点也较高，在后续使用中可以会混用，所以认为是17个
%原因是只有一个点没有匹配，其他的都配对了
sdata = [zeros(1,100),sync_word1_time_cp,zeros(1,100)];
for d = 1:length(sdata)-63
data_corr(d) = sum(conj(sdata(d:d+31)).*sdata(d+32:d+63));
end
N=N+1;figure(N);hold  on
plot(abs(data_corr),'b');
scatter([101,117],abs(data_corr([101,117])),"filled")
title("Correlation of Synchronic Words 1 with Circuit Prefix")

%%
%接下来开始继续处理接收信号
%利用同步字1及其CP来时间同步
for d = 1:length(Rx_Sig)-63
data_corr(d) = sum(conj(Rx_Sig(d:d+31)).*Rx_Sig(d+32:d+63));
end
N=N+1;figure(N);hold  on
plot(abs(data_corr),'b');
title("Correlation of Received Signal")
index = find(abs(data_corr) > max(abs(data_corr))*9/10);
start_data = index(1)%这里就是数据的起始点
scatter(start_data,abs(data_corr(start_data)),"filled");
Rx_Sig_frame = Rx_Sig(start_data:start_data+799);%读取一帧数据

%%
% --------------------频偏估计-------------------------------
x_ofdm1 = Rx_Sig_frame(17:80);
phase  = angle(conj(x_ofdm1(1:32)).*x_ofdm1(33:64));
delta_f = angle(sum(conj(x_ofdm1(1:32)).*x_ofdm1(33:64)));
if delta_f > pi/2
    delta_f = delta_f - pi;
    elseif delta_f <-pi/2
        delta_f = delta_f + pi;
    else
end
for k = 1:length(phase)
        if phase(k) > pi/2
        phase(k) = phase(k) - pi;
        elseif phase(k) <-pi/2
            phase(k) = phase(k) + pi;
        else
        end
end
N=N+1;figure(N);hold  on
plot(phase*180/pi,'b');
title("freq offset estimation")

% --------------------频偏补偿-------------------------------
% Rx_Sig_frame = Rx_Sig_frame.*exp(-2j*delta_f*(0:799)'/64); % 频偏补偿(为什么是一点点补偿
Rx_Sig_frame = Rx_Sig_frame.*exp(-j*delta_f/2);
%%
% FFT
N=N+1;figure(N);grid on%;hold  on;
for d = 0
x_frame_s2p = reshape(Rx_Sig_frame,80,10);
x_frame_s2p(:,1:2) = []; % 去掉同步字共2个OFDM符号
z_frame = fftshift(fft(x_frame_s2p(17-d:80-d,:),64,1),1);

%星座图
N=N+1;figure(N);hold  on;grid on
scatter(real(z_frame),imag(z_frame));
title("Constellation Diagram");
pause(1);

end
%%
% --------------------信道估计-------------------------------
data_pilot = [1;1;1;-1];
h_est = ones(64,8);
index_pilot = [-21,-7,7,21];
h_est_pilot = z_frame(33 + index_pilot, :) ./ (data_pilot*ones(1,8));
h_est(((-26:26)+33),:) = interp1(index_pilot,h_est_pilot,(-26:26),'linear','extrap');
 
% channel equalization
z_frame = z_frame./h_est;

%星座图
N=N+1;figure(N);grid on;hold  on;
scatter(real(z_frame),imag(z_frame));
title("Constellation Diagram");

%%
% de-mapping
z_frame(33+[-32:-27,-21,-7,0,7,21,27:31],:) = []; % 去掉非数据子载波
data_complex = z_frame(:).';
%%

pause(1);
data_bin = [imag(data_complex)>0;real(data_complex)>0];
c = dec2bin(0:95,8)-48;c = c.';c = c(:);
N=N+1;figure(N);hold  on;
plot(c-data_bin(:),'ro-');hold on;
xlabel("Time");ylabel("Error")
title("Error Conditions")
ylim([-1.1 1.1])
BER = 1-mean(data_bin(:) == c)
%星座图
N=N+1;figure(N);hold  on;grid on
scatter(real(data_complex),imag(data_complex));
error_index = find(c-data_bin(:))
scatter(real(data_complex(error_index)),imag(data_complex(error_index)),'r','filled')
title("Constellation Diagram");

以上代码中会使用到“sync_word_data.mat”文件，需要的uu可以从我的博客的资源库里自提。

有什么疑惑，也可以来询问我，互相讨论切磋~

参考文献

[1] Timothy M. Schmidl,Donald C. Cox.Robust Frequency and Timing Synchronization for OFDM[J].IEEE TRANSACTIONS ON COMMUNICATIONS,VOL.45,NO.12,DECEMBER 1997.