source: ResearchApps/PHY/WARPLAB/WARPLab7/M_Code_Examples/wl_example_mimo_ofdm_txrx.m

Last change on this file was 4985, checked in by welsh, 7 years ago

updating so that commands are only issued to the appropriate node. made it easier to change RF interfaces.

File size: 31.8 KB
Line 
1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2% wl_example_mimo_ofdm_txrx.m
3% 2x2 MIMO OFDM Example
4% A detailed write-up of this example is available on the wiki:
5% http://warpproject.org/trac/wiki/WARPLab/Examples/MIMO_OFDM
6%
7% Copyright (c) 2015 Mango Communications - All Rights Reserved
8% Distributed under the WARP License (http://warpproject.org/license)
9%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
10clear
11
12% Params:
13USE_WARPLAB_TXRX        = 1;           % Enable WARPLab-in-the-loop (otherwise sim-only)
14WRITE_PNG_FILES         = 0;           % Enable writing plots to PNG
15CHANNEL                 = 11;          % Channel to tune Tx and Rx radios
16
17% Waveform params
18N_OFDM_SYMS             = 1000;        % Number of OFDM symbols (must be even valued)
19MOD_ORDER               = 16;          % Modulation order (2/4/16 = BSPK/QPSK/16-QAM)
20TX_SCALE                = 1.0;         % Scale for Tx waveform ([0:1])
21INTERP_RATE             = 2;           % Interpolation rate (must be 2)
22TX_SPATIAL_STREAM_SHIFT = 3;           % Number of samples to shift the transmission from RFB
23
24% OFDM params
25SC_IND_PILOTS           = [8 22 44 58];                           % Pilot subcarrier indices
26SC_IND_DATA             = [2:7 9:21 23:27 39:43 45:57 59:64];     % Data subcarrier indices
27N_SC                    = 64;                                     % Number of subcarriers
28CP_LEN                  = 16;                                     % Cyclic prefix length
29N_DATA_SYMS             = N_OFDM_SYMS * length(SC_IND_DATA);      % Number of data symbols (one per data-bearing subcarrier per OFDM symbol)
30
31% Rx processing params
32FFT_OFFSET                    = 4;           % Number of CP samples to use in FFT (on average)
33LTS_CORR_THRESH               = 0.8;         % Normalized threshold for LTS correlation
34DO_APPLY_CFO_CORRECTION       = 1;           % Enable CFO estimation/correction
35DO_APPLY_PHASE_ERR_CORRECTION = 1;           % Enable Residual CFO estimation/correction
36DO_APPLY_SFO_CORRECTION       = 1;           % Enable SFO estimation/correction
37DECIMATE_RATE                 = INTERP_RATE;
38
39% WARPLab experiment params
40USE_AGC                 = true;        % Use the AGC if running on WARP hardware
41MAX_TX_LEN              = 2^19;        % Maximum number of samples to use for this experiment
42SAMP_PADDING            = 100;         % Extra samples to receive to ensure both start and end of waveform visible
43
44
45if(USE_WARPLAB_TXRX)
46    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
47    % Set up the WARPLab experiment
48    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
49
50    NUMNODES = 2;
51
52    % Create a vector of node objects
53    nodes   = wl_initNodes(NUMNODES);
54    node_tx = nodes(1);
55    node_rx = nodes(2);
56
57    % Create a UDP broadcast trigger and tell each node to be ready for it
58    eth_trig = wl_trigger_eth_udp_broadcast;
59    wl_triggerManagerCmd(nodes, 'add_ethernet_trigger', [eth_trig]);
60
61    % Read Trigger IDs into workspace
62    trig_in_ids  = wl_getTriggerInputIDs(nodes(1));
63    trig_out_ids = wl_getTriggerOutputIDs(nodes(1));
64
65    % For both nodes, we will allow Ethernet to trigger the buffer baseband and the AGC
66    wl_triggerManagerCmd(nodes, 'output_config_input_selection', [trig_out_ids.BASEBAND, trig_out_ids.AGC], [trig_in_ids.ETH_A]);
67
68    % Set the trigger output delays.
69    nodes.wl_triggerManagerCmd('output_config_delay', [trig_out_ids.BASEBAND], 0);
70    nodes.wl_triggerManagerCmd('output_config_delay', [trig_out_ids.AGC], 3000);     %3000 ns delay before starting the AGC
71
72    % Get IDs for the interfaces on the boards.
73    ifc_ids_TX = wl_getInterfaceIDs(node_tx);
74    ifc_ids_RX = wl_getInterfaceIDs(node_rx);
75
76    % Set up the TX / RX nodes and RF interfaces
77    TX_RF     = ifc_ids_TX.RF_ON_BOARD;
78    TX_RF_VEC = ifc_ids_TX.RF_ON_BOARD_VEC;
79    TX_RF_ALL = ifc_ids_TX.RF_ALL;
80   
81    RX_RF     = ifc_ids_RX.RF_ON_BOARD;
82    RX_RF_VEC = ifc_ids_RX.RF_ON_BOARD_VEC;
83    RX_RF_ALL = ifc_ids_RX.RF_ALL;
84
85    % Set up the interface for the experiment
86    wl_interfaceCmd(node_tx, TX_RF_ALL, 'channel', 2.4, CHANNEL);
87    wl_interfaceCmd(node_rx, RX_RF_ALL, 'channel', 2.4, CHANNEL);
88
89    wl_interfaceCmd(node_tx, TX_RF_ALL, 'tx_gains', 3, 30);
90   
91    if(USE_AGC)
92        wl_interfaceCmd(node_rx, RX_RF_ALL, 'rx_gain_mode', 'automatic');
93        wl_basebandCmd(node_rx, 'agc_target', -13); %-13
94    else
95        wl_interfaceCmd(node_rx, RX_RF_ALL, 'rx_gain_mode', 'manual');
96        RxGainRF = 2;                  % Rx RF Gain in [1:3]
97        RxGainBB = 12;                 % Rx Baseband Gain in [0:31]
98        wl_interfaceCmd(node_rx, RX_RF_ALL, 'rx_gains', RxGainRF, RxGainBB);
99    end
100
101    % Get parameters from the node
102    SAMP_FREQ    = wl_basebandCmd(nodes(1), 'tx_buff_clk_freq');
103    Ts           = 1/SAMP_FREQ;
104
105    % We will read the transmitter's maximum I/Q buffer length
106    % and assign that value to a temporary variable.
107    %
108    % NOTE:  We assume that the buffers sizes are the same for all interfaces
109
110    maximum_buffer_len = min(MAX_TX_LEN, wl_basebandCmd(node_tx, TX_RF_VEC, 'tx_buff_max_num_samples'));
111    example_mode_string = 'hw';
112else
113    % Use sane defaults for hardware-dependent params in sim-only version
114    maximum_buffer_len  = min(MAX_TX_LEN, 2^20);
115    SAMP_FREQ           = 40e6;
116    example_mode_string = 'sim';
117end
118
119%% Define a half-band 2x interpolation filter response
120interp_filt2 = zeros(1,43);
121interp_filt2([1 3 5 7 9 11 13 15 17 19 21]) = [12 -32 72 -140 252 -422 682 -1086 1778 -3284 10364];
122interp_filt2([23 25 27 29 31 33 35 37 39 41 43]) = interp_filt2(fliplr([1 3 5 7 9 11 13 15 17 19 21]));
123interp_filt2(22) = 16384;
124interp_filt2 = interp_filt2./max(abs(interp_filt2));
125
126% Define the preamble
127% Note: The STS symbols in the preamble meet the requirements needed by the
128% AGC core at the receiver. Details on the operation of the AGC are
129% available on the wiki: http://warpproject.org/trac/wiki/WARPLab/AGC
130sts_f = zeros(1,64);
131sts_f(1:27) = [0 0 0 0 -1-1i 0 0 0 -1-1i 0 0 0 1+1i 0 0 0 1+1i 0 0 0 1+1i 0 0 0 1+1i 0 0];
132sts_f(39:64) = [0 0 1+1i 0 0 0 -1-1i 0 0 0 1+1i 0 0 0 -1-1i 0 0 0 -1-1i 0 0 0 1+1i 0 0 0];
133sts_t = ifft(sqrt(13/6).*sts_f, 64);
134sts_t = sts_t(1:16);
135
136% LTS for CFO and channel estimation
137lts_f = [0 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 0 0 0 0 0 0 0 0 0 0 0 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];
138lts_t = ifft(lts_f, 64);
139
140% We break the construction of our preamble into two pieces. First, the
141% legacy portion, is used for CFO recovery and timing synchronization at
142% the receiver. The processing of this portion of the preamble is SISO.
143% Second, we include explicit MIMO channel training symbols.
144
145% Legacy Preamble
146
147% Use 30 copies of the 16-sample STS for extra AGC settling margin
148% To avoid accidentally beamforming the preamble transmissions, we will
149% let RFA be dominant and handle the STS and first set of LTS. We will
150% append an extra LTS sequence from RFB so that we can build out the
151% channel matrix at the receiver
152
153sts_t_rep = repmat(sts_t, 1, 30);
154
155preamble_legacy_A = [sts_t_rep, lts_t(33:64), lts_t, lts_t];
156preamble_legacy_B = [circshift(sts_t_rep, [0, TX_SPATIAL_STREAM_SHIFT]), zeros(1, 160)];
157
158% MIMO Preamble
159
160% There are many strategies for training MIMO channels. Here, we will use
161% the LTS sequence defined before and orthogonalize over time. First we
162% will send the sequence on stream A and then we will send it on stream B
163
164preamble_mimo_A = [lts_t(33:64), lts_t, zeros(1,96)];
165preamble_mimo_B = [zeros(1,96), lts_t(33:64), lts_t];
166
167preamble_A = [preamble_legacy_A, preamble_mimo_A];
168preamble_B = [preamble_legacy_B, preamble_mimo_B];
169
170% Sanity check variables that affect the number of Tx samples
171if(SAMP_PADDING + INTERP_RATE*((N_OFDM_SYMS/2 * (N_SC + CP_LEN)) + length(preamble_A) + 100) > maximum_buffer_len)
172    fprintf('Too many OFDM symbols for TX_NUM_SAMPS!\n');
173    fprintf('Raise TX_NUM_SAMPS to %d, or \n', SAMP_PADDING + INTERP_RATE*((N_OFDM_SYMS/2 * (N_SC + CP_LEN)) + length(preamble_A) + 100));
174    fprintf('Reduce N_OFDM_SYMS to %d\n',  2*(floor(( (maximum_buffer_len/INTERP_RATE)-length(preamble_A)-100-SAMP_PADDING )/( N_SC + CP_LEN )) - 1));
175    return;
176end
177
178%% Generate a payload of random integers
179tx_data = randi(MOD_ORDER, 1, N_DATA_SYMS) - 1;
180
181% Functions for data -> complex symbol mapping (like qammod, avoids comm toolbox requirement)
182% These anonymous functions implement the modulation mapping from IEEE 802.11-2012 Section 18.3.5.8
183modvec_bpsk   =  (1/sqrt(2))  .* [-1 1];
184modvec_16qam  = (1/sqrt(10))  .* [-3 -1 +3 +1];
185
186mod_fcn_bpsk  = @(x) complex(modvec_bpsk(1+x),0);
187mod_fcn_qpsk  = @(x) complex(modvec_bpsk(1+bitshift(x, -1)), modvec_bpsk(1+mod(x, 2)));
188mod_fcn_16qam = @(x) complex(modvec_16qam(1+bitshift(x, -2)), modvec_16qam(1+mod(x,4)));
189
190% Map the data values on to complex symbols
191switch MOD_ORDER
192    case 2         % BPSK
193        tx_syms = arrayfun(mod_fcn_bpsk, tx_data);
194    case 4         % QPSK
195        tx_syms = arrayfun(mod_fcn_qpsk, tx_data);
196    case 16        % 16-QAM
197        tx_syms = arrayfun(mod_fcn_16qam, tx_data);
198    otherwise
199        fprintf('Invalid MOD_ORDER (%d)!  Must be in [2, 4, 16]\n', MOD_ORDER);
200        return;
201end
202
203% Reshape the symbol vector into two different spatial streams
204tx_syms_space_mat = reshape(tx_syms, 2, length(tx_syms)/2);
205
206% Break up the matrix into a vector for each antenna
207tx_syms_A = tx_syms_space_mat(1,:);
208tx_syms_B = tx_syms_space_mat(2,:);
209
210
211% Reshape the symbol vector to a matrix with one column per OFDM symbol
212tx_syms_mat_A = reshape(tx_syms_A, length(SC_IND_DATA), N_OFDM_SYMS/2);
213tx_syms_mat_B = reshape(tx_syms_B, length(SC_IND_DATA), N_OFDM_SYMS/2);
214
215% Define the pilot tone values as BPSK symbols
216%  We will transmit pilots only on RF A
217pilots_A = [1 1 -1 1].';
218pilots_B = [0 0 0 0].';
219
220% Repeat the pilots across all OFDM symbols
221pilots_mat_A = repmat(pilots_A, 1, N_OFDM_SYMS/2);
222pilots_mat_B = repmat(pilots_B, 1, N_OFDM_SYMS/2);
223
224%% IFFT
225
226% Construct the IFFT input matrix
227ifft_in_mat_A = zeros(N_SC, N_OFDM_SYMS/2);
228ifft_in_mat_B = zeros(N_SC, N_OFDM_SYMS/2);
229
230% Insert the data and pilot values; other subcarriers will remain at 0
231ifft_in_mat_A(SC_IND_DATA, :)   = tx_syms_mat_A;
232ifft_in_mat_A(SC_IND_PILOTS, :) = pilots_mat_A;
233
234ifft_in_mat_B(SC_IND_DATA, :)   = tx_syms_mat_B;
235ifft_in_mat_B(SC_IND_PILOTS, :) = pilots_mat_B;
236
237%Perform the IFFT
238tx_payload_mat_A = ifft(ifft_in_mat_A, N_SC, 1);
239tx_payload_mat_B = ifft(ifft_in_mat_B, N_SC, 1);
240
241% Insert the cyclic prefix
242if(CP_LEN > 0)
243    tx_cp = tx_payload_mat_A((end-CP_LEN+1 : end), :);
244    tx_payload_mat_A = [tx_cp; tx_payload_mat_A];
245
246    tx_cp = tx_payload_mat_B((end-CP_LEN+1 : end), :);
247    tx_payload_mat_B = [tx_cp; tx_payload_mat_B];
248end
249
250% Reshape to a vector
251tx_payload_vec_A = reshape(tx_payload_mat_A, 1, numel(tx_payload_mat_A));
252tx_payload_vec_B = reshape(tx_payload_mat_B, 1, numel(tx_payload_mat_B));
253
254% Construct the full time-domain OFDM waveform
255tx_vec_A = [preamble_A tx_payload_vec_A];
256tx_vec_B = [preamble_B tx_payload_vec_B];
257
258% Pad with zeros for transmission
259tx_vec_padded_A = [tx_vec_A zeros(1,50)];
260tx_vec_padded_B = [tx_vec_B zeros(1,50)];
261
262%% Interpolate
263if(INTERP_RATE == 1)
264    tx_vec_air_A = tx_vec_padded_A;
265    tx_vec_air_B = tx_vec_padded_B;
266elseif(INTERP_RATE == 2)
267    % Zero pad then filter (same as interp or upfirdn without signal processing toolbox)
268    tx_vec_2x_A = zeros(1, 2*numel(tx_vec_padded_A));
269    tx_vec_2x_A(1:2:end) = tx_vec_padded_A;
270    tx_vec_air_A = filter(interp_filt2, 1, tx_vec_2x_A);
271    tx_vec_2x_B = zeros(1, 2*numel(tx_vec_padded_B));
272    tx_vec_2x_B(1:2:end) = tx_vec_padded_B;
273    tx_vec_air_B = filter(interp_filt2, 1, tx_vec_2x_B);
274end
275
276% Scale the Tx vector to +/- 1
277tx_vec_air_A = TX_SCALE .* tx_vec_air_A ./ max(abs(tx_vec_air_A));
278tx_vec_air_B = TX_SCALE .* tx_vec_air_B ./ max(abs(tx_vec_air_B));
279
280
281TX_NUM_SAMPS = length(tx_vec_air_A);
282
283if(USE_WARPLAB_TXRX)
284    wl_basebandCmd(nodes, 'tx_delay', 0);
285    wl_basebandCmd(nodes, 'tx_length', TX_NUM_SAMPS+100);                   % Number of samples to send
286    wl_basebandCmd(nodes, 'rx_length', TX_NUM_SAMPS+SAMP_PADDING);      % Number of samples to receive
287end
288
289if(USE_WARPLAB_TXRX)
290    %% WARPLab Tx/Rx
291
292    tx_mat_air = [tx_vec_air_A(:) , tx_vec_air_B(:)];
293    %tx_mat_air = [tx_vec_air_B(:) , tx_vec_air_A(:)];
294
295    % Write the Tx waveform to the Tx node
296    wl_basebandCmd(node_tx, TX_RF_VEC, 'write_IQ', tx_mat_air);
297
298    % Enable the Tx and Rx radios
299    wl_interfaceCmd(node_tx, TX_RF, 'tx_en');
300    wl_interfaceCmd(node_rx, RX_RF, 'rx_en');
301
302    % Enable the Tx and Rx buffers
303    wl_basebandCmd(node_tx, TX_RF, 'tx_buff_en');
304    wl_basebandCmd(node_rx, RX_RF, 'rx_buff_en');
305
306    % Trigger the Tx/Rx cycle at both nodes
307    eth_trig.send();
308
309    % Retrieve the received waveform from the Rx node
310    rx_mat_air = wl_basebandCmd(node_rx, RX_RF_VEC, 'read_IQ', 0, TX_NUM_SAMPS+SAMP_PADDING);
311
312    rx_vec_air_A = rx_mat_air(:,1).';
313    rx_vec_air_B = rx_mat_air(:,2).';
314
315    % Disable the Tx/Rx radios and buffers
316    wl_basebandCmd(node_tx, TX_RF_ALL, 'tx_rx_buff_dis');
317    wl_basebandCmd(node_rx, RX_RF_ALL, 'tx_rx_buff_dis');
318   
319    wl_interfaceCmd(node_tx, TX_RF_ALL, 'tx_rx_dis');
320    wl_interfaceCmd(node_rx, RX_RF_ALL, 'tx_rx_dis');
321else
322
323    rx_vec_air_A = tx_vec_air_A + .5 * tx_vec_air_B;
324    rx_vec_air_B = tx_vec_air_B + .5 * tx_vec_air_A;
325
326    rx_vec_air_A = [rx_vec_air_A, zeros(1,SAMP_PADDING)];
327    rx_vec_air_B = [rx_vec_air_B, zeros(1,SAMP_PADDING)];
328
329    rx_vec_air_A = rx_vec_air_A + 1e-2*complex(randn(1,length(rx_vec_air_A)), randn(1,length(rx_vec_air_A)));
330    rx_vec_air_B = rx_vec_air_B + 1e-2*complex(randn(1,length(rx_vec_air_B)), randn(1,length(rx_vec_air_B)));
331
332end
333
334%% Decimate
335if(DECIMATE_RATE == 1)
336    raw_rx_dec_A = rx_vec_air_A;
337    raw_rx_dec_B = rx_vec_air_B;
338elseif(DECIMATE_RATE == 2)
339    raw_rx_dec_A = filter(interp_filt2, 1, rx_vec_air_A);
340    raw_rx_dec_A = raw_rx_dec_A(1:2:end);
341    raw_rx_dec_B = filter(interp_filt2, 1, rx_vec_air_B);
342    raw_rx_dec_B = raw_rx_dec_B(1:2:end);
343end
344
345%% Correlate for LTS
346
347% For simplicity, we'll only use RFA for LTS correlation and peak
348% discovery. A straightforward addition would be to repeat this process for
349% RFB and combine the results for detection diversity.
350
351% Complex cross correlation of Rx waveform with time-domain LTS
352lts_corr = abs(conv(conj(fliplr(lts_t)), sign(raw_rx_dec_A)));
353
354% Skip early and late samples - avoids occasional false positives from pre-AGC samples
355lts_corr = lts_corr(32:end-32);
356
357% Find all correlation peaks
358lts_peaks = find(lts_corr > LTS_CORR_THRESH*max(lts_corr));
359
360% Select best candidate correlation peak as LTS-payload boundary
361% In this MIMO example, we actually have 3 LTS symbols sent in a row.
362% The first two are sent by RFA on the TX node and the last one was sent
363% by RFB. We will actually look for the separation between the first and the
364% last for synchronizing our starting index.
365
366[LTS1, LTS2] = meshgrid(lts_peaks,lts_peaks);
367[lts_last_peak_index,y] = find(LTS2-LTS1 == length(lts_t));
368
369% Stop if no valid correlation peak was found
370if(isempty(lts_last_peak_index))
371    fprintf('No LTS Correlation Peaks Found!\n');
372    return;
373end
374
375% Set the sample indices of the payload symbols and preamble
376% The "+32" here corresponds to the 32-sample cyclic prefix on the preamble LTS
377% The "+192" corresponds to the length of the extra training symbols for MIMO channel estimation
378mimo_training_ind = lts_peaks(max(lts_last_peak_index)) + 32;
379payload_ind = mimo_training_ind + 192;
380
381% Subtract of 2 full LTS sequences and one cyclic prefixes
382% The "-160" corresponds to the length of the preamble LTS (2.5 copies of 64-sample LTS)
383lts_ind = mimo_training_ind-160;
384
385if(DO_APPLY_CFO_CORRECTION)
386    %Extract LTS (not yet CFO corrected)
387    rx_lts = raw_rx_dec_A(lts_ind : lts_ind+159); %Extract the first two LTS for CFO
388    rx_lts1 = rx_lts(-64+-FFT_OFFSET + [97:160]);
389    rx_lts2 = rx_lts(-FFT_OFFSET + [97:160]);
390
391    %Calculate coarse CFO est
392    rx_cfo_est_lts = mean(unwrap(angle(rx_lts2 .* conj(rx_lts1))));
393    rx_cfo_est_lts = rx_cfo_est_lts/(2*pi*64);
394else
395    rx_cfo_est_lts = 0;
396end
397
398% Apply CFO correction to raw Rx waveforms
399rx_cfo_corr_t = exp(-1i*2*pi*rx_cfo_est_lts*[0:length(raw_rx_dec_A)-1]);
400rx_dec_cfo_corr_A = raw_rx_dec_A .* rx_cfo_corr_t;
401rx_dec_cfo_corr_B = raw_rx_dec_B .* rx_cfo_corr_t;
402
403% MIMO Channel Estimatation
404lts_ind_TXA_start = mimo_training_ind + 32 - FFT_OFFSET;
405lts_ind_TXA_end = lts_ind_TXA_start + 64 - 1;
406
407lts_ind_TXB_start = mimo_training_ind + 32 + 64 + 32 - FFT_OFFSET;
408lts_ind_TXB_end = lts_ind_TXB_start + 64 - 1;
409
410rx_lts_AA = rx_dec_cfo_corr_A( lts_ind_TXA_start:lts_ind_TXA_end );
411rx_lts_BA = rx_dec_cfo_corr_A( lts_ind_TXB_start:lts_ind_TXB_end );
412
413rx_lts_AB = rx_dec_cfo_corr_B( lts_ind_TXA_start:lts_ind_TXA_end );
414rx_lts_BB = rx_dec_cfo_corr_B( lts_ind_TXB_start:lts_ind_TXB_end );
415
416rx_lts_AA_f = fft(rx_lts_AA);
417rx_lts_BA_f = fft(rx_lts_BA);
418
419rx_lts_AB_f = fft(rx_lts_AB);
420rx_lts_BB_f = fft(rx_lts_BB);
421
422% Calculate channel estimate
423rx_H_est_AA = lts_f .* rx_lts_AA_f;
424rx_H_est_BA = lts_f .* rx_lts_BA_f;
425
426rx_H_est_AB = lts_f .* rx_lts_AB_f;
427rx_H_est_BB = lts_f .* rx_lts_BB_f;
428
429%% Rx payload processing
430
431% Extract the payload samples (integral number of OFDM symbols following preamble)
432payload_vec_A = rx_dec_cfo_corr_A(payload_ind : payload_ind+(N_OFDM_SYMS/2)*(N_SC+CP_LEN)-1);
433payload_mat_A = reshape(payload_vec_A, (N_SC+CP_LEN), (N_OFDM_SYMS/2));
434
435payload_vec_B = rx_dec_cfo_corr_B(payload_ind : payload_ind+(N_OFDM_SYMS/2)*(N_SC+CP_LEN)-1);
436payload_mat_B = reshape(payload_vec_B, (N_SC+CP_LEN), (N_OFDM_SYMS/2));
437
438% Remove the cyclic prefix, keeping FFT_OFFSET samples of CP (on average)
439payload_mat_noCP_A = payload_mat_A(CP_LEN-FFT_OFFSET+[1:N_SC], :);
440payload_mat_noCP_B = payload_mat_B(CP_LEN-FFT_OFFSET+[1:N_SC], :);
441
442% Take the FFT
443syms_f_mat_A = fft(payload_mat_noCP_A, N_SC, 1);
444syms_f_mat_B = fft(payload_mat_noCP_B, N_SC, 1);
445
446% Equalize pilots
447% Because we only used Tx RFA to send pilots, we can do SISO equalization
448% here. This is zero-forcing (just divide by chan estimates)
449syms_eq_mat_pilots = syms_f_mat_A ./ repmat(rx_H_est_AA.', 1, N_OFDM_SYMS/2);
450
451if DO_APPLY_SFO_CORRECTION
452    % SFO manifests as a frequency-dependent phase whose slope increases
453    % over time as the Tx and Rx sample streams drift apart from one
454    % another. To correct for this effect, we calculate this phase slope at
455    % each OFDM symbol using the pilot tones and use this slope to
456    % interpolate a phase correction for each data-bearing subcarrier.
457
458    % Extract the pilot tones and "equalize" them by their nominal Tx values
459    pilots_f_mat = syms_eq_mat_pilots(SC_IND_PILOTS,:);
460    pilots_f_mat_comp = pilots_f_mat.*pilots_mat_A;
461
462    % Calculate the phases of every Rx pilot tone
463    pilot_phases = unwrap(angle(fftshift(pilots_f_mat_comp, 1)), [], 1);
464
465    pilot_spacing_mat = repmat(mod(diff(fftshift(SC_IND_PILOTS)), 64).', 1, N_OFDM_SYMS/2);
466    pilot_slope_mat = mean(diff(pilot_phases) ./ pilot_spacing_mat);
467
468    % Calculate the SFO correction phases for each OFDM symbol
469    pilot_phase_sfo_corr = fftshift((-32:31).' * pilot_slope_mat, 1);
470    pilot_phase_corr = exp(-1i*(pilot_phase_sfo_corr));
471
472    % Apply the pilot phase correction per symbol
473    syms_f_mat_A = syms_f_mat_A .* pilot_phase_corr;
474    syms_f_mat_B = syms_f_mat_B .* pilot_phase_corr;
475else
476    % Define an empty SFO correction matrix (used by plotting code below)
477    pilot_phase_sfo_corr = zeros(N_SC, N_OFDM_SYMS);
478end
479
480% Extract the pilots and calculate per-symbol phase error
481if DO_APPLY_PHASE_ERR_CORRECTION
482    pilots_f_mat = syms_eq_mat_pilots(SC_IND_PILOTS, :);
483    pilot_phase_err = angle(mean(pilots_f_mat.*pilots_mat_A));
484else
485    % Define an empty phase correction vector (used by plotting code below)
486    pilot_phase_err = zeros(1, N_OFDM_SYMS/2);
487end
488pilot_phase_corr = repmat(exp(-1i*pilot_phase_err), N_SC, 1);
489
490% Apply pilot phase correction to both received streams
491syms_f_mat_pc_A = syms_f_mat_A .* pilot_phase_corr;
492syms_f_mat_pc_B = syms_f_mat_B .* pilot_phase_corr;
493
494% MIMO Equalization
495% We need to apply the MIMO equalization to each subcarrier separately.
496% There, unfortunately, is no great vector-y solution to do this, so we
497% reluctantly employ a FOR loop.
498
499syms_eq_mat_A = zeros(N_SC, N_OFDM_SYMS/2);
500syms_eq_mat_B = zeros(N_SC, N_OFDM_SYMS/2);
501channel_condition_mat = zeros(1,N_SC);
502for sc_idx = [SC_IND_DATA, SC_IND_PILOTS]
503   y = [syms_f_mat_pc_A(sc_idx,:) ; syms_f_mat_pc_B(sc_idx,:)];
504   H = [rx_H_est_AA(sc_idx), rx_H_est_BA(sc_idx) ; rx_H_est_AB(sc_idx), rx_H_est_BB(sc_idx)];
505   x = inv(H)*y;
506   syms_eq_mat_A(sc_idx, :) = x(1,:);
507   syms_eq_mat_B(sc_idx, :) = x(2,:);
508   channel_condition_mat(sc_idx) = rcond(H);
509end
510
511%subplot(2,1,1)
512%plot(syms_eq_mat_A,'.')
513%subplot(2,1,2)
514%plot(syms_eq_mat_B,'.')
515
516payload_syms_mat_A = syms_eq_mat_A(SC_IND_DATA, :);
517payload_syms_mat_B = syms_eq_mat_B(SC_IND_DATA, :);
518
519%% Demodulate
520rx_syms_A = reshape(payload_syms_mat_A, 1, N_DATA_SYMS/2);
521rx_syms_B = reshape(payload_syms_mat_B, 1, N_DATA_SYMS/2);
522
523% Combine both streams to a single vector of symbols
524rx_syms_space_mat = [rx_syms_A; rx_syms_B];
525rx_syms = reshape(rx_syms_space_mat, 1, length(rx_syms_A)*2);
526
527demod_fcn_bpsk = @(x) double(real(x)>0);
528demod_fcn_qpsk = @(x) double(2*(real(x)>0) + 1*(imag(x)>0));
529demod_fcn_16qam = @(x) (8*(real(x)>0)) + (4*(abs(real(x))<0.6325)) + (2*(imag(x)>0)) + (1*(abs(imag(x))<0.6325));
530
531switch(MOD_ORDER)
532    case 2         % BPSK
533        rx_data = arrayfun(demod_fcn_bpsk, rx_syms);
534    case 4         % QPSK
535        rx_data = arrayfun(demod_fcn_qpsk, rx_syms);
536    case 16        % 16-QAM
537        rx_data = arrayfun(demod_fcn_16qam, rx_syms);
538end
539
540
541%% Plot Results
542cf = 0;
543
544% Tx signal
545cf = cf + 1;
546figure(cf); clf;
547
548subplot(2,2,1);
549plot(real(tx_vec_air_A), 'b');
550axis([0 length(tx_vec_air_A) -TX_SCALE TX_SCALE])
551grid on;
552title('RFA Tx Waveform (I)');
553
554subplot(2,2,2);
555plot(imag(tx_vec_air_A), 'r');
556axis([0 length(tx_vec_air_A) -TX_SCALE TX_SCALE])
557grid on;
558title('RFA Tx Waveform (Q)');
559
560subplot(2,2,3);
561plot(real(tx_vec_air_B), 'b');
562axis([0 length(tx_vec_air_B) -TX_SCALE TX_SCALE])
563grid on;
564title('RFB Tx Waveform (I)');
565
566subplot(2,2,4);
567plot(imag(tx_vec_air_B), 'r');
568axis([0 length(tx_vec_air_B) -TX_SCALE TX_SCALE])
569grid on;
570title('RFB Tx Waveform (Q)');
571
572if(WRITE_PNG_FILES)
573    print(gcf,sprintf('wl_mimo_ofdm_plots_%s_txIQ', example_mode_string), '-dpng', '-r96', '-painters')
574end
575
576% Rx signal
577cf = cf + 1;
578figure(cf); clf;
579subplot(2,2,1);
580plot(real(rx_vec_air_A), 'b');
581axis([0 length(rx_vec_air_A) -TX_SCALE TX_SCALE])
582grid on;
583title('RFA Rx Waveform (I)');
584
585subplot(2,2,2);
586plot(imag(rx_vec_air_A), 'r');
587axis([0 length(rx_vec_air_A) -TX_SCALE TX_SCALE])
588grid on;
589title('RFA Rx Waveform (Q)');
590
591subplot(2,2,3);
592plot(real(rx_vec_air_B), 'b');
593axis([0 length(rx_vec_air_B) -TX_SCALE TX_SCALE])
594grid on;
595title('RFB Rx Waveform (I)');
596
597subplot(2,2,4);
598plot(imag(rx_vec_air_B), 'r');
599axis([0 length(rx_vec_air_B) -TX_SCALE TX_SCALE])
600grid on;
601title('RFB Rx Waveform (Q)');
602
603if(WRITE_PNG_FILES)
604    print(gcf,sprintf('wl_mimo_ofdm_plots_%s_rxIQ', example_mode_string), '-dpng', '-r96', '-painters')
605end
606
607% Rx LTS correlation
608cf = cf + 1;
609figure(cf); clf;
610lts_to_plot = lts_corr;
611plot(lts_to_plot, '.-b', 'LineWidth', 1);
612hold on;
613grid on;
614line([1 length(lts_to_plot)], LTS_CORR_THRESH*max(lts_to_plot)*[1 1], 'LineStyle', '--', 'Color', 'r', 'LineWidth', 2);
615title('LTS Correlation and Threshold')
616xlabel('Sample Index')
617myAxis = axis();
618axis([1, 1000, myAxis(3), myAxis(4)])
619
620if(WRITE_PNG_FILES)
621    print(gcf,sprintf('wl_mimo_ofdm_plots_%s_ltsCorr', example_mode_string), '-dpng', '-r96', '-painters')
622end
623
624% Channel Estimates
625cf = cf + 1;
626
627rx_H_est_plot_AA = rx_H_est_AA;
628rx_H_est_plot_AB = rx_H_est_AB;
629rx_H_est_plot_BA = rx_H_est_BA;
630rx_H_est_plot_BB = rx_H_est_BB;
631
632x = (20/N_SC) * (-(N_SC/2):(N_SC/2 - 1));
633
634figure(cf); clf;
635
636subplot(3,2,1);
637bh = bar(x, fftshift(abs(rx_H_est_plot_AA)),1,'LineWidth', 1);
638shading flat
639set(bh,'FaceColor',[0 0 1])
640axis([min(x) max(x) 0 1.1*max(abs(rx_H_est_plot_AA))])
641grid on;
642title('A->A Channel Estimates (Magnitude)')
643xlabel('Baseband Frequency (MHz)')
644
645subplot(3,2,2);
646bh = bar(x, fftshift(abs(rx_H_est_plot_AB)),1,'LineWidth', 1);
647shading flat
648set(bh,'FaceColor',[0 0 1])
649axis([min(x) max(x) 0 1.1*max(abs(rx_H_est_plot_AB))])
650grid on;
651title('A->B Channel Estimates (Magnitude)')
652xlabel('Baseband Frequency (MHz)')
653
654subplot(3,2,3);
655bh = bar(x, fftshift(abs(rx_H_est_plot_BA)),1,'LineWidth', 1);
656shading flat
657set(bh,'FaceColor',[0 0 1])
658axis([min(x) max(x) 0 1.1*max(abs(rx_H_est_plot_BA))])
659grid on;
660title('B->A Channel Estimates (Magnitude)')
661xlabel('Baseband Frequency (MHz)')
662
663subplot(3,2,4);
664bh = bar(x, fftshift(abs(rx_H_est_plot_BB)),1,'LineWidth', 1);
665shading flat
666set(bh,'FaceColor',[0 0 1])
667axis([min(x) max(x) 0 1.1*max(abs(rx_H_est_plot_BB))])
668grid on;
669title('B->B Channel Estimates (Magnitude)')
670xlabel('Baseband Frequency (MHz)')
671
672
673subplot(3,1,3);
674bh = bar(x, fftshift(channel_condition_mat) ,1,'LineWidth', 1);
675shading flat
676set(bh,'FaceColor',[1 0 0])
677axis([min(x) max(x) 0 1.1])
678grid on;
679title('Channel Condition')
680xlabel('Baseband Frequency (MHz)')
681
682if(WRITE_PNG_FILES)
683    print(gcf,sprintf('wl_mimo_ofdm_plots_%s_chanEst', example_mode_string), '-dpng', '-r96', '-painters')
684end
685
686% Pilot phase error estimate
687cf = cf + 1;
688figure(cf); clf;
689subplot(2,1,1)
690plot(pilot_phase_err, 'b', 'LineWidth', 2);
691title('Phase Error Estimates')
692xlabel('OFDM Symbol Index')
693ylabel('Radians')
694axis([1 N_OFDM_SYMS/2 -3.2 3.2])
695grid on
696
697h = colorbar;
698set(h,'Visible','off');
699
700subplot(2,1,2)
701imagesc(1:N_OFDM_SYMS, (SC_IND_DATA - N_SC/2), fftshift(pilot_phase_sfo_corr,1))
702xlabel('OFDM Symbol Index')
703ylabel('Subcarrier Index')
704title('Phase Correction for SFO')
705colorbar
706myAxis = caxis();
707if(myAxis(2)-myAxis(1) < (pi))
708   caxis([-pi/2 pi/2])
709end
710
711if(WRITE_PNG_FILES)
712    print(gcf,sprintf('wl_mimo_ofdm_plots_%s_phaseError', example_mode_string), '-dpng', '-r96', '-painters')
713end
714
715% Symbol constellation
716cf = cf + 1;
717figure(cf); clf;
718
719subplot(2,2,1)
720rx_syms_plot = syms_f_mat_pc_A(SC_IND_DATA,:);
721plot(rx_syms_plot(:),'go','MarkerSize',1)
722axis square
723grid on
724title('Unequalized Rx Symbols RFA')
725myAxis = [];
726myAxis(1,:) = axis();
727
728subplot(2,2,2)
729rx_syms_plot = syms_f_mat_pc_B(SC_IND_DATA,:);
730plot(rx_syms_plot(:),'go','MarkerSize',1)
731axis square
732grid on
733title('Unequalized Rx Symbols RFB')
734myAxis(2,:) = axis();
735
736subplot(2,2,1); axis([-max(abs(myAxis(:))), max(abs(myAxis(:))), ...
737                      -max(abs(myAxis(:))), max(abs(myAxis(:)))])
738subplot(2,2,2); axis([-max(abs(myAxis(:))), max(abs(myAxis(:))), ...
739                      -max(abs(myAxis(:))), max(abs(myAxis(:)))])
740
741
742
743subplot(2,2,3)
744rx_syms_plot = syms_eq_mat_A(SC_IND_DATA,:);
745plot(rx_syms_plot(:),'ro','MarkerSize',1)
746axis square; axis(1.5*[-1 1 -1 1]);
747grid on
748title('Equalized Rx Symbols RFA')
749hold on;
750plot(tx_syms_mat_A(:),'bo');
751hold off;
752
753subplot(2,2,4)
754rx_syms_plot = syms_eq_mat_B(SC_IND_DATA,:);
755plot(rx_syms_plot(:),'ro','MarkerSize',1)
756axis square; axis(1.5*[-1 1 -1 1]);
757grid on
758title('Equalized Rx Symbols RFB')
759hold on;
760plot(tx_syms_mat_B(:),'bo');
761hold off;
762
763if(WRITE_PNG_FILES)
764    print(gcf,sprintf('wl_mimo_ofdm_plots_%s_constellations', example_mode_string), '-dpng', '-r96', '-painters')
765end
766
767
768
769% EVM & SNR
770cf = cf + 1;
771figure(cf); clf;
772
773evm_mat_A = abs(payload_syms_mat_A - tx_syms_mat_A).^2;
774evm_mat_B = abs(payload_syms_mat_B - tx_syms_mat_B).^2;
775aevms_A = mean(evm_mat_A(:));
776aevms_B = mean(evm_mat_B(:));
777snr_A = 10*log10(1./aevms_A);
778snr_B = 10*log10(1./aevms_B);
779
780subplot(2,2,1)
781plot(100*evm_mat_A(:),'o','MarkerSize',1)
782axis tight
783hold on
784plot([1 length(evm_mat_A(:))], 100*[aevms_A, aevms_A],'r','LineWidth',4)
785hold off
786xlabel('Data Symbol Index')
787ylabel('EVM (%)');
788legend('Per-Symbol EVM','Average EVM','Location','NorthWest');
789
790title('Stream A')
791grid on
792myAxis_mat(1,:) = axis();
793
794subplot(2,2,2)
795plot(100*evm_mat_B(:),'o','MarkerSize',1)
796axis tight
797hold on
798plot([1 length(evm_mat_B(:))], 100*[aevms_B, aevms_B],'r','LineWidth',4)
799hold off
800xlabel('Data Symbol Index')
801ylabel('EVM (%)');
802legend('Per-Symbol EVM','Average EVM','Location','NorthWest');
803title('Stream B')
804grid on
805myAxis_mat(2,:) = axis();
806
807subplot(2,2,1); axis([min(myAxis_mat(:,1)), max(myAxis_mat(:,2)) ...
808                      min(myAxis_mat(:,3)), min(myAxis_mat(:,4))]);
809myAxis = axis;
810h = text(round(.05*length(evm_mat_A(:))), 100*aevms_A+ .1*(myAxis(4)-myAxis(3)), sprintf('Effective SINR: %.1f dB', snr_A));
811set(h,'Color',[1 0 0])
812set(h,'FontWeight','bold')
813set(h,'FontSize',10)
814set(h,'EdgeColor',[1 0 0])
815set(h,'BackgroundColor',[1 1 1])
816subplot(2,2,2); axis([min(myAxis_mat(:,1)), max(myAxis_mat(:,2)) ...
817                      min(myAxis_mat(:,3)), min(myAxis_mat(:,4))]);
818myAxis = axis;
819h = text(round(.05*length(evm_mat_B(:))), 100*aevms_B+ .1*(myAxis(4)-myAxis(3)), sprintf('Effective SINR: %.1f dB', snr_B));
820set(h,'Color',[1 0 0])
821set(h,'FontWeight','bold')
822set(h,'FontSize',10)
823set(h,'EdgeColor',[1 0 0])
824set(h,'BackgroundColor',[1 1 1])
825
826
827subplot(2,2,3)
828imagesc(1:N_OFDM_SYMS, (SC_IND_DATA - N_SC/2), 100*fftshift(evm_mat_A,1))
829
830hold on
831h = line([1,N_OFDM_SYMS],[0,0]);
832set(h,'Color',[1 0 0])
833set(h,'LineStyle',':')
834hold off
835%myAxis = gca;
836%set(myAxis,'YTickLabel', fftshift(get(myAxis,'YTickLabel'),1))
837grid on
838xlabel('OFDM Symbol Index')
839ylabel('Subcarrier Index')
840title('Stream A')
841h = colorbar;
842set(get(h,'title'),'string','EVM (%)');
843myAxis = caxis();
844if (myAxis(2)-myAxis(1)) < 5
845    caxis([myAxis(1), myAxis(1)+5])
846end
847
848myAxis_mat = [];
849myAxis_mat(1,:) = caxis();
850
851subplot(2,2,4)
852imagesc(1:N_OFDM_SYMS, (SC_IND_DATA - N_SC/2), 100*fftshift(evm_mat_B,1))
853
854hold on
855h = line([1,N_OFDM_SYMS],[0,0]);
856set(h,'Color',[1 0 0])
857set(h,'LineStyle',':')
858hold off
859%myAxis = gca;
860%set(myAxis,'YTickLabel', fftshift(get(myAxis,'YTickLabel'),1))
861grid on
862xlabel('OFDM Symbol Index')
863ylabel('Subcarrier Index')
864title('Stream B')
865h = colorbar;
866set(get(h,'title'),'string','EVM (%)');
867myAxis = caxis();
868if (myAxis(2)-myAxis(1)) < 5
869    caxis([myAxis(1), myAxis(1)+5])
870end
871
872myAxis_mat(2,:) = caxis();
873
874subplot(2,2,3); caxis([min(myAxis_mat(:,1)), max(myAxis_mat(:,2))])
875
876if(WRITE_PNG_FILES)
877    print(gcf,sprintf('wl_mimo_ofdm_plots_%s_evm', example_mode_string), '-dpng', '-r96', '-painters')
878end
879
880%% Calculate Rx stats
881
882sym_errs = sum(tx_data ~= rx_data);
883bit_errs = length(find(dec2bin(bitxor(tx_data, rx_data),8) == '1'));
884rx_evm   = sqrt(sum((real(rx_syms) - real(tx_syms)).^2 + (imag(rx_syms) - imag(tx_syms)).^2)/(length(SC_IND_DATA) * N_OFDM_SYMS));
885
886fprintf('\nResults:\n');
887fprintf('Num Bytes:   %d\n', N_DATA_SYMS * log2(MOD_ORDER) / 8);
888fprintf('Sym Errors:  %d (of %d total symbols)\n', sym_errs, N_DATA_SYMS);
889fprintf('Bit Errors:  %d (of %d total bits)\n', bit_errs, N_DATA_SYMS * log2(MOD_ORDER));
890
891cfo_est_lts = rx_cfo_est_lts*(SAMP_FREQ/INTERP_RATE);
892cfo_est_phaseErr = mean(diff(unwrap(pilot_phase_err)))/(4e-6*2*pi);
893cfo_total_ppm = ((cfo_est_lts + cfo_est_phaseErr) /  ((2.412+(.005*(CHANNEL-1)))*1e9)) * 1e6;
894
895fprintf('CFO Est:     %3.2f kHz (%3.2f ppm)\n', (cfo_est_lts + cfo_est_phaseErr)*1e-3, cfo_total_ppm);
896fprintf('     LTS CFO Est:                  %3.2f kHz\n', cfo_est_lts*1e-3);
897fprintf('     Phase Error Residual CFO Est: %3.2f kHz\n', cfo_est_phaseErr*1e-3);
898
899if DO_APPLY_SFO_CORRECTION
900    drift_sec = pilot_slope_mat / (2*pi*312500);
901    sfo_est_ppm =  1e6*mean((diff(drift_sec) / 4e-6));
902    sfo_est = sfo_est_ppm*20;
903    fprintf('SFO Est:     %3.2f Hz (%3.2f ppm)\n', sfo_est, sfo_est_ppm);
904
905end
906
907
Note: See TracBrowser for help on using the repository browser.