1% cohpsk_dev.m 2% David Rowe Mar 2015 3% 4% Coherent PSK modem development and testing functions 5% 6 7cohpsk_lib; 8 9% Init HF channel model from stored sample files of spreading signal ---------------------------------- 10 11function [spread spread_2ms hf_gain] = init_hf_model(Fs, nsam) 12 13 % convert "spreading" samples from 1kHz carrier at Fss to complex 14 % baseband, generated by passing a 1kHz sine wave through PathSim 15 % with the ccir-poor model, enabling one path at a time. 16 17 Fc = 1000; Fss = 8000; 18 fspread = fopen("../raw/sine1k_2Hz_spread.raw","rb"); 19 spread1k = fread(fspread, "int16")/10000; 20 fclose(fspread); 21 fspread = fopen("../raw/sine1k_2ms_delay_2Hz_spread.raw","rb"); 22 spread1k_2ms = fread(fspread, "int16")/10000; 23 fclose(fspread); 24 25 % down convert to complex baseband 26 spreadbb = spread1k.*exp(-j*(2*pi*Fc/Fss)*(1:length(spread1k))'); 27 spreadbb_2ms = spread1k_2ms.*exp(-j*(2*pi*Fc/Fss)*(1:length(spread1k_2ms))'); 28 29 % remove -2000 Hz image 30 b = fir1(50, 5/Fss); 31 spread = filter(b,1,spreadbb); 32 spread_2ms = filter(b,1,spreadbb_2ms); 33 34 % discard first 1000 samples as these were near 0, probably as 35 % PathSim states were ramping up 36 37 spread = spread(1000:length(spread)); 38 spread_2ms = spread_2ms(1000:length(spread_2ms)); 39 40 % change output samples so they are at rate Fs reqd by caller 41 42 spread = resample(spread, Fs, Fss); 43 spread_2ms = resample(spread_2ms, Fs, Fss); 44 45 % Determine "gain" of HF channel model, so we can normalise 46 % carrier power during HF channel sim to calibrate SNR. I imagine 47 % different implementations of ccir-poor would do this in 48 % different ways, leading to different BER results. Oh Well! 49 50 hf_gain = 1.0/sqrt(var(spread(1:nsam))+var(spread_2ms(1:nsam))); 51endfunction 52 53 54function write_pilot_file(pilot, Nsymbrowpilot, Ns, Nsymrow, Npilotsframe, Nc); 55 56 filename = sprintf("../src/cohpsk_defs.h", Npilotsframe, Nc); 57 f=fopen(filename,"wt"); 58 fprintf(f,"/* Generated by write_pilot_file() Octave function */\n\n"); 59 fprintf(f,"#define NSYMROW %d /* number of data symbols on each row (i.e. each carrier) */\n", Nsymrow); 60 fprintf(f,"#define NS %d /* number of data symbols between pilots */\n", Ns); 61 fprintf(f,"#define NPILOTSFRAME %d /* number of pilot symbols on each row */\n", Npilotsframe); 62 fprintf(f,"#define PILOTS_NC %d /* number of carriers */\n\n", Nc); 63 fprintf(f,"#define NSYMROWPILOT %d /* length of row after pilots inserted */\n\n", Nsymbrowpilot); 64 fclose(f); 65 66 filename = sprintf("../src/pilots_coh.h", Npilotsframe, Nc); 67 f=fopen(filename,"wt"); 68 fprintf(f,"/* Generated by write_pilot_file() Octave function */\n\n"); 69 fprintf(f,"float pilots_coh[][PILOTS_NC]={\n"); 70 for r=1:Npilotsframe 71 fprintf(f, " {"); 72 for c=1:Nc-1 73 fprintf(f, " %f,", pilot(r, c)); 74 end 75 if r < Npilotsframe 76 fprintf(f, " %f},\n", pilot(r, Nc)); 77 else 78 fprintf(f, " %f}\n};", pilot(r, Nc)); 79 end 80 end 81 fclose(f); 82endfunction 83 84 85% Save test bits frame to a text file in the form of a C array 86 87function test_bits_coh_file(test_bits_coh) 88 89 f=fopen("../src/test_bits_coh.h","wt"); 90 fprintf(f,"/* Generated by test_bits_coh_file() Octave function */\n\n"); 91 fprintf(f,"const int test_bits_coh[]={\n"); 92 for m=1:length(test_bits_coh)-1 93 fprintf(f," %d,\n",test_bits_coh(m)); 94 endfor 95 fprintf(f," %d\n};\n",test_bits_coh(length(test_bits_coh))); 96 fclose(f); 97 98endfunction 99 100 101% Rate Rs BER tests ------------------------------------------------------------------------------ 102 103function sim_out = ber_test(sim_in) 104 sim_in = symbol_rate_init(sim_in); 105 106 Fs = sim_in.Fs; 107 Rs = sim_in.Rs; 108 Ntrials = sim_in.Ntrials; 109 verbose = sim_in.verbose; 110 plot_scatter = sim_in.plot_scatter; 111 framesize = sim_in.framesize; 112 bps = sim_in.bps; 113 114 Esvec = sim_in.Esvec; 115 ldpc_code = sim_in.ldpc_code; 116 rate = sim_in.ldpc_code_rate; 117 code_param = sim_in.code_param; 118 tx_bits_buf = sim_in.tx_bits_buf; 119 Nsymb = sim_in.Nsymb; 120 Nsymbrow = sim_in.Nsymbrow; 121 Nsymbrowpilot = sim_in.Nsymbrowpilot; 122 Nc = sim_in.Nc; 123 Npilotsframe = sim_in.Npilotsframe; 124 Ns = sim_in.Ns; 125 Np = sim_in.Np; 126 Nd = sim_in.Nd; 127 modulation = sim_in.modulation; 128 pilot = sim_in.pilot; 129 prev_sym_tx = sim_in.prev_sym_tx; 130 prev_sym_rx = sim_in.prev_sym_rx; 131 rx_symb_buf = sim_in.rx_symb_buf; 132 tx_pilot_buf = sim_in.tx_pilot_buf; 133 rx_pilot_buf = sim_in.rx_pilot_buf; 134 135 hf_sim = sim_in.hf_sim; 136 nhfdelay = sim_in.hf_delay_ms*Rs/1000; 137 hf_mag_only = sim_in.hf_mag_only; 138 f_off = sim_in.f_off; 139 div_time_shift = sim_in.div_timeshift; 140 141 [spread spread_2ms hf_gain] = init_hf_model(Rs, Nsymbrowpilot*(Ntrials+2)); 142 143 if strcmp(modulation,'dqpsk') 144 Nsymbrowpilot = Nsymbrow; 145 end 146 147 % Start Simulation ---------------------------------------------------------------- 148 149 for ne = 1:length(Esvec) 150 EsNodB = Esvec(ne); 151 EsNo = 10^(EsNodB/10); 152 153 variance = 1/EsNo; 154 if verbose > 1 155 printf("EsNo (dB): %f EsNo: %f variance: %f\n", EsNodB, EsNo, variance); 156 end 157 158 Terrs = 0; Tbits = 0; 159 160 s_ch_tx_log = []; 161 rx_symb_log = []; 162 noise_log = []; 163 errors_log = []; 164 Nerrs_log = []; 165 phi_log = []; 166 amp_log = []; 167 EsNo__log = []; 168 169 ldpc_errors_log = []; ldpc_Nerrs_log = []; 170 171 Terrsldpc = Tbitsldpc = Ferrsldpc = 0; 172 173 % init HF channel 174 175 hf_n = 1; 176 177 phase_offset_rect = 1; 178 w_offset = 2*pi*f_off/Rs; 179 w_offset_rect = exp(j*w_offset); 180 181 ct_symb_buf = zeros(2*Nsymbrowpilot, Nc*Nd); 182 prev_tx_symb = prev_rx_symb = ones(1, Nc*Nd); 183 184 % simulation starts here----------------------------------- 185 186 for nn = 1:Ntrials+2 187 188 if ldpc_code 189 tx_bits = round(rand(1,framesize*rate)); 190 else 191 tx_bits = round(rand(1,framesize)); 192 end 193 194 if strcmp(modulation,'qpsk') 195 [tx_symb tx_bits] = bits_to_qpsk_symbols(sim_in, tx_bits, code_param); 196 197 % one frame delay on bits for qpsk 198 199 tx_bits_buf(1:framesize) = tx_bits_buf(framesize+1:2*framesize); 200 tx_bits_buf(framesize+1:2*framesize) = tx_bits; 201 202 end 203 if strcmp(modulation, 'dqpsk') 204 [tx_symb prev_tx_symb] = bits_to_dqpsk_symbols(sim_in, tx_bits, prev_tx_symb); 205 tx_bits_buf(1:framesize) = tx_bits; 206 end 207 208 s_ch = tx_symb; 209 210 % HF channel simulation ------------------------------------ 211 212 hf_fading = ones(1,Nsymb); 213 if hf_sim 214 215 % separation between carriers. Note this effectively 216 % under samples at Rs, I dont think this matters. 217 % Equivalent to doing freq shift at Fs, then 218 % decimating to Rs. 219 220 wsep = 2*pi*(1+0.5); % e.g. 75Hz spacing at Rs=50Hz, alpha=0.5 filters 221 222 hf_model(hf_n, :) = zeros(1,Nc*Nd); 223 224 for r=1:Nsymbrowpilot 225 for c=1:Nd*Nc 226 if c > Nc 227 time_shift = sim_in.div_timeshift; 228 else 229 time_shift = 1; 230 end 231 ahf_model = hf_gain*(spread(hf_n+time_shift) + exp(-j*c*wsep*nhfdelay)*spread_2ms(hf_n+time_shift)); 232 233 if hf_mag_only 234 s_ch(r,c) *= abs(ahf_model); 235 else 236 s_ch(r,c) *= ahf_model; 237 end 238 hf_model(hf_n, c) = ahf_model; 239 end 240 hf_n++; 241 end 242 end 243 244 % keep a record of each tx symbol so we can check average power 245 246 for r=1:Nsymbrow 247 for c=1:Nd*Nc 248 s_ch_tx_log = [s_ch_tx_log s_ch(r,c)]; 249 end 250 end 251 252 % AWGN noise and phase/freq offset channel simulation 253 % 0.5 factor ensures var(noise) == variance , i.e. splits power between Re & Im 254 255 noise = sqrt(variance*0.5)*(randn(Nsymbrowpilot,Nc*Nd) + j*randn(Nsymbrowpilot,Nc*Nd)); 256 noise_log = [noise_log noise]; 257 258 for r=1:Nsymbrowpilot 259 s_ch(r,:) *= phase_offset_rect; 260 phase_offset_rect *= w_offset_rect; 261 end 262 s_ch += noise; 263 264 ct_symb_buf(1:Nsymbrowpilot,:) = ct_symb_buf(Nsymbrowpilot+1:2*Nsymbrowpilot,:); 265 ct_symb_buf(Nsymbrowpilot+1:2*Nsymbrowpilot,:) = s_ch; 266 267 if strcmp(modulation,'qpsk') 268 [rx_symb rx_bits rx_symb_linear amp_ phi_ sig_rms noise_rms sim_in] = qpsk_symbols_to_bits(sim_in, ct_symb_buf(1:Nsymbrowpilot+Npilotsframe,:)); 269 phi_log = [phi_log; phi_]; 270 amp_log = [amp_log; amp_]; 271 end 272 if strcmp(modulation,'dqpsk') 273 [rx_symb rx_bits rx_symb_linear prev_rx_symb] = dqpsk_symbols_to_bits(sim_in, s_ch, prev_rx_symb); 274 end 275 276 % Wait until we have enough frames to do pilot assisted phase estimation 277 278 if nn > 1 279 rx_symb_log = [rx_symb_log rx_symb_linear]; 280 %EsNo__log = [EsNo__log EsNo_]; 281 282 % Measure BER 283 284 error_positions = xor(rx_bits, tx_bits_buf(1:framesize)); 285 Nerrs = sum(error_positions); 286 Terrs += Nerrs; 287 Tbits += length(tx_bits); 288 errors_log = [errors_log error_positions]; 289 Nerrs_log = [Nerrs_log Nerrs]; 290 291 % Optionally LDPC decode 292 293 if ldpc_code 294 detected_data = ldpc_dec(code_param, sim_in.max_iterations, sim_in.demod_type, sim_in.decoder_type, ... 295 rx_symb_linear, min(100,EsNo_), amp_linear); 296 error_positions = xor( detected_data(1:framesize*rate), tx_bits_buf(1:framesize*rate) ); 297 Nerrs = sum(error_positions); 298 ldpc_Nerrs_log = [ldpc_Nerrs_log Nerrs]; 299 ldpc_errors_log = [ldpc_errors_log error_positions]; 300 if Nerrs 301 Ferrsldpc++; 302 end 303 Terrsldpc += Nerrs; 304 Tbitsldpc += framesize*rate; 305 end 306 end 307 end 308 309 TERvec(ne) = Terrs; 310 BERvec(ne) = Terrs/Tbits; 311 312 if verbose 313 av_tx_pwr = (s_ch_tx_log * s_ch_tx_log')/length(s_ch_tx_log); 314 315 printf("EsNo (dB): %3.1f Terrs: %d Tbits: %d BER %5.3f QPSK BER theory %5.3f av_tx_pwr: %3.2f", 316 EsNodB, Terrs, Tbits, 317 Terrs/Tbits, 0.5*erfc(sqrt(EsNo/2)), av_tx_pwr); 318 if ldpc_code 319 printf("\n LDPC: Terrs: %d BER: %4.2f Ferrs: %d FER: %4.2f", 320 Terrsldpc, Terrsldpc/Tbitsldpc, Ferrsldpc, Ferrsldpc/Ntrials); 321 end 322 printf("\n"); 323 end 324 end 325 326 Ebvec = Esvec - 10*log10(bps); 327 sim_out.BERvec = BERvec; 328 sim_out.Ebvec = Ebvec; 329 sim_out.TERvec = TERvec; 330 sim_out.errors_log = errors_log; 331 sim_out.ldpc_errors_log = ldpc_errors_log; 332 333 if plot_scatter 334 figure(2); 335 clf; 336 scat = rx_symb_log .* exp(j*pi/4); 337 plot(real(scat), imag(scat),'+'); 338 title('Scatter plot'); 339 a = 1.5*max(real(scat)); b = 1.5*max(imag(scat)); 340 axis([-a a -b b]); 341 342 if hf_sim 343 figure(3); 344 clf; 345 346 y = 1:(hf_n-1); 347 x = 1:Nc*Nd; 348 EsNodBSurface = 20*log10(abs(hf_model(y,:))) - 10*log10(variance); 349 EsNodBSurface(find(EsNodBSurface < -5)) = -5; 350 EsNodBSurface(find(EsNodBSurface > 25)) = 25; 351 mesh(x,y,EsNodBSurface); 352 grid 353 axis([1 Nc*Nd 1 Rs*5 -5 25]) 354 title('HF Channel Es/No'); 355 356 if verbose 357 [m n] = size(hf_model); 358 av_hf_pwr = sum(sum(abs(hf_model(:,:)).^2))/(m*n); 359 printf("average HF power: %3.2f over %d symbols\n", av_hf_pwr, m*n); 360 end 361 362 end 363 364 if strcmp(modulation,'qpsk') 365 % set up time axis to include gaps for pilots 366 367 [m1 n1] = size(phi_log); 368 phi_x = []; 369 phi_x_counter = 1; 370 p = Ns; 371 for r=1:m1 372 if p == Ns 373 phi_x_counter += Npilotsframe; 374 p = 0; 375 end 376 p++; 377 phi_x = [phi_x phi_x_counter++]; 378 end 379 380 phi_x -= Nsymbrowpilot; % account for delay in pilot buffer 381 382 figure(5); 383 clf 384 subplot(211) 385 [m n] = size(phi_log); 386 plot(phi_x, phi_log(:,2),'r+;Estimated HF channel phase;') 387 if hf_sim 388 hold on; 389 [m n] = size(hf_model); 390 plot(angle(hf_model(1:m,2)),'g;HF channel phase;') 391 hold off; 392 end 393 ylabel('Phase (rads)'); 394 legend('boxoff'); 395 axis([1 m -1.1*pi 1.1*pi]) 396 397 subplot(212) 398 plot(phi_x, amp_log(:,2),'r+;Estimated HF channel amp;') 399 if hf_sim 400 hold on; 401 plot(abs(hf_model(1:m,2))) 402 hold off; 403 end 404 ylabel('Amplitude'); 405 xlabel('Time (symbols)'); 406 legend('boxoff'); 407 axis([1 m 0 3]) 408 end 409 410 figure(4) 411 clf 412 stem(Nerrs_log) 413 axis([1 length(Nerrs_log) 0 max(Nerrs_log)+1]) 414 end 415 416endfunction 417 418function sim_in = standard_init 419 sim_in.verbose = 1; 420 sim_in.do_write_pilot_file = 0; 421 sim_in.plot_scatter = 0; 422 423 sim_in.Esvec = 50; 424 sim_in.Ntrials = 30; 425 sim_in.framesize = 2; 426 sim_in.Rs = 50; 427 428 sim_in.phase_offset = 0; 429 sim_in.w_offset = 0; 430 sim_in.phase_noise_amp = 0; 431 432 sim_in.hf_delay_ms = 2; 433 sim_in.hf_sim = 0; 434 sim_in.hf_mag_only = 0; 435 436 sim_in.Nd = 1; 437endfunction 438 439 440