OFDM系统设计 Matlab实现

通信系统集成电路设计
专业课
发布日期:   2021-06-02
更新日期:   2021-09-19
文章字数:   6k
阅读时长:   27 分

OFDM系统设计 Matlab实现

OFDM介绍

在时域上看，如果发送的符号/码元周期T远小于时延扩展Td，那么接收端由于多径时延，接收信号中的一个符号的波形会扩展到其他符号当中， 造成符号间干扰 （InterSymbol Interference，ISI），即码间串扰。从频域上看，符号周期 等效于 ， 即信号带宽或者信号速率远大于相干带宽时， 信号通过无线信道后某些频率成分信号的幅值可以增强（相位相同，幅值叠加），而另外一些频率成分信号的幅值会被削弱（相位相反，幅值相消），引起信号波形的失真，此时就认为发生了频率选择性衰落。

相反，在时域上看，符号周期远大于时延扩展 时，多径时延的影响就很小，没有符号间干扰。从频域上看，信号带宽远小于相干带宽，接收端的信号通过无线信道后各频率分量都受到相同的衰落（由于各频率成分比较接近，经多条路径到达接收端后相位也很接近，各频率成分各自叠加后的幅度也很接近），因而衰落波形不会失真，则认为信号只是经历了平坦衰落。

为了对抗频率选择性衰落，人们采用了正交频分复用（OFDM）技术，该技术将宽带信号分成很多子带，频域上分成很多子载波发送出去，每个子带的信号带宽由于小于相干带宽，从而减少甚至避免了频率选择性衰落。

OFDM: 正交频分复用技术。载波中心处在正交频率上，子载波间隔为1/T（T是符号持续时间），且相互正交，意味着子信道之间的串扰被消除，并且不需要载波间保护带。子载波具有保持其相应时域波形正交性所需的最小频率间隔。

调制和多路复用技术

OFDM系统的设计原理关键如下图。

OFDM原理

采用OFDM将频域选择性衰落转变为平坦衰落， 并通过交织分散深衰落， 再通过编码恢复深衰落，即通过频率分离解决频域选择性衰落问题。

OFDM系统设计

OFDM系统设计的流程如下图。

OFDM系统的发送/接收过程

主要参数

本实验用到的主要参数如下表。

主要参数

训练序列

基于802.11a的OFDM帧结构

OFDM的前导训练序列包含了十个相同的短训练序列和两个相同的长训练序列，中间利用长度为32的保护间隔隔开。前导训练序列是用来做OFDM系统的包检测、时频同步、信道估计和补偿等。下面就短训练序列和长训练序列进行具体介绍。

短训练序列

$t_1\sim t_{10}$代表10个短训练序列，其调制因子为：

乘以$\sqrt{13/6}$是为了将52个子载波中的12个进行能量归一化。这十二个子带为$[-24,-20,-16,-12,-8,-4,4,8,12,16,20,24]$。

每个短训练序列长度为16，则总的段训练时域长度为160。短训练序列的具体作用是进行包检测、频偏(CFO)估计和纠正。程序体现如下。

% 训练定义
% 用于CFO估计的短训练序列 (符号数量NumSymbols = 52)
ShortTrain = sqrt(13/6)*[0 0 1+j 0 0 0 -1-j 0 0 0 1+j 0 0 0 -1-j 0 0 0 -1-j 0 ...
0 0 1+j 0 0 0 0 0 0 -1-j 0 0 0 -1-j 0 0 0 1+j 0 0 0 1+j 0 0 0 1+j 0 0 0 1+j 0 0].';
NumShortTrainBlks = 10;

short_train = tx_freqd_to_timed(ShortTrain);
short_train_blk = short_train(1:16);
short_train_blks = repmat(short_train_blk,NumShortTrainBlks,1);

长训练序列

$T_1$和$T_2$代表两个长训练序列，其调制因子为：

每个长训练序列的时域长度为64，由调制过的52个子载波构成。其中，在长短训练序列间插入了长度为32的保护间隔 ，来防止长、短训练序列间相互干扰。所以加上保护间隔后，长训练序列的总时域长度为160。

​ 综上，前导训练序列的总时域长度为短训练序列和长训练序列的和，即总时域长度为320。

调制与解调

本实验使用4QAM对子载波进行调制与解调，阶数为2。在程序中对应如下。

% 调制
paradata = reshape(CodedSeq,length(CodedSeq)/ml,ml);
ModedSeq = qammod(bi2de(paradata),2^ml)/NormFactor;

% 解调
DemodSeq = de2bi(qamdemod(Data_seq*NormFactor,2^ml),ml);
deint_bits = reshape(DemodSeq,size(DemodSeq,1)*ml,1).';

信道编码与译码

在信道编码方面，采用802.11a的（2,1,7）卷积码对DATA部分进行卷积编码。其约束长度为7，卷积编码的码率为$1/2$，对应编码器如下图。

卷积编码器

信道编码在程序中体现为：

inf_bits = randn(1,NumBitsPerPkt)>0; % 生成信息比特
CodedSeq = convenc(inf_bits,trellis); % 对输入数据进行卷积编码

在信道译码方面，采用的是维特比译码。维特比译码算法是基于卷积码网格图的最大似然译码算法。译码部分在程序中体现为：

% 维特比译码
DecodedBits = vitdec(deint_bits(1:length(CodedSeq)),trellis,tb,'trunc','hard');

插入导频

将导频插入到标号为-21、-7、7、21的子载波上，作用是防止频谱偏移并增强自相关的检测性能。在其它的48个子载波上插入的是有效数据。

交织与解交织

一个QAM中的比特在编码比特序列中不是连续的，即把QAM比特分散开，这个过程就叫交织。解交织即是交织的相反过程。错误的字节在时间和频率上分散，可以更好地进行纠错。交织/解交织可以恢复非突发错误比特。

交织/解交织

映射与解映射

交织之后有映射(Mapping)这个步骤。各子载波的标号分配如下图。

子载波映射标号表

将子载波标号（id）为1~26的符号映射IFFT的1~26作为输入，将-26~-1的符号映射到IFFT的38~63作为输入，其余位置无输入（设为0）。其中数据子载波标号分别1~5、7~19、21~26、27~32、34~46、48~52。导频子载波标号为6、20、33、47，如下图所示。

映射后IFFT的输入输出

循环前缀CP

通过引入循环前缀充当保护间隔，在时域内把OFDM符号的后面部分插入到符号的开始部分，用来消除多径时延扩展引起的码间串扰（ISI）和载波间的干扰（ICI）。引入的循环前缀使得在OFDM符号的持续期间，每个子载波总是存在整数个周期。本次选择的循环前缀长度为16。

​ 通过在发送端加入CP，在接收端去除CP，将多径信道的线性卷积变换为循环卷积，从而简化了频域均衡。

在保护间隔加入CP

OFDM关键技术流程及结果

包检测

对在接收端获取接收的信号进行包检测，通过信号中的短训练序列，检测收到数据包的起始时间点，便于后续处理。

实现原理

包检测

对信号做归一化自相关，通过确定的一个门限值$T_h=0.75$来确定是否检测到包，相关公式如下：

将归一化自相关得到的$m_n$与门限值$T_h$比较，如果大于门限值则表示收到了包，然后便将此时刻作为数据的起始点。

程序实现

Pkt detection.m

%% ************** 准备部分 ********************
clear all; clc;
% 系统参数
gi = 1/4;            % 保护间隔: 最大时延为2e-6s，相当于16个点. 
fftlen = 64;
gilen = gi*fftlen;   % 保护间隔长度（点）

% 训练序列
% 短训练序列
ShortTrain = sqrt(13/6)*[0 0 1+j 0 0 0 -1-j 0 0 0 1+j 0 0 0 -1-j 0 0 0 -1-j 0 ...
0 0 1+j 0 0 0 0 0 0 -1-j 0 0 0 -1-j 0 0 0 1+j 0 0 0 1+j 0 0 0 1+j 0 0 0 1+j 0 0].';
NumShortTrainBlks = 10; % 短训练序列个数

short_train = tx_freqd_to_timed(ShortTrain); % 短训练序列转化到时域
figure
plot(abs(short_train));
short_train_blk = short_train(1:16); % 只保留前16个
short_train_blks = repmat(short_train_blk,NumShortTrainBlks,1); % 把十个短训练序列接在一起

% 长训练序列
LongTrain = [1 1 -1 -1 1 1 -1 1 -1 1 1 1 1 1 1 -1 -1 1 1 -1 1 -1 1 1 1 1 ...
      1 -1 -1 1 1 -1 1 -1 1 -1 -1 -1 -1 -1 1 1 -1 -1 1 -1 1 -1 1 1 1 1].';
NumLongTrainBlks = 2; % 长训练序列个数

long_train = tx_freqd_to_timed(LongTrain); % 长训练序列转换到时域
long_train_syms = [long_train(fftlen-2*gilen+1:fftlen,:); long_train; long_train];
% 将两个长训练序列接在一起，并在前面加上保护间隔，隔开长短训练序列
%% ************** 信道特性 ***************************
h = zeros(gilen,1);
h(1) = 1;
h(5) = 0.5;
h(10) = 0.3; 
start_noise_len = 500;
snr = 20;

%% ************** 发送 ***************************
tx = [short_train_blks; long_train_syms]; % 将所有的训练序列连在一起

%% ************** 经过信道 ***************************
rx_signal = filter(h,1,tx); % 经过信道后接收
noise_var = 1/(10^(snr/10))/2;
len = length(rx_signal); % 接收信号的长度，即总训练序列的长度16*10+(32+64*2)=320
noise = sqrt(noise_var) * (randn(len,1) + j*randn(len,1)); % 设置噪声
% 添加噪声
rx_signal = rx_signal + noise;
start_noise = sqrt(noise_var) * (randn(start_noise_len,1) + j*randn(start_noise_len,1));
rx_signal = [start_noise; rx_signal]; % 添加噪声后的接收信号


%% ************** 接收端 ***************************
search_win = 700; 
D = 16;

% 计算时延自相关
delay_xcorr = rx_signal(1:search_win+2*D).*conj(rx_signal(1*D+1:search_win+3*D));

% 时延自相关的移动平均值
ma_delay_xcorr = abs(filter(ones(1,2*D), 1, delay_xcorr));

% 接收能量的移动平均值
ma_rx_pwr = filter(ones(1,2*D), 1, abs(rx_signal(1*D+1:search_win+3*D)).^2);

% 决策变量
delay_len = length(ma_delay_xcorr);
ma_M = ma_delay_xcorr(1:delay_len)./ma_rx_pwr(1:delay_len);
figure 

% 删除延迟样本
ma_M(1:2*D) = [];
plot(ma_M);

threshold = 0.75;
thres_idx = find(ma_M > threshold); % 与阈值比较，比阈值大说明检测到包

if isempty(thres_idx)
  thres_idx = 1;
else
  thres_idx = thres_idx(1);
end

detected_packet = rx_signal(thres_idx:length(rx_signal));

函数tx_freqd_to_timed.m

function time_syms = tx_freqd_to_timed(mod_ofdm_syms) 
num_symbols = size(mod_ofdm_syms,2);

UsedSubcIdx = [7:32 34:59];
resample_patt=[33:64 1:32];

syms_into_ifft = zeros(64, num_symbols);
syms_into_ifft(UsedSubcIdx,:) = mod_ofdm_syms;

syms_into_ifft(resample_patt,:) = syms_into_ifft;

% Convert to time domain
time_syms = sqrt(64)*ifft(sqrt(64/52)*syms_into_ifft);

仿真结果

仿真结果

由上图可以看出，当检测到数据包时，对应取样表示的序列号处产生上升沿。

频率同步

OFDM信号在通过信道之后有可能会产生载波频偏，通过数据中的短训练序列来进行相关运算提取出频偏的估计，并进行校正达到同步的目的。

实现过程

频率同步_过程

首先进行相关运算：

载波频偏(CFO)估计：

载波频偏(CFO)补偿：

程序实现

Freq Sync.m

%% ************** 准备部分 ********************
clear all; clc;
% 系统参数
fs = 20e6;           % 采用频率
gi = 1/4;            % 保护间隔: 最大时延为2e-6s，相当于16个点. 
fftlen = 64;
gilen = gi*fftlen;   % 保护间隔长度（点）

% 训练序列
% 短训练序列
ShortTrain = sqrt(13/6)*[0 0 1+j 0 0 0 -1-j 0 0 0 1+j 0 0 0 -1-j 0 0 0 -1-j 0 ...
0 0 1+j 0 0 0 0 0 0 -1-j 0 0 0 -1-j 0 0 0 1+j 0 0 0 1+j 0 0 0 1+j 0 0 0 1+j 0 0].';
NumShortTrainBlks = 10; % 短训练序列个数

short_train = tx_freqd_to_timed(ShortTrain); % 短训练序列转化到时域
figure
plot(abs(short_train));
short_train_blk = short_train(1:16); % 只保留前16个
short_train_blks = repmat(short_train_blk,NumShortTrainBlks,1); % 把十个短训练序列接在一起

% 长训练序列
LongTrain = [1 1 -1 -1 1 1 -1 1 -1 1 1 1 1 1 1 -1 -1 1 1 -1 1 -1 1 1 1 1 ...
      1 -1 -1 1 1 -1 1 -1 1 -1 -1 -1 -1 -1 1 1 -1 -1 1 -1 1 -1 1 1 1 1].';
NumLongTrainBlks = 2; % 长训练序列个数

long_train = tx_freqd_to_timed(LongTrain); % 长训练序列转换到时域
long_train_syms = [long_train(fftlen-2*gilen+1:fftlen,:); long_train; long_train];
% 将两个长训练序列接在一起，并在前面加上保护间隔，隔开长短训练序列
%% ************** 信道 ***************************
h = zeros(gilen,1);
h(1) = 1;
h(5) = 0.5;
h(10) = 0.3;
cfo = 0.1*fs/fftlen % 设置载波频偏
%% ************** 循环启动 ***************************
snr = 2:2:20;
mse = zeros(1,length(snr));
pkt_num = 2000; % 包的数量
for snr_idx = 1:length(snr) 
    est_err = zeros(1,pkt_num); 
    for pkt_idx = 1:pkt_num % 对每个包进行检测
        %% 发送
        tx = [short_train_blks; long_train_syms]; % 长短训练序列连接在一起

        %% 信道
        rx_signal = filter(h,1,tx);
        noise_var = 1/(10^(snr(snr_idx)/10))/2;
        len = length(rx_signal);
        noise = sqrt(noise_var) * (randn(len,1) + j*randn(len,1));
        % 加入噪声
        rx_signal = rx_signal + noise;
        % 加入CFO(后面进行频率同步)
        total_length = length(rx_signal);
        t = [0:total_length-1]/fs;
        phase_shift = exp(j*2*pi*cfo*t).';
        rx_signal = rx_signal.*phase_shift;

        %% 接收端
        % 对于开始时的生样本（以及实际系统中的同步错误）
        pkt_det_offset = 30;
        % 平均长度
        rlen = length(short_train_blks) - pkt_det_offset;
        % 短训练序列符号周期
        D = 16;

        phase = rx_signal(pkt_det_offset:pkt_det_offset+rlen-D).* ...
          conj(rx_signal(pkt_det_offset+D:pkt_det_offset+rlen));

        % 加上所有估计
        phase = sum(phase);

        freq_est = -angle(phase) / (2*D*pi/fs);
        est_err(pkt_idx) = (freq_est - cfo)/cfo; % 估计错误

    end
    mse(snr_idx) = mean(abs(est_err).^2); %MSE: 均方误差，Mean Squared Error
end
semilogy(snr,mse);
xlabel('SNR/dB');
ylabel('MSE');

函数tx_freqd_to_timed.m

function time_syms = tx_freqd_to_timed(mod_ofdm_syms) 
num_symbols = size(mod_ofdm_syms,2);

UsedSubcIdx = [7:32 34:59];
resample_patt=[33:64 1:32];

syms_into_ifft = zeros(64, num_symbols);
syms_into_ifft(UsedSubcIdx,:) = mod_ofdm_syms;

syms_into_ifft(resample_patt,:) = syms_into_ifft;

% Convert to time domain
time_syms = sqrt(64)*ifft(sqrt(64/52)*syms_into_ifft);

仿真结果

频率同步_SNR-MSE曲线

由上图可以看出，随着信噪比（SNR）的增加，均方误差（MSE）逐渐下降，估计的结果越来越精，与预想结果一致。

时间精同步

时间精同步的目的是为了找到OFDM符号的起始位并定位长训练序列，删除端训练序列，补偿CFO。

实现过程

时间精同步

接收到的序列与长训练序列共轭相乘，找到最大值的位置即为长训练序列的起始位置。

程序实现

Time Sync.m

%% ************** 准备部分 ********************
clear all; clc;
% 系统参数
gi = 1/4;                   % 保护间隔: 最大时延为2e-6s，相当于16个点. 
fftlen = 64;
gilen = gi*fftlen;           % 保护间隔长度(点)

% 训练序列
% 短训练序列
ShortTrain = sqrt(13/6)*[0 0 1+j 0 0 0 -1-j 0 0 0 1+j 0 0 0 -1-j 0 0 0 -1-j 0 ...
0 0 1+j 0 0 0 0 0 0 -1-j 0 0 0 -1-j 0 0 0 1+j 0 0 0 1+j 0 0 0 1+j 0 0 0 1+j 0 0].';
NumShortTrainBlks = 10; % 短训练序列个数

short_train = tx_freqd_to_timed(ShortTrain); % 短训练序列转化到时域
figure
plot(abs(short_train));
short_train_blk = short_train(1:16); % 只保留前16个
short_train_blks = repmat(short_train_blk,NumShortTrainBlks,1); % 把十个短训练序列接在一起

% 长训练序列
LongTrain = [1 1 -1 -1 1 1 -1 1 -1 1 1 1 1 1 1 -1 -1 1 1 -1 1 -1 1 1 1 1 ...
      1 -1 -1 1 1 -1 1 -1 1 -1 -1 -1 -1 -1 1 1 -1 -1 1 -1 1 -1 1 1 1 1].';
NumLongTrainBlks = 2; % 长训练序列个数

long_train = tx_freqd_to_timed(LongTrain); % 长训练序列转换到时域
long_train_syms = [long_train(fftlen-2*gilen+1:fftlen,:); long_train; long_train];
% 将两个长训练序列接在一起，并在前面加上保护间隔，隔开长短训练序列
%% ************** 信道 ***************************
h = zeros(gilen,1);
h(1) = 1;
h(5) = 0.5;
h(10) = 0.3;
%% ************** 循环启动 ***************************
snr = -10:5:0;
mse = zeros(1,length(snr));
pkt_num = 1000; % 包的数量
ideal_start = 193;
for snr_idx = 1:length(snr)
    snr_idx
    est_err = zeros(1,pkt_num);
    for pkt_idx = 1:pkt_num % 对每个包进行检测
        %% 发送
        tx = [short_train_blks; long_train_syms]; % 长短训练序列连接在一起

        %% 信道
        rx_signal = filter(h,1,tx);
        noise_var = 1/(10^(snr(snr_idx)/10))/2;
        len = length(rx_signal);
        noise = sqrt(noise_var) * (randn(len,1) + j*randn(len,1));
        % 添加噪声
        rx_signal = rx_signal + noise;

        %% 接收端
        % 时域搜索窗口大小
        start_search=150;
        end_search=200;
        time_corr_long = zeros(1,end_search-start_search+1);

        for idx=start_search:end_search
            time_corr_long(idx-start_search+1) = sum((rx_signal(idx:idx+63).*conj(long_train)));
        end

        [max_corr,long_search_idx] = max(abs(time_corr_long));

        fine_time_est = start_search-1 + long_search_idx;
        err_est(pkt_idx) = fine_time_est - ideal_start;
    end
    mse(snr_idx) = mean(abs(err_est).^2);
end
semilogy(snr,mse,'-b.');hold on;
xlabel('SNR/dB');
ylabel('MSE');

函数tx_freqd_to_timed.m

function time_syms = tx_freqd_to_timed(mod_ofdm_syms) 
num_symbols = size(mod_ofdm_syms,2);

UsedSubcIdx = [7:32 34:59];
resample_patt=[33:64 1:32];

syms_into_ifft = zeros(64, num_symbols);
syms_into_ifft(UsedSubcIdx,:) = mod_ofdm_syms;

syms_into_ifft(resample_patt,:) = syms_into_ifft;

% Convert to time domain
time_syms = sqrt(64)*ifft(sqrt(64/52)*syms_into_ifft);

仿真结果

时间精同步_SNR/MSE曲线

在信噪比较低时，随着信噪比（SNR）的增加，均方误差（MSE）不断降低。在信噪比较高时，MSE将一直为0。

信道估计

实现原理

采用迫零检测法，利用长训练序列进行信道估计。

程序实现

ChaEst.m

%% ************** 准备部分 ********************
clear all; clc;
% 系统参数
fs =20e6;
ml = 2;                      % 调制模式: 2--4QAM; 4--16QAM; 6--64QAM
NormFactor = sqrt(2/3*(ml.^2-1));
gi = 1/4;                   % 保护间隔: 最大时延为2e-6s，相当于16个点. 
fftlen = 64;
gilen = gi*fftlen;           % 保护间隔长度 (点)
blocklen = fftlen + gilen;   % CP块的总长度

% 参数定义

DataSubcPatt = [1:5 7:19 21:26 27:32 34:46 48:52];
PilotSubcPatt = [6 20 33 47];
UsedSubcIdx = [7:32 34:59];

% 长训练序列
LongTrain = [1 1 -1 -1 1 1 -1 1 -1 1 1 1 1 1 1 -1 -1 1 1 -1 1 -1 1 1 1 1 ...
      1 -1 -1 1 1 -1 1 -1 1 -1 -1 -1 -1 -1 1 1 -1 -1 1 -1 1 -1 1 1 1 1].';
NumLongTrainBlks = 2; % 长训练序列个数

long_train = tx_freqd_to_timed(LongTrain); % 长训练序列转换到时域
long_train_syms = [long_train(fftlen-2*gilen+1:fftlen,:); long_train; long_train];
% 将两个长训练序列接在一起，并在前面加上保护间隔，隔开长短训练序列
%% ************** 信道 ***************************
h = zeros(gilen,1);
h(1) = 1;
h(5) = 0.5;
h(10) = 0.3;
reorder = [33:64 1:32];
channel = fft(h, 64);
channel(reorder) = channel;
channel = channel(UsedSubcIdx);

%% ************** 循环启动 ***************************
snr = 0:5:20;
mse = zeros(1,length(snr));
pkt_num = 1000; % 包的数量
for snr_idx = 1:length(snr) 
    err_est = zeros(1, pkt_num);
    for pkt_idx = 1:pkt_num % 对每个包进行检测
        tx = long_train_syms; % 短序列在时间精同步的时候已经被去除

        rx_signal = filter(h,1,tx);
        noise_var = 1/(10^(snr(snr_idx)/10))/2;
        len = length(rx_signal);
        noise = sqrt(noise_var) * (randn(len,1) + 1j*randn(len,1));
        % 添加噪声
        rx_signal = rx_signal + noise;

        long_tr_syms = rx_signal(33:end); % 去掉长训练前的保护间隔
        long_tr_syms = reshape(long_tr_syms, 64, 2); % 将两个长训练序列分开在两列

        % 转换到频域
        freq_long_tr = fft(long_tr_syms)/(64/sqrt(52));
        freq_long_tr(reorder,:) = freq_long_tr;
        freq_tr_syms = freq_long_tr(UsedSubcIdx,:);

        channel_est = mean(freq_tr_syms,2).*conj(LongTrain); 
        err_est(pkt_idx) = mean(abs(channel_est-channel).^2)/mean(abs(channel).^2);
    end
    mse(snr_idx) = mean(err_est);
end
semilogy(snr,mse,'-b.');hold on;
xlabel('SNR/dB');
ylabel('MSE');

函数tx_freqd_to_timed.m

function time_syms = tx_freqd_to_timed(mod_ofdm_syms) 
num_symbols = size(mod_ofdm_syms,2);

UsedSubcIdx = [7:32 34:59];
resample_patt=[33:64 1:32];

syms_into_ifft = zeros(64, num_symbols);
syms_into_ifft(UsedSubcIdx,:) = mod_ofdm_syms;

syms_into_ifft(resample_patt,:) = syms_into_ifft;

% Convert to time domain
time_syms = sqrt(64)*ifft(sqrt(64/52)*syms_into_ifft);

仿真结果

信道估计_SNR/MSE曲线

随着信噪比SNR的增大，均方误差MSE不断减小，符合预期。

完整OFDM系统流程

完整仿真流程

仿真结果

OFDM完整系统_仿真曲线

上图中，红色曲线为误包率（PER）随着信噪比（SNR）变化的曲线，绿色曲线为误码率（BER）随着信噪比（SNR）变化的曲线。可以看出，随着SNR的增加，PER和BER都是在渐渐减小。而且PER总是大于BER，在信噪比较低时，PER能够达到百分之百。

程序实现

CompleteOFDM.m

%% ************** 准备部分 ********************
clear all; clc;
% 系统参数
fs = 8e6;
ml = 2;                      % 调制方式: 2--4QAM; 4--16QAM; 6--64QAM
NormFactor = sqrt(2/3*(ml.^2-1));
gi = 1/4;                   % 保护间隔: 最大时延为2e-6s，相当于16个点. 
fftlen = 64;
gilen = gi*fftlen;           % 保护间隔长度 (点)
blocklen = fftlen + gilen;   % CP块的总长度

% 参数定义

DataSubcPatt = [1:5 7:19 21:26 27:32 34:46 48:52]';
PilotSubcPatt = [6 20 33 47];
UsedSubcIdx = [7:32 34:59];

% 信道编码参数
%trellis = poly2trellis([4 3],[4 5 17;7 4 2]); % 2/3
trellis = poly2trellis(7,[133 171]); % 卷积码1/2编码
% tb是指定反馈深度的正整数标量
% 如果码率是1/2, tb的一个典型的值是5次。
% 编码的约束长度是K (这里K = 7). 
tb = 7*5;
ConvCodeRate = 1/2;       % 卷积码编码码率：1/2, 2/3

% 训练定义
% 用于CFO估计的短训练序列 (符号数量NumSymbols = 52)
ShortTrain = sqrt(13/6)*[0 0 1+j 0 0 0 -1-j 0 0 0 1+j 0 0 0 -1-j 0 0 0 -1-j 0 ...
0 0 1+j 0 0 0 0 0 0 -1-j 0 0 0 -1-j 0 0 0 1+j 0 0 0 1+j 0 0 0 1+j 0 0 0 1+j 0 0].';
NumShortTrainBlks = 10;
NumShortComBlks = 16*NumShortTrainBlks/blocklen;

% 用于时间精同步和信道估计的长训练序列 (NumSymbols = 52)
LongTrain = [1 1 -1 -1 1 1 -1 1 -1 1 1 1 1 1 1 -1 -1 1 1 -1 1 -1 1 1 1 1 ...
      1 -1 -1 1 1 -1 1 -1 1 -1 -1 -1 -1 -1 1 1 -1 -1 1 -1 1 -1 1 1 1 1].';
NumLongTrainBlks = 2;
NumTrainBlks = NumShortComBlks + NumLongTrainBlks; 

% 前导码生成
short_train = tx_freqd_to_timed(ShortTrain);
%plot(abs(short_train));
short_train_blk = short_train(1:16);
short_train_blks = repmat(short_train_blk,NumShortTrainBlks,1);

long_train = tx_freqd_to_timed(LongTrain);
long_train_syms = [long_train(fftlen-2*gilen+1:fftlen,:); long_train; long_train];
preamble = [short_train_blks; long_train_syms];

% 包信息
NumBitsPerBlk = 48*ml*ConvCodeRate;
NumBlksPerPkt = 50;
NumBitsPerPkt = NumBitsPerBlk*NumBlksPerPkt;
NumPkts = 250;

% 信道和CFO
h = zeros(gilen,1);
h(1) = 1; h(3) = 0.5; h(5) = 0.2;
h = h/norm(h);
CFO = 0.1*fs/fftlen;

% 时域参数
ExtraNoiseSamples = 500;

%% ************** 循环启动 ***************************
snr = 1:2:20;
ber = zeros(1,length(snr));
per = zeros(1,length(snr));
for snr_index = 1:length(snr)
    num_err = 0;
    err = zeros(1,NumPkts);
    for pkt_index = 1:NumPkts
        [snr_index pkt_index]
%% *********************** 发送 ******************************
        % 生成信息比特
        inf_bits = randn(1,NumBitsPerPkt)>0;
        CodedSeq = convenc(inf_bits,trellis); % 对输入数据进行卷积编码

        % 调制
        paradata = reshape(CodedSeq,length(CodedSeq)/ml,ml);
        ModedSeq = qammod(bi2de(paradata),2^ml)/NormFactor;

        mod_ofdm_syms = zeros(52, NumBlksPerPkt);
        mod_ofdm_syms(DataSubcPatt,:) = reshape(ModedSeq, 48, NumBlksPerPkt);
        mod_ofdm_syms(PilotSubcPatt,:) = 1;

        tx_blks = tx_freqd_to_timed(mod_ofdm_syms);

        % 保护间隔插入
        tx_frames = [tx_blks(fftlen-gilen+1:fftlen,:); tx_blks];
        % 并串转换
        tx_seq = reshape(tx_frames,NumBlksPerPkt*blocklen,1);
        tx = [preamble;tx_seq]; % 前导训练序列和有用信息连在一起

%% ****************************** 信道 ****************************
        FadedSignal = filter(h,1,tx);
        len = length(FadedSignal);
        noise_var = 1/(10^(snr(snr_index)/10))/2;
        noise = sqrt(noise_var) * (randn(len,1) + j*randn(len,1));
        % 添加噪声
        rx_signal = FadedSignal + noise;

        % 提取的噪声样本在包前被插入，来检测包搜索算法
        extra_noise = sqrt(noise_var) * (randn(ExtraNoiseSamples,1) + j*randn(ExtraNoiseSamples,1));
        % 增加端部噪声，防止模拟因接收器中的错误定时而崩溃
        end_noise = sqrt(noise_var) * (randn(170,1) + j*randn(170,1));

        rx = [extra_noise; rx_signal; end_noise];

        % 引入 CFO
        total_length = length(rx);
        t = [0:total_length-1]/fs;
        phase_shift = exp(j*2*pi*CFO*t).';
        rx = rx.*phase_shift;

%% *************************  接收端  ****************************
        % 包搜索
        rx_signal = pkt_dect(rx);

        % CFO 估计和纠正
        rx_signal = frequency_sync(rx_signal,fs);

        % 时间精同步
        fine_time_est = fine_time_sync(rx_signal, long_train);

        % 时间同步后的信号
        expected_length = 64*2+80*NumBlksPerPkt;
        sync_time_signal = rx_signal(fine_time_est:fine_time_est+expected_length-1);

        [freq_tr_syms, freq_data_syms, freq_pilot_syms] = timed_to_freqd(sync_time_signal);     

        % 信道估计和均衡
        channel_est = mean(freq_tr_syms,2).*conj(LongTrain);         

        % 数据符号信道纠正
        chan_corr_mat = repmat(channel_est(DataSubcPatt), 1, size(freq_data_syms,2));
        freq_data_syms = freq_data_syms.*conj(chan_corr_mat);
        chan_corr_mat = repmat(channel_est(PilotSubcPatt), 1, size(freq_pilot_syms,2));
        freq_pilot_syms = freq_pilot_syms.*conj(chan_corr_mat);

        % 幅度归一化
        chan_sq_amplitude = sum(abs(channel_est(DataSubcPatt,:)).^2, 2);
        chan_sq_amplitude_mtx = repmat(chan_sq_amplitude,1, size(freq_data_syms,2));
        data_syms_out = freq_data_syms./chan_sq_amplitude_mtx;

        chan_sq_amplitude = sum(abs(channel_est(PilotSubcPatt,:)).^2, 2);
        chan_sq_amplitude_mtx = repmat(chan_sq_amplitude,1, size(freq_pilot_syms,2));
        pilot_syms_out = freq_pilot_syms./chan_sq_amplitude_mtx;

        phase_est = angle(sum(pilot_syms_out));
        phase_comp = exp(-j*phase_est);
        data_syms_out = data_syms_out.*repmat(phase_comp,48,1);

        Data_seq = reshape(data_syms_out,48*NumBlksPerPkt,1);

        % 解调
        DemodSeq = de2bi(qamdemod(Data_seq*NormFactor,2^ml),ml);
        deint_bits = reshape(DemodSeq,size(DemodSeq,1)*ml,1).';

        % 维特比译码
        DecodedBits = vitdec(deint_bits(1:length(CodedSeq)),trellis,tb,'trunc','hard');
        % 错误计算
        err(pkt_index) = sum(abs(DecodedBits-inf_bits));
        num_err = num_err + err(pkt_index);
    end
    ber(snr_index) = num_err/(NumPkts*NumBitsPerPkt);
    per(snr_index) = length(find(err~=0))/NumPkts;
end

%% display SNR-BER
semilogy(snr,ber,'-g.');hold on; %SNR/BER: 误码率
semilogy(snr,per,'-r'); %SNR/PER: 误包率
xlabel('SNR (dB)');
ylabel('BER and PER');

函数tx_freqd_to_timed.m

function time_syms = tx_freqd_to_timed(mod_ofdm_syms) 
num_symbols = size(mod_ofdm_syms,2);

UsedSubcIdx = [7:32 34:59];
resample_patt=[33:64 1:32];

syms_into_ifft = zeros(64, num_symbols);
syms_into_ifft(UsedSubcIdx,:) = mod_ofdm_syms;

syms_into_ifft(resample_patt,:) = syms_into_ifft;

% Convert to time domain
time_syms = sqrt(64)*ifft(sqrt(64/52)*syms_into_ifft);

函数timed_to_freqd.m

function [freq_tr_syms, freq_data_syms, freq_pilot_syms] = rx_timed_to_freqd(time_signal)

UsedSubcIdx = [7:32 34:59];
DataSubcIdx = [7:11 13:25 27:32 34:39 41:53 55:59];
PilotSubcIdx = [12 26 40 54];

% Long Training symbols
long_tr_syms = time_signal(1:2*64);
long_tr_syms = reshape(long_tr_syms, 64, 2);

% to frequency domain
freq_long_tr = fft(long_tr_syms)/(64/sqrt(52));
reorder = [33:64 1:32];
freq_long_tr(reorder,:) = freq_long_tr;

% Select training carriers
freq_tr_syms = freq_long_tr(UsedSubcIdx,:);

% Take data symbols
data_syms = time_signal(129:length(time_signal));

data_sig_len = length(data_syms);
n_data_syms = floor(data_sig_len/80);

% Cut to multiple of symbol period
data_syms = data_syms(1:n_data_syms*80);
data_syms = reshape(data_syms, 80, n_data_syms);
% remove guard intervals
data_syms(1:16,:) = [];

% perform fft
freq_data = fft(data_syms)/(64/sqrt(52));

%Reorder pattern is [33:64 1:32]
freq_data(reorder,:) = freq_data;

%Select data carriers
freq_data_syms = freq_data(DataSubcIdx,:);

%Select the pilot carriers
freq_pilot_syms = freq_data(PilotSubcIdx,:);

函数pkt_dect.m

function detected_packet = rx_find_packet_edge(rx_signal)

search_win = 800;
D = 16;

% Calculate the delayd correlation
delay_xcorr = rx_signal(1:search_win+2*D).*conj(rx_signal(1*D+1:search_win+3*D));

% Moving average of the delayed correlation
ma_delay_xcorr = abs(filter(ones(1,2*D), 1, delay_xcorr));

% Moving average of received power
ma_rx_pwr = filter(ones(1,2*D), 1, abs(rx_signal(1*D+1:search_win+3*D)).^2);

% The decision variable
delay_len = length(ma_delay_xcorr);
ma_M = ma_delay_xcorr(1:delay_len)./ma_rx_pwr(1:delay_len);

% remove delay samples
ma_M(1:2*D) = [];

threshold = 0.75;

thres_idx = find(ma_M > threshold);
if isempty(thres_idx)
  thres_idx = 1;
else
  thres_idx = thres_idx(1);
end


detected_packet = rx_signal(thres_idx:length(rx_signal));

函数frequency_sync.m

function out_signal = rx_frequency_sync(rx_signal,fs)

pkt_det_offset = 20;

% averaging length
rlen = 160-pkt_det_offset;

% short training symbol periodicity
D = 16;

phase = rx_signal(pkt_det_offset:pkt_det_offset+rlen-D).* ...
  conj(rx_signal(pkt_det_offset+D:pkt_det_offset+rlen));

% add all estimates 
phase = sum(phase);

freq_est = -angle(phase) / (2*D*pi/fs);
radians_per_sample = 2*pi*freq_est/fs;
time_base = 0:length(rx_signal)-1;
correction = exp(-1j*(radians_per_sample)*time_base);             
out_signal = rx_signal.*correction.';

函数fine_time_sync.m

function fine_time_est = rx_fine_time_sync(input_signal,long_train);

%timing search window size
start_search=150;
end_search=210;
time_corr_long = zeros(1,end_search-start_search+1);

for idx=start_search:end_search
    time_corr_long(idx-start_search+1) = sum((input_signal(idx:idx+63).*conj(long_train)));
end

[max_corr,long_search_idx] = max(abs(time_corr_long));

fine_time_est = start_search-1 + long_search_idx;
文章作者: Mat Jenin
文章链接: http://matjenin.xyz
版权声明: 本博客所有文章除特別声明外，均采用 CC BY 4.0 许可协议。转载请注明来源 Mat Jenin !
通信系统集成电路设计
赏
  目录