1% cohpsk_lib.m 2% David Rowe Mar 2015 3% 4% Coherent PSK modem functions 5% 6 71; 8 9% Gray coded QPSK modulation function 10 11function symbol = qpsk_mod(two_bits) 12 two_bits_decimal = sum(two_bits .* [2 1]); 13 switch(two_bits_decimal) 14 case (0) symbol = 1; 15 case (1) symbol = j; 16 case (2) symbol = -j; 17 case (3) symbol = -1; 18 endswitch 19endfunction 20 21 22% Gray coded QPSK demodulation function 23 24function two_bits = qpsk_demod(symbol) 25 if isscalar(symbol) == 0 26 printf("only works with scalars\n"); 27 return; 28 end 29 bit0 = real(symbol*exp(j*pi/4)) < 0; 30 bit1 = imag(symbol*exp(j*pi/4)) < 0; 31 two_bits = [bit1 bit0]; 32endfunction 33 34% init function for symbol rate processing -------------------------------------------------------- 35 36function sim_in = symbol_rate_init(sim_in) 37 sim_in.Fs = Fs = 8000; 38 39 modulation = sim_in.modulation; 40 verbose = sim_in.verbose; 41 framesize = sim_in.framesize; 42 Ntrials = sim_in.Ntrials; 43 Esvec = sim_in.Esvec; 44 phase_offset = sim_in.phase_offset; 45 w_offset = sim_in.w_offset; 46 plot_scatter = sim_in.plot_scatter; 47 48 Rs = sim_in.Rs; 49 Nc = sim_in.Nc; 50 51 hf_sim = sim_in.hf_sim; 52 nhfdelay = sim_in.hf_delay_ms*Rs/1000; 53 hf_mag_only = sim_in.hf_mag_only; 54 55 Nd = sim_in.Nd; % diveristy 56 Ns = sim_in.Ns; % step size between pilots 57 ldpc_code = sim_in.ldpc_code; 58 rate = sim_in.ldpc_code_rate; 59 60 sim_in.bps = bps = 2; 61 62 sim_in.Nsymb = Nsymb = framesize/bps; 63 sim_in.Nsymbrow = Nsymbrow = Nsymb/Nc; 64 sim_in.Npilotsframe = Npilotsframe = 2; 65 sim_in.Nsymbrowpilot = Nsymbrowpilot = Nsymbrow + Npilotsframe; 66 67 if verbose == 2 68 printf("Each frame contains %d data bits or %d data symbols, transmitted as %d symbols by %d carriers.", framesize, Nsymb, Nsymbrow, Nc); 69 printf(" There are %d pilot symbols in each carrier together at the start of each frame, then %d data symbols.", Npilotsframe, Ns); 70 printf(" Including pilots, the frame is %d symbols long by %d carriers.\n\n", Nsymbrowpilot, Nc); 71 end 72 73 sim_in.prev_sym_tx = qpsk_mod([0 0])*ones(1,Nc*Nd); 74 sim_in.prev_sym_rx = qpsk_mod([0 0])*ones(1,Nc*Nd); 75 76 sim_in.rx_symb_buf = zeros(3*Nsymbrow, Nc*Nd); 77 sim_in.rx_pilot_buf = zeros(3*Npilotsframe,Nc*Nd); 78 sim_in.tx_bits_buf = zeros(1,2*framesize); 79 80 % pilot sequence is used for phase and amplitude estimation, and frame sync 81 82 pilot = 1 - 2*(rand(Npilotsframe,Nc) > 0.5); 83 sim_in.pilot = pilot; 84 sim_in.tx_pilot_buf = [pilot; pilot; pilot]; 85 86 if sim_in.do_write_pilot_file 87 write_pilot_file(pilot, Nsymbrowpilot, Ns, Nsymbrow, Npilotsframe, Nc); 88 end 89 90 % we use first 2 pilots of next frame to help with frame sync and fine freq 91 92 sim_in.Nct_sym_buf = 2*Nsymbrowpilot + 2; 93 sim_in.ct_symb_buf = zeros(sim_in.Nct_sym_buf, Nc*Nd); 94 95 sim_in.ff_phase = 1; 96 97 sim_in.ct_symb_ff_buf = zeros(Nsymbrowpilot + 2, Nc*Nd); 98 99 % Init LDPC -------------------------------------------------------------------- 100 101 if ldpc_code 102 % Start CML library 103 104 currentdir = pwd; 105 addpath '~/cml/mat' % assume the source files stored here 106 cd ~/cml 107 CmlStartup % note that this is not in the cml path! 108 cd(currentdir) 109 110 % Our LDPC library 111 112 ldpc; 113 114 mod_order = 4; 115 modulation2 = 'QPSK'; 116 mapping = 'gray'; 117 118 sim_in.demod_type = 0; 119 sim_in.decoder_type = 0; 120 sim_in.max_iterations = 100; 121 122 code_param = ldpc_init(rate, framesize, modulation2, mod_order, mapping); 123 code_param.code_bits_per_frame = framesize; 124 code_param.symbols_per_frame = framesize/bps; 125 sim_in.code_param = code_param; 126 else 127 sim_in.rate = 1; 128 sim_in.code_param = []; 129 end 130endfunction 131 132 133% Symbol rate processing for tx side (modulator) ------------------------------------------------------- 134 135% legacy DQPSK mod for comparative testing 136 137function [tx_symb prev_tx_symb] = bits_to_dqpsk_symbols(sim_in, tx_bits, prev_tx_symb) 138 Nc = sim_in.Nc; 139 Nsymbrow = sim_in.Nsymbrow; 140 141 tx_symb = zeros(Nsymbrow,Nc); 142 143 for c=1:Nc 144 for r=1:Nsymbrow 145 i = (c-1)*Nsymbrow + r; 146 tx_symb(r,c) = qpsk_mod(tx_bits(2*(i-1)+1:2*i)); 147 tx_symb(r,c) *= prev_tx_symb(c); 148 prev_tx_symb(c) = tx_symb(r,c); 149 end 150 end 151 152endfunction 153 154 155% legacy DQPSK demod for comparative testing 156 157function [rx_symb rx_bits rx_symb_linear prev_rx_symb] = dqpsk_symbols_to_bits(sim_in, rx_symb, prev_rx_symb) 158 Nc = sim_in.Nc; 159 Nsymbrow = sim_in.Nsymbrow; 160 161 tx_symb = zeros(Nsymbrow,Nc); 162 163 for c=1:Nc 164 for r=1:Nsymbrow 165 tmp = rx_symb(r,c); 166 rx_symb(r,c) *= conj(prev_rx_symb(c))/abs(prev_rx_symb(c)); 167 prev_rx_symb(c) = tmp; 168 i = (c-1)*Nsymbrow + r; 169 rx_symb_linear(i) = rx_symb(r,c); 170 rx_bits((2*(i-1)+1):(2*i)) = qpsk_demod(rx_symb(r,c)); 171 end 172 end 173 174endfunction 175 176 177function [tx_symb tx_bits] = bits_to_qpsk_symbols(sim_in, tx_bits, code_param) 178 ldpc_code = sim_in.ldpc_code; 179 rate = sim_in.ldpc_code_rate; 180 framesize = sim_in.framesize; 181 Nsymbrow = sim_in.Nsymbrow; 182 Nsymbrowpilot = sim_in.Nsymbrowpilot; 183 Nc = sim_in.Nc; 184 Npilotsframe = sim_in.Npilotsframe; 185 Ns = sim_in.Ns; 186 modulation = sim_in.modulation; 187 pilot = sim_in.pilot; 188 Nd = sim_in.Nd; 189 190 if ldpc_code 191 [tx_bits, tmp] = ldpc_enc(tx_bits, code_param); 192 end 193 194 % modulate -------------------------------------------- 195 196 % organise symbols into a Nsymbrow rows by Nc cols 197 % data and parity bits are on separate carriers 198 199 tx_symb = zeros(Nsymbrow,Nc); 200 201 for c=1:Nc 202 for r=1:Nsymbrow 203 i = (c-1)*Nsymbrow + r; 204 tx_symb(r,c) = qpsk_mod(tx_bits(2*(i-1)+1:2*i)); 205 end 206 end 207 208 % insert pilots at start of frame 209 210 tx_symb = [pilot(1,:); pilot(2,:); tx_symb;]; 211 212 % copy to other carriers (diversity) 213 214 tmp = tx_symb; 215 for d=1:Nd-1 216 tmp = [tmp tx_symb]; 217 end 218 tx_symb = tmp; 219 220 % ensures energy/symbol is normalised with diversity 221 222 tx_symb = tx_symb/sqrt(Nd); 223endfunction 224 225 226% Symbol rate processing for rx side (demodulator) ------------------------------------------------------- 227 228function [rx_symb rx_bits rx_symb_linear amp_ phi_ sig_rms noise_rms cohpsk] = qpsk_symbols_to_bits(cohpsk, ct_symb_buf) 229 framesize = cohpsk.framesize; 230 Nsymb = cohpsk.Nsymb; 231 Nsymbrow = cohpsk.Nsymbrow; 232 Nsymbrowpilot = cohpsk.Nsymbrowpilot; 233 Nc = cohpsk.Nc; 234 Nd = cohpsk.Nd; 235 Npilotsframe = cohpsk.Npilotsframe; 236 pilot = cohpsk.pilot; 237 verbose = cohpsk.verbose; 238 coh_en = cohpsk.coh_en; 239 240 % Use pilots to get phase and amplitude estimates We assume there 241 % are two samples at the start of each frame and two at the end 242 % Note: correlation (averging) method was used initially, but was 243 % poor at tracking fast phase changes that we experience on fading 244 % channels. Linear regression (fitting a straight line) works 245 % better on fading channels, but increases BER slightly for AWGN 246 % channels. 247 248 sampling_points = [1 2 cohpsk.Nsymbrow+3 cohpsk.Nsymbrow+4]; 249 pilot2 = [cohpsk.pilot(1,:); cohpsk.pilot(2,:); cohpsk.pilot(1,:); cohpsk.pilot(2,:);]; 250 phi_ = zeros(Nsymbrow, Nc*Nd); 251 amp_ = zeros(Nsymbrow, Nc*Nd); 252 253 for c=1:Nc*Nd 254 y = ct_symb_buf(sampling_points,c) .* pilot2(:,c-Nc*floor((c-1)/Nc)); 255 [m b] = linreg(sampling_points, y, length(sampling_points)); 256 yfit = m*[3 4 5 6] + b; 257 phi_(:, c) = angle(yfit); 258 259 mag = sum(abs(ct_symb_buf(sampling_points,c))); 260 amp_(:, c) = mag/length(sampling_points); 261 end 262 263 % now correct phase of data symbols 264 265 rx_symb = zeros(Nsymbrow, Nc); 266 rx_symb_linear = zeros(1, Nsymbrow*Nc*Nd); 267 rx_bits = zeros(1, framesize); 268 for c=1:Nc*Nd 269 for r=1:Nsymbrow 270 if coh_en == 1 271 rx_symb(r,c) = ct_symb_buf(2+r,c)*exp(-j*phi_(r,c)); 272 else 273 rx_symb(r,c) = ct_symb_buf(2+r,c); 274 end 275 i = (c-1)*Nsymbrow + r; 276 rx_symb_linear(i) = rx_symb(r,c); 277 end 278 end 279 280 % and finally optional diversity combination and make decn on bits 281 282 for c=1:Nc 283 for r=1:Nsymbrow 284 i = (c-1)*Nsymbrow + r; 285 div_symb = rx_symb(r,c); 286 for d=1:Nd-1 287 div_symb += rx_symb(r,c + Nc*d); 288 end 289 rx_bits((2*(i-1)+1):(2*i)) = qpsk_demod(div_symb); 290 end 291 end 292 293 % Estimate noise power from demodulated symbols. One method is to 294 % calculate the distance of each symbol from the average symbol 295 % position. However this is complicated by fading, which means the 296 % amplitude of the symbols is constantly changing. 297 298 % Now the scatter diagram in a fading channel is a X or cross 299 % shape. The noise can be resolved into two components at right 300 % angles to each other. The component along the the "thickness" 301 % of the arms is proportional to the noise power and not affected 302 % by fading. We only use points further along the real axis than 303 % the mean amplitude so we keep out of the central nosiey blob 304 305 sig_rms = mean(abs(rx_symb_linear)); 306 307 sum_x = 0; 308 sum_xx = 0; 309 n = 0; 310 for i=1:Nsymb*Nd 311 s = rx_symb_linear(i); 312 if abs(real(s)) > sig_rms 313 sum_x += imag(s); 314 sum_xx += imag(s)*imag(s); 315 n++; 316 end 317 end 318 319 noise_var = 0; 320 if n > 1 321 noise_var = (n*sum_xx - sum_x*sum_x)/(n*(n-1)); 322 end 323 noise_rms = sqrt(noise_var); 324endfunction 325 326function [ch_symb rx_timing rx_filt rx_baseband afdmdv f_est] = rate_Fs_rx_processing(afdmdv, ch_fdm_frame, f_est, nsymb, nin, freq_track) 327 M = afdmdv.M; 328 329 rx_baseband = []; 330 rx_filt = []; 331 rx_timing = []; 332 333 ch_fdm_frame_index = 1; 334 335 for r=1:nsymb 336 % shift signal to nominal baseband, this will put Nc/2 carriers either side of 0 Hz 337 338 [rx_fdm_frame_bb afdmdv.fbb_phase_rx] = freq_shift(ch_fdm_frame(ch_fdm_frame_index:ch_fdm_frame_index + nin - 1), -f_est, afdmdv.Fs, afdmdv.fbb_phase_rx); 339 ch_fdm_frame_index += nin; 340 341 % downconvert each FDM carrier to Nc separate baseband signals 342 343 [arx_baseband afdmdv] = fdm_downconvert(afdmdv, rx_fdm_frame_bb, nin); 344 [arx_filt afdmdv] = rx_filter(afdmdv, arx_baseband, nin); 345 [rx_onesym arx_timing env afdmdv] = rx_est_timing(afdmdv, arx_filt, nin); 346 347 rx_baseband = [rx_baseband arx_baseband]; 348 rx_filt = [rx_filt arx_filt]; 349 rx_timing = [rx_timing arx_timing]; 350 351 ch_symb(r,:) = rx_onesym; 352 353 % we only allow a timing shift on one symbol per frame 354 355 if nin != M 356 nin = M; 357 end 358 359 % freq tracking, see test_ftrack.m for unit test. Placed in 360 % this function as it needs to work on a symbol by symbol basis 361 % rather than frame by frame. This means the control loop 362 % operates at a sample rate of Rs = 50Hz for say 1 Hz/s drift. 363 364 if freq_track 365 beta = 0.005; 366 g = 0.2; 367 368 % combine difference on phase from last symbol over Nc carriers 369 370 mod_strip = 0; 371 for c=1:afdmdv.Nc+1 372 adiff = rx_onesym(c) .* conj(afdmdv.prev_rx_symb(c)); 373 afdmdv.prev_rx_symb(c) = rx_onesym(c); 374 375 % 4th power strips QPSK modulation, by multiplying phase by 4 376 % Using the abs value of the real coord was found to help 377 % non-linear issues when noise power was large 378 379 amod_strip = adiff.^4; 380 amod_strip = abs(real(amod_strip)) + j*imag(amod_strip); 381 mod_strip += amod_strip; 382 end 383 384 % loop filter made up of 1st order IIR plus integrator. Integerator 385 % was found to be reqd 386 387 afdmdv.filt = (1-beta)*afdmdv.filt + beta*angle(mod_strip); 388 f_est += g*afdmdv.filt; 389 end 390 end 391endfunction 392 393 394function ct_symb_buf = update_ct_symb_buf(ct_symb_buf, ch_symb, Nct_sym_buf, Nsymbrowpilot) 395 396 % update memory in symbol buffer 397 398 for r=1:Nct_sym_buf-Nsymbrowpilot 399 ct_symb_buf(r,:) = ct_symb_buf(r+Nsymbrowpilot,:); 400 end 401 i = 1; 402 for r=Nct_sym_buf-Nsymbrowpilot+1:Nct_sym_buf 403 ct_symb_buf(r,:) = ch_symb(i,:); 404 i++; 405 end 406endfunction 407 408 409% returns index of start of frame and fine freq offset 410 411function [next_sync cohpsk] = frame_sync_fine_freq_est(cohpsk, ch_symb, sync, next_sync) 412 ct_symb_buf = cohpsk.ct_symb_buf; 413 Nct_sym_buf = cohpsk.Nct_sym_buf; 414 Rs = cohpsk.Rs; 415 Nsymbrowpilot = cohpsk.Nsymbrowpilot; 416 Nc = cohpsk.Nc; 417 Nd = cohpsk.Nd; 418 419 ct_symb_buf = update_ct_symb_buf(ct_symb_buf, ch_symb, Nct_sym_buf, Nsymbrowpilot); 420 cohpsk.ct_symb_buf = ct_symb_buf; 421 422 % sample pilots at start of this frame and start of next frame 423 424 sampling_points = [1 2 cohpsk.Nsymbrow+3 cohpsk.Nsymbrow+4]; 425 pilot2 = [ cohpsk.pilot(1,:); cohpsk.pilot(2,:); cohpsk.pilot(1,:); cohpsk.pilot(2,:);]; 426 427 if sync == 0 428 429 % sample correlation over 2D grid of time and fine freq points 430 431 max_corr = 0; 432 for f_fine=-20:0.25:20 433 f_fine_rect = exp(-j*f_fine*2*pi*sampling_points/Rs)'; % note: this could be pre-computed at init or compile time 434 for t=0:cohpsk.Nsymbrowpilot-1 435 corr = 0; mag = 0; 436 for c=1:Nc*Nd 437 f_corr_vec = f_fine_rect .* ct_symb_buf(t+sampling_points,c); % note: this could be pre-computed at init or compile time 438 acorr = 0.0; 439 for p=1:length(sampling_points) 440 acorr += pilot2(p,c-Nc*floor((c-1)/Nc)) * f_corr_vec(p); 441 mag += abs(f_corr_vec(p)); 442 end 443 corr += abs(acorr); 444 end 445 446 if corr >= max_corr 447 max_corr = corr; 448 max_mag = mag; 449 cohpsk.ct = t; 450 cohpsk.f_fine_est = f_fine; 451 cohpsk.ff_rect = exp(-j*f_fine*2*pi/Rs); 452 end 453 end 454 end 455 456 printf(" [%d] fine freq f: %f max_ratio: %f ct: %d\n", cohpsk.frame, cohpsk.f_fine_est, abs(max_corr)/max_mag, cohpsk.ct); 457 if abs(max_corr/max_mag) > 0.9 458 printf(" [%d] encouraging sync word! ratio: %f\n", cohpsk.frame, abs(max_corr/max_mag)); 459 cohpsk.sync_timer = 0; 460 next_sync = 1; 461 else 462 next_sync = 0; 463 end 464 cohpsk.ratio = abs(max_corr/max_mag); 465 end 466 467 % single point correlation just to see if we are still in sync 468 469 if sync == 1 470 corr = 0; mag = 0; 471 f_fine_rect = exp(-j*cohpsk.f_fine_est*2*pi*sampling_points/Rs)'; 472 for c=1:Nc*Nd 473 f_corr_vec = f_fine_rect .* ct_symb_buf(cohpsk.ct+sampling_points,c); 474 acorr = 0; 475 for p=1:length(sampling_points) 476 acorr += pilot2(p, c-Nc*floor((c-1)/Nc)) * f_corr_vec(p); 477 mag += abs(f_corr_vec(p)); 478 end 479 corr += abs(acorr); 480 end 481 cohpsk.ratio = abs(corr)/mag; 482 end 483 484endfunction 485 486 487% misc sync state machine code, just wanted it in a function 488 489function [sync cohpsk] = sync_state_machine(cohpsk, sync, next_sync) 490 491 if sync == 1 492 493 % check that sync is still good, fall out of sync on consecutive bad frames */ 494 495 if cohpsk.ratio < 0.8 496 cohpsk.sync_timer++; 497 else 498 cohpsk.sync_timer = 0; 499 end 500 501 if cohpsk.sync_timer == 10 502 printf(" [%d] lost sync ....\n", cohpsk.frame); 503 next_sync = 0; 504 end 505 end 506 507 sync = next_sync; 508endfunction 509 510