OTFS基带通信系统（脉冲导频，信道估计，MP解调算法）
Embedded Pilot-Aided Channel Estimation for OTFS in Delay–Doppler Channels | IEEE Journals & Magazine | IEEE Xplore
一、OTFS通信系统
如下图简要概括了OTFS基带通信系统过程，废话不多说给出完整系统详细代码。
以下仿真结果基于四抽头信道
估计信道 BER曲线
完美信道BER曲线
主函数部分代码
% Author: [YZ ]
% Date: [2024/9/8 ]
% Description: [Brief description of the code’s functionality]
% Version: [Version number, if applicable]clc
clear all
close all
tic
%% OTFS parameters%%%%%%%%%%
% number of symbol
N = 16; %码元数为8
% number of subcarriers
M = 32; %载波数为8
Nifft=N*M;% size of constellation
M_mod = 4; %进制数为4，那么一码元就包含2比特
M_bits = log2(M_mod);
% average energy per data symbol
eng_sqrt = (M_mod==2)+(M_mod~=2)*sqrt((M_mod-1)/6*(2^2)); %符号的能量
% number of symbols per frame
Protect_SubBlock_L=8;
Protect_Front=8;
data_row_range=4:12;
data_col_range=9:24;
Pilot_row_index=8;
Pilot_col_index=1;DATA_BLOCK_N=max(data_row_range)-min(data_row_range)+1;
DATA_BLOCK_M=max(data_col_range)-min(data_col_range)+1;
N_syms_perfram = DATA_BLOCK_N*DATA_BLOCK_M; %一帧中包含的符号数=单个载波包含的码元数*载波个数\
% number of bits per frame
N_bits_perfram = N_syms_perfram*2; %%信噪比范围
SNR_dB =0:5:20; %信噪比dB从0~15，间隔5dB取样，length为4
SNR = 10.^(SNR_dB/10); %dB和10进制换算
noise_var_sqrt = sqrt(1./SNR); %噪声
sigma_2 = abs(eng_sqrt*noise_var_sqrt).^2;baud_rate=6400;
tau = [0 2e-3 3e-3 4e-3 5e-3]; % 五条多径的时延
pdb = [0 -4 -8 -12 -16]; % 五条多径的平均路径增益
fdmax=6.5;
channel = comm.RayleighChannel(...
'SampleRate',baud_rate, ...
'PathDelays',tau, ... % 五条多径的时延
'AveragePathGains',pdb, ... % 五条多径的平均路径增益
'MaximumDopplerShift',fdmax, ...
'NormalizePathGains',true,...
'PathGainsOutputPort',true,...
'RandomStream','mt19937ar with seed',...
'Seed',2);%%N_fram = 10;
PAPR_dB_SUM=0;
err_ber = zeros(length(SNR_dB),1); %生成0矩阵，length返回 X 中最大数组维度的长度，此时为4,即此命令生成4*1维的0矩阵fprintf('Start to do OTFS simulation. Please wait...\n');
for iesn0 = 1:length(SNR_dB) %1~4SNR_temp = SNR_dB(iesn0); %SNR_temp取0,5,10,15SNR_temp for ifram = 1:N_fram rng(3);%% random input bits generation and 4QAM modulation%%%%%data_info_bit = randi([0,1],N_bits_perfram,1); %生成0或1的随机数，且维数是128*1，也就是生成了一个帧结构中包含的bit数data_bit = reshape(data_info_bit,N_syms_perfram,M_bits); data_temp = bi2de(data_bit);%二进制数转换成10进制数，reshape为重组数组，语法为reshape(A,a,b)将矩阵A重构为a*b的矩阵x_temp = qammod(data_temp,M_mod, 'gray'); %输出使用正交幅度调制4-QAM消息信号X的复包络,gray格雷码编码x_temp = reshape(x_temp,DATA_BLOCK_N,DATA_BLOCK_M); %将64*1的x重组成8*8数组x=zeros(N,M);pilot=zeros(N,Protect_Front);pilot(Pilot_row_index,Pilot_col_index)=10;x(:,1:Protect_Front)=pilot;x(data_row_range,data_col_range)=x_temp;%% OTFS modulation%%%%s = OTFS_modulation(N,M,x);% 使用 imagesc 函数绘制二维图% figure(3)% magnitude=abs(x);% imagesc(magnitude);% % 添加颜色条，显示颜色映射% colorbar;% % 添加标题和轴标签% title('2D Plot of Array y');% xlabel('X Axis');% ylabel('Y Axis');[taps,delay_taps,Doppler_taps,chan_coef] = OTFS_channel_gen(N,M);%% OTFS channel output%%%%%r = OTFS_channel_output(N,M,taps,delay_taps,Doppler_taps,chan_coef,sigma_2(iesn0),s);%% add noise 根据过信道信道能量 以及信噪比 计算所加噪声功率 %%%% 信道估计信道均衡%%%% y0_down_duopuley = OTFS_demodulation(N,M,r);% figure(2);% % 使用 imagesc 函数绘制二维图% magnitude=abs(y);% imagesc(magnitude);% % 添加颜色条，显示颜色映射% colorbar;% % 添加标题和轴标签% title('2D Plot of Array y');% xlabel('X Axis');% ylabel('Y Axis');threshold=0.01;k=1;for i=1:Nfor j=1:Protect_Frontif(abs(y(i,j))>max(abs(y(:)))*0.1)est_delay_tap(k)=j-1;est_doppler_tap(k)=i-Pilot_row_index;est_chan_coef(k) =y(i,j)*exp(-1i*est_doppler_tap(k));est_chan_coef(k) =y(i,j);k=k+1;endendendest_chan_coef=est_chan_coef./10;tic%% 估计信道mp检测x_est = OTFS_mp_detector(N,M,M_mod,length(est_delay_tap),est_delay_tap,est_doppler_tap,est_chan_coef,sigma_2(iesn0),y);%% 完美信道mp检测% x_est = OTFS_mp_detector(N,M,M_mod,taps,delay_taps,Doppler_taps,chan_coef,sigma_2(iesn0),y);toc;% figure(3);% % 使用 imagesc 函数绘制二维图% magnitude=abs(x_est);% imagesc(magnitude);% % 添加颜色条，显示颜色映射% colorbar;% % 添加标题和轴标签% title('2D Plot of Array y');% xlabel('X Axis');% ylabel('Y Axis');%% message passing detector%%%% %% output bits and errors count%%%%%data_demapping_temp = qamdemod(x_est,M_mod,'gray');data_demapping=zeros(DATA_BLOCK_N,DATA_BLOCK_M);data_demapping=data_demapping_temp(data_row_range,data_col_range);data_info_est = reshape(de2bi(data_demapping,M_bits),N_bits_perfram,1);errors = sum(xor(data_info_est,data_info_bit));err_rate=errors/N_bits_perfram;err_ber(iesn0) = errors + err_ber(iesn0); endPAPR=PAPR_dB_SUM/N_framPAPR_dB_SUM=0;err_ber_fram_temp = err_ber(iesn0) / N_bits_perfram / N_fram;err_ber_fram_temp %表达式加;可以隐藏输出err_ber_fram(iesn0) = err_ber_fram_temp;
endfigure(1)err_ber_fram;
semilogy(SNR_dB, err_ber_fram,'-*');
title(sprintf('OTFS vs COTFS in doppler ch'))
ylabel('BER'); xlabel('SNR in dB');
grid on
toc其它函数
function s = OTFS_modulation(N,M,x)
%% OTFS Modulation: 1. ISFFT, 2. Heisenberg transform
X = fft(ifft(x).').'/sqrt(M/N); %%%ISFFT,其中.'表示转置，
s_mat = ifft(X.')*sqrt(M); % Heisenberg transform
s = s_mat(:);%(:)表示将矩阵重构为列向量：
endfunction y = OTFS_demodulation(N,M,r)
%% OTFS demodulation: 1. Wiegner transform, 2. SFFT
r_mat = reshape(r,M,N); %重新生成8*8的接收符号
Y = fft(r_mat)/sqrt(M); % Wigner transform
Y = Y.';
y = ifft(fft(Y).').'/sqrt(N/M); % SFFT
endfunction [taps,delay_taps,Doppler_taps,chan_coef] = OTFS_channel_gen(N,M)
%% Channel for testing%%%%%
%channel wih 4 taps of uniform power%%%
taps = 4; %抽头数设置为4
delay_taps = [0 1 2 3]; %时延抽头设置为4
Doppler_taps = [0 1 2 3 ]; %多普勒域抽头设置为4
pow_prof = (1/taps) * (ones(1,taps)); %每一个抽头都分到1/4的功率
chan_coef = sqrt(pow_prof).*(sqrt(1/2) * (randn(1,taps)+1i*randn(1,taps))); %每一个抽头生成一个信道系数，一共四个信道系数
function r = OTFS_channel_output(N,M,taps,delay_taps,Doppler_taps,chan_coef,sigma_2,s)
%% wireless channel and noise
L = max(delay_taps);%最大的时延
s = [s(N*M-L+1:N*M);s];%add one cp %加入循环前缀编码（将OTFS符号数组的后四位信号复制到头部构成的一组（64+4）*1的数组）
s_chan = 0; %信道输入信号初始化
for itao = 1:tapss_chan = s_chan+chan_coef(itao)*circshift([s.*exp(1j*2*pi/M *(-L:-L+length(s)-1)*Doppler_taps(itao)/N).';zeros(delay_taps(end),1)],delay_taps(itao));%%s_chan是71*1维数组
end
noise = sqrt(sigma_2/2)*(randn(size(s_chan)) + 1i*randn(size(s_chan)));%信道噪声
r = s_chan+noise;
r = r(L+1:L+(N*M));%discard cp(去掉循环前缀，也就是输出68*1维数组的后64位）
end