-
Notifications
You must be signed in to change notification settings - Fork 19
/
simRate.m
356 lines (280 loc) · 10.8 KB
/
simRate.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
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
function [RSOFT, RHARD] = simRate(snrdb, N, K, T, Pk, b, M, sig, rec, plotout)
% Title : Acheivable rate computation: one-bit massive MIMO uplink
% File : simRate.m
% -------------------------------------------------------------------------
% Description :
% Computes the acheivable rate for PSK/QAM signalling in a one-bit
% massive MIMO uplink system. The receiver employs MRC or ZF to map the
% received signal over all antennas to a scalar output. When the channel
% fading coefficients that are required for MRC/ZF are unknown to the
% receiver, LS channel estimation is used to estimate the channel based
% on pilot symbols.
% The probabilities that are required to compute the acheivable rates are
% approximated by Monte Carlo analysis. The MRC/ZF output domain is
% discretized to a rectangular grid. For details see:
% [S. Jacobsson, "Throughput analysis of massive MIMO uplink with one-bit
% ADCs", Master's thesis, Chalmers University of Technology, Gothenburg,
% Sweden, Jan. 2015].
% If the function is called without any inputs it will run the simulation
% for some default parameters and display the receiver output.
%
% Output parameters:
% RSOFT = Rate per user [b/s/Hz] with soft decoding
% RHARD = Rate per user [b/s/Hz] with hard decoding
%
% Input parameters:
% snrdb = SNR in dB
% N = number of antennas
% K = number of users
% T = cohrence interval measured in channel uses
% Pk = number of pilots per user (0 = perfect CSIR)
% b = number of quantization bits (0 = inf.prec., 1 = one-bit)
% M = constellation size
% sig = signalling scheme (PSK or QAM)
% rec = receiver architechture (MRC or ZF)
% plotout = plot receiver output
%
% Simulation parameters
% numC = number of channel realizations (default = 100).
% numU = number of noise/interference realizations per channel realization(default = 100).
% numbin = number of output bins in discrete grid to compute soft probabilitie
%
% Display parameters:
% showAxes : display x/y-axis (defauilt = 1)
% showGrid : display MRC/ZF output discretization grid (default = 0)
% showInpt : display input alphabet (default = 0)
%
% -------------------------------------------------------------------------
% Revisions :
% Date Version Author Description
% 26-jan-15 1.0 svenja created the file
% 29-jan-15 1.1 svenja added description and comments
% -------------------------------------------------------------------------
% Copyright (c) 2015, Sven Jacobsson
% All rights reserved.
% =========================================================================
% INITIALIZATION PHASE
%-------------------------------------------------------------------------%
% include mtimesx directory for fast computation of 3-dim matrices.
addpath mtimesx;
% number of channel realisations and noise/interference realisations
numC = 200;
numU = 200;
% number of output bins per dimension for the discrete grid on which the
% soft probabilities are approximated
numbin = 100;
% display axes/grid/inputs
showAxes = 1;
showGrid = 0;
showInpt = 0;
% total number of inputs/outputs for MC simulation.
MAX = numC*numU;
tic;
fprintf('Running simRate.m for %i channel realisations and %i noise/interference realisations...\n ',numC,numU);
% set default parameters
if nargin == 0
snrdb = 0;
N = 400;
K = 1;
T = 1:1000;
Pk = 20;
b = 1;
M = 16;
sig = 'QAM';
rec = 'MRC';
% display receiver output signal
plotout = 1;
% figure handle
fig_out = figure(99); clf; set(fig_out,'name','Receiver output');
elseif nargin == 9
plotout = 0; % do not display receiver output signal
end
% SNR in dB to linear units
snr = db2pow(snrdb);
% determine input modulation and cardinality (C = input alphabet)
if strcmp(sig,'PSK')
C = pskmod(0:M-1,M,pi/M)';
elseif strcmp(sig,'QAM')
C = qammod(0:M-1,M)';
QAMnorm = sqrt(mean(abs(C).^2));
C =C/QAMnorm;
else
error('Only PSK or QAM is supported.');
end
% Rayleigh fading channel
H = (randn(N,K,numC) + 1i*randn(N,K,numC))/sqrt(2);
% TRAINING & CHANNEL ESTIMATION PHASE
%-------------------------------------------------------------------------%
% total number of pilots
Pt = K*Pk;
% channel estimatiion: if P = 0 then assume perfect CSIR, otherwise perform
% LS channel estimation based on transmitted pilots.
if Pt == 0
He = H;
else
% pilots: each user sends at different time to ensure orthogonality. The
% pilots are choosen at random from the input alphabet.
Xp = zeros(K,Pt,numC);
for k = 1:K
Xp(k,1+(k-1)*Pk:k*Pk,:) = ones(1,Pk,numC);
end
Xp = sqrt(snr) * sqrt(K) * C(randi([1,M],K,Pt,numC)) .* Xp;
% additive noise acting on pilots
Wp = (randn(N,Pt,numC) + 1i*randn(N,Pt,numC))/sqrt(2);
% received pilot sequence (unquantised)
Yp = mtimesx(H,Xp) + Wp;
% quantised pilot sequence
Rp = quantise(Yp,b);
% LS estimation
RXp = mtimesx(Rp,Xp,'c');
XXp = mtimesx(Xp,Xp,'c');
He = mtimesx(RXp,multinv(XXp));
end
% DATA TRANSMISSION & DETECTION PHASE
%-------------------------------------------------------------------------%
% soft estimate for the first user after MRC/ZF, used to approximate the
% necessary probabilities required to compute acheivable rate.
xsoft = nan(M,MAX);
for m = 1:1:M
% data symbols: all interfering users have their inputs choosen
% randomly from the input alphabet.
Xd = C(randi([1,M],K,numU,numC));
Xd(1,:,:) = C(m);
Xd = sqrt(snr)*Xd;
% additive noise acting on data symbols
Wd = (randn(N,numU,numC) + 1i*randn(N,numU,numC))/sqrt(2);
% received data sequence (unquantised)
Yd = mtimesx(H,Xd) + Wd;
% quantised data sequence
Rd = quantise(Yd,b);
% soft decoding with MRC/ZF filter. The acquired channel estimate is
% used to map the received signal vector to a scalar channel.
if strcmp(rec,'MRC')
he = He(:,1,:);
xsoft(m,:) = reshape(mtimesx(1./sum(conj(he).*he),mtimesx(he,'c',Rd)),1,[]);
elseif strcmp(rec,'ZF')
HeHe = mtimesx(He,'c',He);
Xsoft = mtimesx(mtimesx(multinv(HeHe),He,'c'),Rd);
xsoft(m,:) = reshape(Xsoft(1,:,:),1,[]);
end
end
% hard decoding
if strcmp(sig,'PSK')
xhard = pskdemod(xsoft+eps,M,pi/M);
elseif strcmp(sig,'QAM')
xhard = qamdemod(QAMnorm/sqrt(snr)*xsoft+eps,M);
end
% PROBABILITIES & RATE COMPUTATION
%-------------------------------------------------------------------------%
% soft channel law and soft output probability. "edges" is used when
% plotting the discretization grid.
[pSOFT, poutSOFT, edges] = softProb(xsoft,numbin);
% hard channel law and hard output probability.
[pHARD, poutHARD] = hardProb(xhard);
% acheivable rate with soft decoding.
RSOFT = rate(pSOFT,poutSOFT,Pt,T);
% acheivable rate with hard decoding.
RHARD = rate(pHARD,poutHARD,Pt,T);
%-------------------------------------------------------------------------%
% PLOT RECEIVER OUTPUT
%-------------------------------------------------------------------------%
% plot 1000 first receiver outputs for each input.
if plotout
figure(99); hold all;
% plot axes
if showAxes
line(sqrt(snr)*[-10, 10],[0, 0],'color','k','linewidth',2);
line([0, 0],sqrt(snr)*[-10, 10],'color','k','linewidth',2);
end
% plot grid
if showGrid
for i = 1:length(edges)
line([edges(i) edges(i)], [edges(1), edges(end)],'color', [0.7 0.7 0.7],'linewidth',0.1);
line([edges(1), edges(end)],[edges(i) edges(i)],'color', [0.7 0.7 0.7],'linewidth',0.1);
end
end
% plot MRC output
if MAX <= 1000
for m = 1:M
plot(xsoft(m,:),'rx');
end
else
for m = 1:M
plot(xsoft(m,1:1000),'rx');
end
end
% plot input constellation
if showInpt
plot(sqrt(snr)*C,'ko','markersize',7);
end
axis equal; box on;
axis(max(edges)*[-1, 1, -1, 1]);
xlabel('In-phase','fontsize',20);
ylabel('Quadrature','fontsize',20);
title([num2str(M),sig,', ',rec,', SNR = ',num2str(snrdb),' dB, K = ',num2str(K),', Pk = ',num2str(Pt),', b = ',num2str(b)],'fontsize',16);
end
elapsedTime = toc;
fprintf('\tSimulation finished after %f seconds.\n',elapsedTime);
end
% FUNCTIONS
%-------------------------------------------------------------------------%
function R = quantise(Y,b)
% perform 1-bit quantisation if b = 1, otherwise return unquantised signal.
if b == 0
R = Y;
elseif b == 1
R = sign(real(Y)) + 1i*sign(imag(Y));
end
end
function [pSOFT, poutSOFT,edges] = softProb(xsoft, numbin)
% approximate soft channel law and output probability by mapping MRC/ZF
% output on a discrete grid and counting the number of instances that the
% MRC/ZF outputs occur in each grid region, for each input.
% extract number of inputs and number of random
% noise/interference/channel realisations
[M,MAX] = size(xsoft);
% find max/min-points of soft outputs to determine the the outer
% boundaries of the discretization grid.
min_edges = min(min(min(real(xsoft))), min(min(imag(xsoft))));
max_edges = max(max(max(real(xsoft))), max(max(imag(xsoft))));
bnd_edges = 1.1*max(abs(min_edges),abs(max_edges));
% edges that specify the grid boundaries per dimension
edges = linspace(-bnd_edges,bnd_edges+eps,numbin);
% channel law: count number of occurences in each bin
pSOFT = nan(M,length(edges)-1,length(edges)-1);
for i = 1:M
pSOFT(i,:,:) = histcn([real(xsoft(i,:))' ,imag(xsoft(i,:))'],edges,edges);
end
pSOFT = reshape(pSOFT,M,(length(edges)-1)^2)/MAX;
% output probability: average the channel law over all inputs
poutSOFT = mean(pSOFT);
end
function [pHARD, poutHARD] = hardProb(xhard)
% approximate hard MRC/ZF channel law and output probability by performing
% hard QAM/PSK decoiding and counting the number of instances that each
% hard output symbol occur, for each input symbol.
% extract number of inputs and number of random
% noise/interference/channel realisations
[M,MAX] = size(xhard);
% channel law
pHARD = nan(M,M);
for i = 1:1:M
for j = 1:1:M
pHARD(i,j) = sum(xhard(i,:) == j-1);
end
end
pHARD = pHARD/MAX;
% output probability
poutHARD = mean(pHARD);
end
function R = rate(p,pout,Pt,T)
% compute the acheivable rate by computing the mutual information and
% normaling with the number of transmitted pilots.
% extract number of inputs
M = size(p,1);
% mutual information
R = mean(sum(p .* log2(p ./ (repmat(pout,M,1)+eps) + eps),2));
% normalisation
R = (T-Pt)./T * R;
R = max(R,0);
end