1% tofdm.m 2% David Rowe and Steve Sampson June 2017 3% 4% Octave script for comparing Octave and C versions of OFDM modem 5% 6% If running from the Octave command line a good idea to clear globals before 7% each run: 8% 9% octave> clear; tofdm; 10 11% ------------------------------------------------------------------ 12 13Nframes = 10; 14sample_clock_offset_ppm = 100; 15foff_hz = 0.5; 16 17more off; format; 18ofdm_lib; 19autotest; 20ldpc 21global passes = 0; 22global fails = 0; 23 24% attempt to start up CML, path will be different on your machine 25 26path_to_cml = '~/cml'; 27addpath(strcat(path_to_cml, "/mex"), strcat(path_to_cml, "/mat")); 28cml_support = 0; 29if exist("Somap") == 0 30 printf("Can't find CML mex directory so we won't run those tests for now...\n"); 31else 32 printf("OK found CML mex directory so will add those tests...\n"); 33 cml_support = 1; 34end 35 36% --------------------------------------------------------------------- 37% Run Octave version 38% --------------------------------------------------------------------- 39 40% useful to test the modem at other Nc's, but if Nc != 17 we aren't set up for 41% LDPC testing so disable 42if getenv("NC") 43 Nc = str2num(getenv("NC")); 44 cml_support = 0; 45else 46 Nc = 17; 47end 48printf("Nc = %d LDPC testing: %d\n", Nc, cml_support); 49 50config = ofdm_init_mode("700D"); 51config.Nc = Nc; 52states = ofdm_init(config); 53states.verbose = 0; 54ofdm_load_const; 55 56printf("Nbitsperframe: %d\n", Nbitsperframe); 57 58if cml_support 59 Nuwtxtsymbolsperframe = (states.Nuwbits+states.Ntxtbits)/bps; 60 S_matrix = [1, j, -j, -1]; 61 EsNo = 10; 62 symbol_likelihood_log = bit_likelihood_log = detected_data_log = []; 63 64 % Set up LDPC code 65 66 mod_order = 4; bps = 2; modulation = 'QPSK'; mapping = 'gray'; 67 demod_type = 0; decoder_type = 0; max_iterations = 100; 68 69 load HRA_112_112.txt 70 [code_param framesize rate] = ldpc_init_user(HRA_112_112, modulation, mod_order, mapping); 71 assert(Nbitsperframe == (code_param.coded_bits_per_frame + states.Nuwbits + states.Ntxtbits)); 72end 73 74tx_bits = zeros(1,Nbitsperframe); 75rand('seed',1); 76 77payload_data_bits = round(rand(1,(Nbitsperframe-Nuwbits-Ntxtbits)/2)); 78states.mean_amp = 1; % start this with something sensible otherwise LDPC decode fails 79if cml_support 80 ibits = payload_data_bits; 81 codeword = LdpcEncode(ibits, code_param.H_rows, code_param.P_matrix); 82 tx_bits(Nuwbits+Ntxtbits+1:end) = codeword; 83 tx_bits(1:Nuwbits+Ntxtbits) = [states.tx_uw zeros(1,Ntxtbits)]; 84else 85 tx_bits = create_ldpc_test_frame(states, coded_frame=0); 86end 87 88% Run tx loop 89 90tx_bits_log = []; tx_log = []; 91for f=1:Nframes 92 tx_bits_log = [tx_bits_log tx_bits]; 93 tx_log = [tx_log ofdm_mod(states, tx_bits)]; 94end 95 96% Channel simulation ---------------------------------------------- 97 98rx_log = sample_clock_offset(tx_log, sample_clock_offset_ppm); 99rx_log = freq_shift(rx_log, foff_hz, Fs); 100 101% Rx --------------------------------------------------------------- 102 103% Init rx with ideal timing so we can test with timing estimation disabled 104 105Nsam = length(rx_log); 106prx = 1; 107nin = Nsamperframe+2*(M+Ncp); 108states.rxbuf(Nrxbuf-nin+1:Nrxbuf) = rx_log(prx:nin); 109prx += nin; 110 111rxbuf_log = []; rxbuf_in_log = []; rx_sym_log = []; foff_hz_log = []; 112timing_est_log = timing_valid_log = timing_mx_log = []; 113coarse_foff_est_hz_log = []; sample_point_log = []; 114phase_est_pilot_log = []; rx_amp_log = []; 115rx_np_log = []; rx_bits_log = []; 116snr_log = []; mean_amp_log = []; 117 118states.timing_en = 1; 119states.foff_est_en = 1; 120states.phase_est_en = 1; 121 122if states.timing_en == 0 123 % manually set ideal timing instant 124 states.sample_point = Ncp; 125end 126 127 128for f=1:Nframes 129 130 % insert samples at end of buffer, set to zero if no samples 131 % available to disable phase estimation on future pilots on last 132 % frame of simulation 133 134 nin = states.nin; 135 lnew = min(Nsam-prx+1,nin); 136 rxbuf_in = zeros(1,nin); 137 %printf("nin: %d prx: %d lnew: %d\n", nin, prx, lnew); 138 if lnew 139 rxbuf_in(1:lnew) = rx_log(prx:prx+lnew-1); 140 end 141 prx += lnew; 142 143 [states rx_bits achannel_est_pilot_log arx_np arx_amp] = ofdm_demod(states, rxbuf_in); 144 145 % log some states for comparison to C 146 147 rxbuf_in_log = [rxbuf_in_log rxbuf_in]; 148 rxbuf_log = [rxbuf_log states.rxbuf]; 149 rx_sym_log = [rx_sym_log; states.rx_sym]; 150 phase_est_pilot_log = [phase_est_pilot_log; angle(achannel_est_pilot_log)]; 151 rx_amp_log = [rx_amp_log arx_amp]; 152 foff_hz_log = [foff_hz_log; states.foff_est_hz]; 153 timing_est_log = [timing_est_log; states.timing_est]; 154 timing_valid_log = [timing_valid_log; states.timing_valid]; 155 timing_mx_log = [timing_mx_log; states.timing_mx]; 156 coarse_foff_est_hz_log = [coarse_foff_est_hz_log; states.coarse_foff_est_hz]; 157 sample_point_log = [sample_point_log; states.sample_point]; 158 rx_np_log = [rx_np_log arx_np]; 159 rx_bits_log = [rx_bits_log rx_bits]; 160 mean_amp_log = [mean_amp_log; states.mean_amp]; 161 EsNo_estdB = esno_est_calc(arx_np); 162 SNR_estdB = snr_from_esno(states, EsNo_estdB); 163 snr_log = [snr_log; SNR_estdB]; 164 165 % Optional testing of LDPC functions 166 167 if cml_support 168 mean_amp = states.mean_amp; 169 %mean_amp = 1; 170 symbol_likelihood = Demod2D(arx_np(Nuwtxtsymbolsperframe+1:end)/mean_amp, S_matrix, EsNo, arx_amp(Nuwtxtsymbolsperframe+1:end)/mean_amp); 171 bit_likelihood = Somap(symbol_likelihood); 172 173 [x_hat paritychecks] = MpDecode(-bit_likelihood(1:code_param.coded_bits_per_frame), code_param.H_rows, code_param.H_cols, max_iterations, decoder_type, 1, 1); 174 [mx mx_ind] = max(paritychecks); 175 detected_data = x_hat(mx_ind,:); 176 177 % make sure LDPC decoding is working OK 178 179 % assert(codeword == detected_data); 180 181 [m n] = size(symbol_likelihood); 182 symbol_likelihood_log = [symbol_likelihood_log; reshape(symbol_likelihood,m*n,1)]; 183 bit_likelihood_log = [bit_likelihood_log; bit_likelihood']; 184 detected_data_log = [detected_data_log detected_data]; 185 end 186 187end 188 189% --------------------------------------------------------------------- 190% Run C version and plot Octave and C states and differences 191% --------------------------------------------------------------------- 192 193printf("\nRunning C version....\n"); 194 195% Override default path by: 196% 1. if running from octave CLI: setting path_to_tofdm = "/your/path/to/tofdm" 197% 2. If running from shell....." set PATH_TO_OFDM = "/your/path/to/tofdm" 198 199if exist("path_to_tofdm", "var") == 0 200 path_to_tofdm = "../build_linux/unittest/tofdm" 201end 202 203if getenv("PATH_TO_TOFDM") 204 path_to_tofdm = getenv("PATH_TO_TOFDM") 205 printf("setting path from env var\n"); 206end 207 208path_to_tofdm = sprintf("%s --nc %d", path_to_tofdm, Nc); % append Nc for variable Nc tests 209 210if cml_support == 0 211 path_to_tofdm = sprintf("%s --noldpc", path_to_tofdm); 212end 213 214system(path_to_tofdm); 215load tofdm_out.txt; 216 217fg = 1; 218 219f = figure(fg++); clf; plot(rx_np_log,'+'); title('Octave Scatter Diagram'); axis([-1.5 1.5 -1.5 1.5]); 220f = figure(fg++); clf; plot(rx_np_log_c,'+'); title('C Scatter Diagram'); axis([-1.5 1.5 -1.5 1.5]); 221 222stem_sig_and_error(fg++, 111, tx_bits_log_c, tx_bits_log - tx_bits_log_c, 'tx bits', [1 length(tx_bits_log) -1.5 1.5]) 223 224stem_sig_and_error(fg, 211, real(tx_log_c), real(tx_log - tx_log_c), 'tx re', [1 length(tx_log_c) -0.1 0.1]) 225stem_sig_and_error(fg++, 212, imag(tx_log_c), imag(tx_log - tx_log_c), 'tx im', [1 length(tx_log_c) -0.1 0.1]) 226 227stem_sig_and_error(fg, 211, real(rx_log_c), real(rx_log - rx_log_c), 'rx re', [1 length(rx_log_c) -0.1 0.1]) 228stem_sig_and_error(fg++, 212, imag(rx_log_c), imag(rx_log - rx_log_c), 'rx im', [1 length(rx_log_c) -0.1 0.1]) 229 230stem_sig_and_error(fg, 211, real(rxbuf_in_log_c), real(rxbuf_in_log - rxbuf_in_log_c), 'rxbuf in re', [1 length(rxbuf_in_log_c) -0.1 0.1]) 231stem_sig_and_error(fg++, 212, imag(rxbuf_in_log_c), imag(rxbuf_in_log - rxbuf_in_log_c), 'rxbuf in im', [1 length(rxbuf_in_log_c) -0.1 0.1]) 232 233stem_sig_and_error(fg, 211, real(rxbuf_log_c), real(rxbuf_log - rxbuf_log_c), 'rxbuf re', [1 length(rxbuf_log_c) -0.1 0.1]) 234stem_sig_and_error(fg++, 212, imag(rxbuf_log_c), imag(rxbuf_log - rxbuf_log_c), 'rxbuf im', [1 length(rxbuf_log_c) -0.1 0.1]) 235 236stem_sig_and_error(fg, 211, real(rx_sym_log_c), real(rx_sym_log - rx_sym_log_c), 'rx sym re', [1 length(rx_sym_log_c) -1.5 1.5]) 237stem_sig_and_error(fg++, 212, imag(rx_sym_log_c), imag(rx_sym_log - rx_sym_log_c), 'rx sym im', [1 length(rx_sym_log_c) -1.5 1.5]) 238 239% for angles pi and -pi are the same 240 241d = phase_est_pilot_log - phase_est_pilot_log_c; d = angle(exp(j*d)); 242 243stem_sig_and_error(fg, 211, phase_est_pilot_log_c, d, 'phase est pilot', [1 length(phase_est_pilot_log_c) -1.5 1.5]) 244stem_sig_and_error(fg++, 212, rx_amp_log_c, rx_amp_log - rx_amp_log_c, 'rx amp', [1 length(rx_amp_log_c) -1.5 1.5]) 245 246stem_sig_and_error(fg , 211, foff_hz_log_c, (foff_hz_log - foff_hz_log_c), 'foff hz', [1 length(foff_hz_log_c) -1.5 1.5]) 247 248stem_sig_and_error(fg++, 212, timing_mx_log_c, (timing_mx_log - timing_mx_log_c), 'timing mx', [1 length(timing_mx_log_c) 0 2]) 249 250stem_sig_and_error(fg, 211, timing_est_log_c, (timing_est_log - timing_est_log_c), 'timing est', [1 length(timing_est_log_c) -1.5 1.5]) 251stem_sig_and_error(fg++, 212, sample_point_log_c, (sample_point_log - sample_point_log_c), 'sample point', [1 length(sample_point_log_c) -1.5 1.5]) 252 253stem_sig_and_error(fg++, 111, rx_bits_log_c, rx_bits_log - rx_bits_log_c, 'rx bits', [1 length(rx_bits_log) -1.5 1.5]) 254 255% Run through checklist ----------------------------- 256 257check(states.rate_fs_pilot_samples, pilot_samples_c, 'pilot_samples'); 258check(tx_bits_log, tx_bits_log_c, 'tx_bits'); 259check(tx_log, tx_log_c, 'tx'); 260check(rx_log, rx_log_c, 'rx'); 261check(rxbuf_in_log, rxbuf_in_log_c, 'rxbuf in'); 262check(rxbuf_log, rxbuf_log_c, 'rxbuf'); 263check(rx_sym_log, rx_sym_log_c, 'rx_sym', tol=10E-3); 264check(phase_est_pilot_log, phase_est_pilot_log_c, 'phase_est_pilot', tol=1E-2, its_an_angle=1); 265check(rx_amp_log, rx_amp_log_c, 'rx_amp'); 266check(timing_est_log, timing_est_log_c, 'timing_est'); 267check(timing_valid_log, timing_valid_log_c, 'timing_valid'); 268check(timing_mx_log, timing_mx_log_c, 'timing_mx'); 269check(coarse_foff_est_hz_log, coarse_foff_est_hz_log_c, 'coarse_foff_est_hz'); 270check(sample_point_log, sample_point_log_c, 'sample_point'); 271check(foff_hz_log, foff_hz_log_c, 'foff_est_hz'); 272check(rx_bits_log, rx_bits_log_c, 'rx_bits'); 273if cml_support 274 check(symbol_likelihood_log, symbol_likelihood_log_c, 'symbol_likelihood_log', tol=1E-2); 275 check(bit_likelihood_log, bit_likelihood_log_c, 'bit_likelihood_log'); 276 check(detected_data_log, detected_data_log_c, 'detected_data'); 277end 278check(mean_amp_log, mean_amp_log_c, 'mean_amp_log'); 279check(snr_log, snr_log_c, 'snr_log'); 280printf("\npasses: %d fails: %d\n", passes, fails); 281 282