forked from bb16177/OTFS-Simulation
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmain.m
276 lines (225 loc) · 14.9 KB
/
main.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
clear;
close all;
%--------------------------------------------------------------------------
%
% This code forms a simulation of a wideband wireless communications system
% in multipath fading. It simulates: OFDM, OTFS, coded-OFDM and coded-OTFS
%
% The following files are required:
% dataGen.m
% multipathChannel.m
% modOFDM.m
% demodOFDM.m
% ISFFT.m
% SFFT.m
% equaliser.m
% plotGraphs.m
%
%--------------------------------------------------------------------------
%
% Author: Bradley Bates
% University of Bristol, UK
% email address: [email protected]
% May 2020
%
% Copyright (c) 2020, Bradley Bates
%
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% Define simulation parameters
%--------------------------------------------------------------------------
M = 16; % Modulation alphabet
k = log2(M); % Bits/symbol
cpSize = 0.07; % OFDM cyclic prefix size
scs = 15e3; % Subcarrier spacing Hz
Bw = 10e6; % System bandwidth Hz
ofdmSym = 14; % No. OFDM symbols / subframe
EbNo = (-3:1:30)'; % Range of energy/bit to noise power ratio
velocity = 120; % Velocity of mobile rx relative tx km/hr
codeRate = 2/4; % FEC code rate used
maxIterations = 25; % Set maximum no. of iterations for LDPC decoder
totalBits = 1e6; % The approx. total no. of bits simulated
repeats = 1; % Number of simulation repetitions
%--------------------------------------------------------------------------
% Initialise Simulation Components
%--------------------------------------------------------------------------
% Initialise OFDM Mod/Demod variables
numSC = pow2(ceil(log2(Bw/scs))); % Calc. nearest base-2 no. of OFDM subcarriers
cpLen = floor(cpSize * numSC); % Calc. cyclic prefix length
numDC = (numSC - 12); % Calc. number of data carriers
% Initialise the AWGN channel
awgnChannel = comm.AWGNChannel('NoiseMethod','Variance', 'VarianceSource','Input port');
errorRate = comm.ErrorRate('ResetInputPort',true);
errorRate1 = comm.ErrorRate('ResetInputPort',true);
% Initialise the LDPC coder/decoder
parityCheck_matrix = dvbs2ldpc(codeRate); % Generate DVB-S.2 paraity check matrix
ldpcEncoder = comm.LDPCEncoder(parityCheck_matrix); % Create encoder system object
ldpcDecoder = comm.LDPCDecoder(parityCheck_matrix); % Create decoder system object
ldpcDecoder.MaximumIterationCount = maxIterations; % Assign decoder's maximum iterations
noCodedbits = size(parityCheck_matrix,2); % Find the Codeword length
% Create Vectors for storing error data
berOFDM = zeros(length(EbNo),3); berCOFDM = zeros(length(EbNo),3); berOTFS = zeros(length(EbNo),3); berCOTFS = zeros(length(EbNo),3);
errorStats_coded = zeros(1,3); errorStats_uncoded = zeros(1,3);
for repetition=1:repeats % Repeat simulation multiple times with a unqique channel for each repeat
% Generate and Encode data
[dataIn, dataBits_in, codedData_in, packetSize, numPackets, numCB] = dataGen(k,numDC,ofdmSym,totalBits,codeRate,ldpcEncoder);
% Generate Rayleigh Fading Channel Impulse Response
txSig_size = zeros((numSC+cpLen),ofdmSym); % Assign the size of the channel
rayChan = multipathChannel(cpSize, scs, txSig_size, velocity); % Create fading channel impulse response
% QAM Modulator
qamTx = qammod(dataIn,M,'InputType','bit','UnitAveragePower',true); % Apply QAM modulation
parallelTx = reshape(qamTx,[numDC,ofdmSym*packetSize]); % Convert to parallel
% Add nulls at index 1
guardbandTx = [zeros(1,ofdmSym*packetSize); parallelTx];
% Add 11 nulls around DC
guardbandTx = [guardbandTx(1:(numDC/2),:); zeros(11,ofdmSym*packetSize); guardbandTx((numDC/2)+1:end,:)];
%--------------------------------------------------------------------------
% OFDM BER Calculation
%--------------------------------------------------------------------------
% Calculate SNR
snr = EbNo + 10*log10(codeRate*k) + 10*log10(numDC/((numSC)));
% Multicarrier Modulation
frameBuffer = guardbandTx; % Create a 'buffer' so subframes can be individually modulated
txframeBuffer = []; % Initilaise matrix
for w = 1:packetSize
ofdmTx = modOFDM(frameBuffer(:,1:ofdmSym),numSC,cpLen,ofdmSym); % Apply OFDM modulation to a subframe of data
frameBuffer(:, 1:ofdmSym) = []; % Delete modulated data from frameBuffer
txframeBuffer = [txframeBuffer;ofdmTx]; % Add modulated subframe to transmit buffer
end
% Loop through different values of EbNo
for m = 1:length(EbNo)
% Loop through the of packets to be transmitted
for j = 1:numPackets
rxframeBuffer = []; % Initialise matrix
% Transmit each subframe individually
for u = 1:packetSize
% Remove next subframe from the transmit buffer
txSig = txframeBuffer( ((u-1)*numel(ofdmTx)+1) : u*numel(ofdmTx) );
% Apply Channel to input signal
fadedSig = zeros(size(txSig)); % Pre-allocate vector size
for i = 1:size(txSig,1) % Perform elementwise...
for j = 1:size(txSig,2) % ...matrix multiplication
fadedSig(i,j) = txSig(i,j).*rayChan(i,j);
end
end
% AWGN Channel
release(awgnChannel);
powerDB = 10*log10(var(fadedSig)); % Calculate Tx signal power
noiseVar = 10.^(0.1*(powerDB-snr(m))); % Calculate the noise variance
rxSig = awgnChannel(fadedSig,noiseVar); % Pass the signal through a noisy channel
% Equalisation
eqSig = equaliser(rxSig,fadedSig,txSig,ofdmSym);
% Demodulation
rxSubframe = demodOFDM(eqSig,cpLen,ofdmSym); % Apply OFDM demodulation
rxframeBuffer = [rxframeBuffer';rxSubframe']'; % Store demodulated subframe in rx buffer
end
% Remove all null carriers
parallelRx = rxframeBuffer;
parallelRx((numDC/2)+1:(numDC/2)+11, :) = []; % Remove nulls around the DC input
parallelRx(1:1, :) = []; % Remove nulls at index 1
qamRx = reshape(parallelRx,[numel(parallelRx),1]);% Convert to serial
% Uncoded demodulation of entire packet
dataOut = qamdemod(qamRx,M,'OutputType','bit','UnitAveragePower',true);% Apply QAM demodulation
codedData_out = randdeintrlv(dataOut,4831); % De-interleave data
codedData_out(numel(codedData_in)+1:end) = []; % Remove pad bits
errorStats_uncoded = errorRate(codedData_in,codedData_out,0); % Collect error statistics
% Coded demodulation of entire packet
powerDB = 10*log10(var(qamRx)); % Calculate Rx signal power
noiseVar = 10.^(0.1*(powerDB-(EbNo(m) + 10*log10(codeRate*k) - 10*log10(sqrt(numDC))))); % Calculate the noise variance
dataOut = qamdemod(qamRx,M,'OutputType', 'approxllr','UnitAveragePower',true,'NoiseVariance',noiseVar);% Apply QAM demodulation
codedData_out1 = randdeintrlv(dataOut,4831); % De-interleave data
codedData_out1(numel(codedData_in)+1:end) = []; % Remove pad bits
% Decode individual code blocks
dataBits_out = []; % Initialise matrix
dataOut_buffer = codedData_out1;
for q = 1:numCB
dataBits_out = [dataBits_out;ldpcDecoder(dataOut_buffer(1:noCodedbits))]; % Decode data & add it to the data bits out matrix
dataOut_buffer(1:noCodedbits) = []; % Delete decoded data from buffer
end
dataBits_out = double(dataBits_out); % Convert to a double compatible w/ errorStats
errorStats_coded = errorRate1(dataBits_in,dataBits_out,0); % Collect error statistics
end
berOFDM(m,:) = errorStats_uncoded; % Save uncoded BER data
berCOFDM(m,:) = errorStats_coded; % Save coded BER data
errorStats_uncoded = errorRate(codedData_in,codedData_out,1); % Reset the error rate calculator
errorStats_coded = errorRate1(dataBits_in,dataBits_out,1); % Reset the error rate calculator
end
%--------------------------------------------------------------------------
% OTFS BER Calculation
%--------------------------------------------------------------------------
% Calculate SNR
snr = EbNo + 10*log10(codeRate*k) + 10*log10(numDC/((numSC))) + 10*log10(sqrt(ofdmSym));
% Multicarrier Modulation
frameBuffer = guardbandTx; % Create a 'buffer' so subframes can be individually modulated
txframeBuffer = []; % Initilaise matrix
for w = 1:packetSize
otfsTx = ISFFT(frameBuffer(:,1:ofdmSym)); % Apply OTFS modulation to a subframe of data
ofdmTx = modOFDM(otfsTx,numSC,cpLen,ofdmSym); % Apply OFDM modulation
frameBuffer(:, 1:ofdmSym) = []; % Delete modulated data from frameBuffer
txframeBuffer = [txframeBuffer;ofdmTx]; % Add modulated subframe to transmit buffer
end
% Loop through different values of EbNo
for m = 1:length(EbNo)
% Loop through the of packets to be transmitted
for j = 1:numPackets
rxframeBuffer = []; % Initialise matrix
% Transmit each subframe individually
for u = 1:packetSize
% Remove next subframe from the transmit buffer
txSig = txframeBuffer( ((u-1)*numel(ofdmTx)+1) : u*numel(ofdmTx) );
% Apply Channel to input signal
fadedSig = zeros(size(txSig)); % Pre-allocate vector size
for i = 1:size(txSig,1) % Perform elementwise...
for j = 1:size(txSig,2) % ...matrix multiplication
fadedSig(i,j) = txSig(i,j).*rayChan(i,j);
end
end
% AWGN Channel
release(awgnChannel);
powerDB = 10*log10(var(fadedSig)); % Calculate Tx signal power
noiseVar = 10.^(0.1*(powerDB-snr(m))); % Calculate the noise variance
rxSig = awgnChannel(fadedSig,noiseVar); % Pass the signal through a noisy channel
% Equalisation
eqSig = equaliser(rxSig,fadedSig,txSig,ofdmSym);
% Demodulation
otfsRx = demodOFDM(eqSig,cpLen,ofdmSym); % Apply OFDM demodulation
rxSubframe = SFFT(otfsRx); % Apply OTFS demodulation
rxframeBuffer = [rxframeBuffer';rxSubframe']'; % Store demodulated subframe in rx buffer
end
% Remove all null carriers
parallelRx = rxframeBuffer;
parallelRx((numDC/2)+1:(numDC/2)+11, :) = []; % Remove nulls around the DC input
parallelRx(1:1, :) = []; % Remove nulls at index 1
qamRx = reshape(parallelRx,[numel(parallelRx),1]); % Convert to serial
% Uncoded demodulation of entire packet
dataOut = qamdemod(qamRx,M,'OutputType','bit','UnitAveragePower',true);% Apply QAM demodulation
codedData_out = randdeintrlv(dataOut,4831); % De-interleave data
codedData_out(numel(codedData_in)+1:end) = []; % Remove pad bits
errorStats_uncoded = errorRate(codedData_in,codedData_out,0); % Collect error statistics
% Coded demodulation of entire packet
powerDB = 10*log10(var(qamRx)); % Calculate Rx signal power
noiseVar = 10.^(0.1*(powerDB-(EbNo(m) + 10*log10(codeRate*k) - 10*log10(sqrt(numDC))))); % Calculate the noise variance
dataOut = qamdemod(qamRx,M,'OutputType', 'approxllr','UnitAveragePower',true,'NoiseVariance',noiseVar);% Apply QAM demodulation
codedData_out1 = randdeintrlv(dataOut,4831); % De-interleave data
codedData_out1(numel(codedData_in)+1:end) = []; % Remove pad bits
% Decode individual code blocks
dataBits_out = []; % Initialise matrix
dataOut_buffer = codedData_out1;
for q = 1:numCB
dataBits_out = [dataBits_out;ldpcDecoder(dataOut_buffer(1:noCodedbits))]; % Decode data & add it to the data bits out matrix
dataOut_buffer(1:noCodedbits) = []; % Delete decoded data from buffer
end
dataBits_out = double(dataBits_out); % Convert to a double compatible w/ errorStats
errorStats_coded = errorRate1(dataBits_in,dataBits_out,0); % Collect error statistics
end
berOTFS(m,:) = errorStats_uncoded; % Save uncoded BER data
berCOTFS(m,:) = errorStats_coded; % Save coded BER data
errorStats_uncoded = errorRate(codedData_in,codedData_out,1); % Reset the error rate calculator
errorStats_coded = errorRate1(dataBits_in,dataBits_out,1); % Reset the error rate calculator
end
end
%--------------------------------------------------------------------------
% Figures
%--------------------------------------------------------------------------
% Plot BER / EbNo curves
plotGraphs(berOFDM, berCOFDM, berOTFS, berCOTFS, M, numSC, EbNo);