文献来源:微信公众号:EW Frontier
OTFS简介
OTFS信道估计
% Clear command window, workspace variables, and close all figures
clc;
clear all;
close all;
% Define Eb values in dB
EbdB = -10:2:10;
% Convert Eb values from dB to linear scale
Eb = 10.^(EbdB/10);
% Define Noise Power
No = 1;
% Calculate Signal-to-Noise Ratio (SNR) in linear scale
SNR = 2*Eb/No;
% Convert SNR to dB scale
SNRdB = 10*log10(SNR);
% Define matrix dimensions and parameters
M = 32;
N = 16;
Ptx = eye(M);
Prx = eye(M);
nTaps = 5;
DelayTaps = [0 1 2 3 4];
DopplerTaps = [0 1 2 3 4];
Ncp = max(DelayTaps);
% Initialize arrays to store Bit Error Rate (BER) for different methods
BER_MMSE = zeros(length(Eb),1);
BER_ZF = zeros(length(Eb),1);
% Number of iterations for Monte Carlo simulation
ITER = 10;
% Precompute matrices for transformation
F_M = 1/sqrt(M)*dftmtx(M);
F_N = 1/sqrt(N)*dftmtx(N);
% Main loop for Monte Carlo simulation
for ite = 1:ITER
ite
% Generate random bits for transmission
XddBits = randi([0,1],M,N);
% Generate random channel taps
h = sqrt(1/2)*(randn(1,nTaps)+ 1j*randn(1,nTaps));
% Construct effective channel matrix
Hmat = zeros(M*N,M*N);
omega = exp(1j*2*pi/(M*N));
for tx = 1:nTaps
Hmat = Hmat + h(tx)*circshift(eye(M*N),DelayTaps(tx))*...
(diag(omega.^((0:M*N-1)*DopplerTaps(tx))));
end
Heff = kron(F_N,Prx)*Hmat*kron(F_N',Ptx);
% Generate Complex Noise
ChNoise = sqrt(No/2)*(randn(1,M*N) + 1j*randn(1,M*N));
% Loop over different Eb/N0 values
for ix = 1:length(Eb)
% Generate modulated symbols
X_DD = sqrt(Eb(ix))*(2*XddBits-1);
X_TF = F_M*X_DD*F_N';
S_mat = Ptx*F_M'*X_TF;
TxSamples = reshape(S_mat,M*N,1).';
TxSamplesCP = [TxSamples(M*N-Ncp+1:M*N) TxSamples];
% Channel filtering
RxsamplesCP = 0;
for tx = 1:nTaps
Doppler = exp(1j*2*pi/M*(-Ncp:M*N-1)*DopplerTaps(tx)/N);
RxsamplesCP = RxsamplesCP + h(tx)*circshift(TxSamplesCP.*Doppler,[1, DelayTaps(tx)]);
end
% Remove cyclic prefix
Rxsamples = RxsamplesCP(Ncp+1:M*N+Ncp) + ChNoise;
R_mat = reshape(Rxsamples.',M, N) ;
Y_TF = F_M*Prx*R_mat;
Y_DD = F_M'*Y_TF*F_N;
y_DD = reshape(Y_DD, M*N, 1) ;
% MMSE Equalization
xhatMMSE = inv(Heff'*Heff + eye(M*N)/Eb(ix))*Heff'*y_DD;
DecodedBits = (real(xhatMMSE) >= 0);
BER_MMSE(ix) = BER_MMSE(ix) + sum(DecodedBits ~= reshape(XddBits,M*N,1));
% Zero Forcing (ZF) Equalization
xhatZF = pinv(Heff)*y_DD;
DecodedBits = (real(xhatZF) >= 0);
BER_ZF(ix) = BER_ZF(ix) + sum(DecodedBits ~= reshape(XddBits,M*N,1));
end
end
% Average BER over iterations and symbols
BER_MMSE = BER_MMSE/M/N/ITER;
BER_ZF = BER_ZF/M/N/ITER;
% Plot BER versus SNR
semilogy(EbdB,BER_MMSE,'b-s','linewidth',3.0,'MarkerFaceColor','b','MarkerSize',9.0);
hold on;
grid on;
semilogy(EbdB,BER_ZF,'r-o','linewidth',3.0,'MarkerFaceColor','r','MarkerSize',9.0);
axis tight;
legend('MMSE','ZF');
title('OTFS BER v/s SNR');
xlabel('SNR(dB)');
ylabel('BER');
OTFS系统建模
% Clear command window, workspace variables, and close all figures
clc;
clear all;
close all;
% Define Eb values in dB
EbdB = -10:2:10;
% Convert Eb values from dB to linear scale
Eb = 10.^(EbdB/10);
% Define Noise Power
No = 1;
% Calculate Signal-to-Noise Ratio (SNR) in linear scale
SNR = 2*Eb/No;
% Convert SNR to dB scale
SNRdB = 10*log10(SNR);
% Define matrix dimensions and parameters
M = 32;
N = 16;
Ptx = eye(M);
Prx = eye(M);
nTaps = 5;
DelayTaps = [0 1 2 3 4];
DopplerTaps = [0 1 2 3 4];
Ncp = max(DelayTaps);
% Initialize arrays to store Bit Error Rate (BER) for different methods
BER_MMSE = zeros(length(Eb),1);
BER_ZF = zeros(length(Eb),1);
% Number of iterations for Monte Carlo simulation
ITER = 10;
% Precompute matrices for transformation
F_M = 1/sqrt(M)*dftmtx(M);
F_N = 1/sqrt(N)*dftmtx(N);
% Main loop for Monte Carlo simulation
for ite = 1:ITER
ite
% Generate random bits for transmission
XddBits = randi([0,1],M,N);
% Generate random channel taps
h = sqrt(1/2)*(randn(1,nTaps)+ 1j*randn(1,nTaps));
% Construct effective channel matrix
Hmat = zeros(M*N,M*N);
omega = exp(1j*2*pi/(M*N));
for tx = 1:nTaps
Hmat = Hmat + h(tx)*circshift(eye(M*N),DelayTaps(tx))*...
(diag(omega.^((0:M*N-1)*DopplerTaps(tx))));
end
Heff = kron(F_N,Prx)*Hmat*kron(F_N',Ptx);
% Generate Complex Noise
ChNoise = sqrt(No/2)*(randn(1,M*N) + 1j*randn(1,M*N));
% Loop over different Eb/N0 values
for ix = 1:length(Eb)
% Generate modulated symbols
X_DD = sqrt(Eb(ix))*(2*XddBits-1);
X_TF = F_M*X_DD*F_N';
S_mat = Ptx*F_M'*X_TF;
TxSamples = reshape(S_mat,M*N,1).';
TxSamplesCP = [TxSamples(M*N-Ncp+1:M*N) TxSamples];
% Channel filtering
RxsamplesCP = 0;
for tx = 1:nTaps
Doppler = exp(1j*2*pi/M*(-Ncp:M*N-1)*DopplerTaps(tx)/N);
RxsamplesCP = RxsamplesCP + h(tx)*circshift(TxSamplesCP.*Doppler,[1, DelayTaps(tx)]);
end
% Remove cyclic prefix
Rxsamples = RxsamplesCP(Ncp+1:M*N+Ncp) + ChNoise;
R_mat = reshape(Rxsamples.',M, N) ;
Y_TF = F_M*Prx*R_mat;
Y_DD = F_M'*Y_TF*F_N;
y_DD = reshape(Y_DD, M*N, 1) ;
% MMSE Equalization
xhatMMSE = inv(Heff'*Heff + eye(M*N)/Eb(ix))*Heff'*y_DD;
DecodedBits = (real(xhatMMSE) >= 0);
BER_MMSE(ix) = BER_MMSE(ix) + sum(DecodedBits ~= reshape(XddBits,M*N,1));
% Zero Forcing (ZF) Equalization
xhatZF = pinv(Heff)*y_DD;
DecodedBits = (real(xhatZF) >= 0);
BER_ZF(ix) = BER_ZF(ix) + sum(DecodedBits ~= reshape(XddBits,M*N,1));
end
end
% Average BER over iterations and symbols
BER_MMSE = BER_MMSE/M/N/ITER;
BER_ZF = BER_ZF/M/N/ITER;
% Plot BER versus SNR
semilogy(EbdB,BER_MMSE,'b-s','linewidth',3.0,'MarkerFaceColor','b','MarkerSize',9.0);
hold on;
grid on;
semilogy(EbdB,BER_ZF,'r-o','linewidth',3.0,'MarkerFaceColor','r','MarkerSize',9.0);
axis tight;
legend('MMSE','ZF');
title('OTFS BER v/s SNR');
xlabel('SNR(dB)');
ylabel('BER');
OTFS系统建模
clc; clear all; close all;
% SIGNAL MODEL OF OTFS
M = 32; N = 32; %M-> no. of subcarriers/delay bins ; N-> no. of symbols/doppler bins
F_M = 1/sqrt(M)*dftmtx(M); %disrete fourier transform matrix (normalised)
F_N = 1/sqrt(N)*dftmtx(N); %forms the grid
Ptx = eye(M);
delta_f = 15e3; %subcarrier BW
T = 1/delta_f; %OFDM symbol Duration
X_DD = zeros(M,N); %Delay Doppler Grid
X_DD(3, 3) = 1; %impulse in DD Domain
X_TF = F_M*X_DD*F_N';
S = Ptx*F_M'
*X_TF;
s = reshape(S,M*N,1);
% figure()
subplot(4,2,[1,2,3,4])
bar3(X_DD);
axis tight;
xlabel('Doppler');
ylabel('Delay');
title('Basis function in DD-domain');
% figure()
subplot(4,2,5)
surf(real(X_TF));
axis tight;
xlabel('Time');
ylabel('Subcarrier');
title('Basis function in TF-domain (Real)');
% figure()
subplot(4,2,6)
surf(imag(X_TF));
axis tight;
xlabel('Time');
ylabel('Subcarrier');
title('Basis function in TF-domain (Imag)');
% figure()
subplot(4,2,7)
plot((0:length(s)-1)*T/M,real(s));
axis tight;
ylim([-0.5 0.5]);
xlabel('Time');
title('Basis function in time-domain (Real)');
% figure()
subplot(4,2,8)
plot((0:length(s)-1)*T/M,imag(s));
axis tight;
ylim([-0.5 0.5]);
xlabel('Time');
title('Basis function in time-domain (Imag)');
OTFS系统性能部分代码
% Clear command window, workspace variables, and close all figures
clc;
clear all;
close all;
% Define Eb values in dB
EbdB = 2;
% Convert Eb values from dB to linear scale
Eb = 10.^(EbdB/10);
% Define Noise Power
No = 1;
% Calculate Signal-to-Noise Ratio (SNR) in linear scale
SNR = 2*Eb/No;
% Convert SNR to dB scale
SNRdB = 10*log10(SNR);
% Define matrix dimensions and parameters
M = 32;
N = 16;
Ptx = eye(M);
Prx = eye(M);
nTaps = 5;
DelayTaps = [5 1 0 3 4];
DopplerTaps = [0 3 2 3 4];
Ncp = max(DelayTaps);
% Precompute matrices for transformation
F_M = 1/sqrt(M)*dftmtx(M);
F_N = 1/sqrt(N)*dftmtx(N);
% Generate random bits for transmission
XddBits = randi([0,1],M,N);
% Generate random channel taps
h = sqrt(1/2)*(randn(1,nTaps)+ 1j*randn(1,nTaps));
% Construct effective channel matrix
Hmat = zeros(M*N,M*N);
omega = exp(1j*2*pi/(M*N));
for tx = 1:nTaps
Hmat = Hmat + h(tx)*circshift(eye(M*N),DelayTaps(tx))*...
(diag(omega.^((0:M*N-1)*DopplerTaps(tx))));
end
Heff = kron(F_N,Prx)*Hmat*kron(F_N',Ptx);
% Generate Complex Noise
ChNoise = sqrt(No/2)*(randn(1,M*N) + 1j*randn(1,M*N));
% Generate modulated symbols
X_DD = sqrt(Eb)*(2*XddBits-1);
X_TF = F_M*X_DD*F_N';
S_mat = Ptx*F_M'*X_TF;
TxSamples = reshape(S_mat,M*N,1).';
TxSamplesCP = [TxSamples(M*N-Ncp+1:M*N) TxSamples];
% Channel filtering
RxsamplesCP = 0;
for tx = 1:nTaps
Doppler = exp(1j*2*pi/M*(-Ncp:M*N-1)*DopplerTaps(tx)/N);
RxsamplesCP = RxsamplesCP + h(tx)*circshift(TxSamplesCP.*Doppler,[1, DelayTaps(tx)]);
end
% Remove cyclic prefix
Rxsamples = RxsamplesCP(Ncp+1:M*N+Ncp) + ChNoise;
R_mat = reshape(Rxsamples.',M, N) ;
Y_TF = F_M*Prx*R_mat;
Y_DD = F_M'*Y_TF*F_N;
y_DD = reshape(Y_DD, M*N, 1) ;
% MMSE Equalization
xhatMMSE = inv(Heff'*Heff + eye(M*N)/Eb)*Heff'*y_DD;
DecodedBits_MMSE = (real(xhatMMSE) >= 0);
DecodedBits_MMSE_reshaped = reshape(DecodedBits_MMSE,M,N);
BER_MMSE_Map = (DecodedBits_MMSE_reshaped ~= XddBits);
BER_MMSE = sum(DecodedBits_MMSE ~= reshape(XddBits,M*N,1));
OTFS模糊函数部分代码
% Clear command window, workspace variables, and close all figures
clc;
clear all;
close all;
% Define Eb values in dB
EbdB = 3;
% Convert Eb values from dB to linear scale
Eb = 10.^(EbdB/10);
% Define Noise Power
No = 1;