1% hf_modem_curves 2% David Rowe Feb 2017 3% 4% Ideal implementations of a bunch of different HF modems, used to 5% generate plots for a blog post. 6 7#{ 8 [X] ideal AWGN/HF curves 9 [X] exp AWGN QPSK curves 10 [X] exp AWGN DQPSK curves 11 [X] exp HF channel model 12 [ ] diversity 13 [ ] COHPSK frames 14 + would require multiple carriers 15 + filtering or OFDM 16#} 17 181; 19 20% Gray coded QPSK modulation function 21 22function symbol = qpsk_mod(two_bits) 23 two_bits_decimal = sum(two_bits .* [2 1]); 24 switch(two_bits_decimal) 25 case (0) symbol = 1; 26 case (1) symbol = j; 27 case (2) symbol = -j; 28 case (3) symbol = -1; 29 endswitch 30endfunction 31 32 33% Gray coded QPSK demodulation function 34 35function two_bits = qpsk_demod(symbol) 36 bit0 = real(symbol*exp(j*pi/4)) < 0; 37 bit1 = imag(symbol*exp(j*pi/4)) < 0; 38 two_bits = [bit1 bit0]; 39endfunction 40 41 42% Rate Rs modem simulation model ------------------------------------------------------- 43 44function sim_out = ber_test(sim_in) 45 bps = 2; % two bits/symbol for QPSK 46 Rs = 50; % symbol rate (needed for HF model) 47 48 verbose = sim_in.verbose; 49 EbNovec = sim_in.EbNovec; 50 hf_en = sim_in.hf_en; 51 52 % user can supply number of bits per point to get good results 53 % at high Eb/No 54 55 if length(sim_in.nbits) > 1 56 nbitsvec = sim_in.nbits; 57 nbitsvec += 100 - mod(nbitsvec,100); % round up to nearest 100 58 else 59 nbitsvec(1:length(EbNovec)) = sim_in.nbits; 60 end 61 62 % init HF model 63 64 if hf_en 65 66 % some typical values 67 68 dopplerSpreadHz = 1.0; path_delay = 1E-3*Rs; 69 70 nsymb = max(nbitsvec)/2; 71 spread1 = doppler_spread(dopplerSpreadHz, Rs, nsymb); 72 spread2 = doppler_spread(dopplerSpreadHz, Rs, nsymb); 73 hf_gain = 1.0/sqrt(var(spread1)+var(spread2)); 74 % printf("nsymb: %d lspread1: %d\n", nsymb, length(spread1)); 75 end 76 77 for ne = 1:length(EbNovec) 78 79 % work out noise power ------------- 80 81 EbNodB = EbNovec(ne); 82 EsNodB = EbNodB + 10*log10(bps); 83 EsNo = 10^(EsNodB/10); 84 variance = 1/EsNo; 85 nbits = nbitsvec(ne); 86 nsymb = nbits/bps; 87 88 % modulator ------------------------ 89 90 tx_bits = rand(1,nbits) > 0.5; 91 tx_symb = []; 92 prev_tx_symb = 1; 93 for s=1:nsymb 94 atx_symb = qpsk_mod(tx_bits(2*s-1:2*s)); 95 if sim_in.dqpsk 96 atx_symb *= prev_tx_symb; 97 prev_tx_symb = atx_symb; 98 end 99 tx_symb = [tx_symb atx_symb]; 100 end 101 102 % channel --------------------------- 103 104 rx_symb = tx_symb; 105 106 if hf_en 107 108 % simplified rate Rs simulation model that doesn't include 109 % ISI, just freq filtering. We assume perfect phase estimation 110 % so it's just amplitude distortion. 111 112 hf_model1 = hf_model2 = zeros(1, nsymb); 113 for s=1:nsymb 114 hf_model1(s) = hf_gain*(spread1(s) + exp(-j*path_delay)*spread2(s)); 115 hf_model = abs(hf_model1(s)); 116 117 if sim_in.diversity 118 % include amplitude information from another frequency in channel model 119 w1 = 7*2*pi; 120 hf_model2(s) = hf_gain*(spread1(s) + exp(-j*w1*path_delay)*spread2(s)); 121 hf_model = 0.5*abs(hf_model1(s)) + 0.5*abs(hf_model2(s)); 122 end 123 124 rx_symb(s) = rx_symb(s).*hf_model; 125 end 126 end 127 128 % variance is noise power, which is divided equally between real and 129 % imag components of noise 130 131 noise = sqrt(variance*0.5)*(randn(1,nsymb) + j*randn(1,nsymb)); 132 rx_symb += noise; 133 134 % demodulator ------------------------------------------ 135 136 % demodulate rx symbols to bits 137 138 rx_bits = []; 139 prev_rx_symb = 1; 140 for s=1:nsymb 141 arx_symb = rx_symb(s); 142 if sim_in.dqpsk 143 tmp = arx_symb; 144 arx_symb *= prev_rx_symb'; 145 prev_rx_symb = tmp; 146 end 147 two_bits = qpsk_demod(arx_symb); 148 rx_bits = [rx_bits two_bits]; 149 end 150 151 % count errors ----------------------------------------- 152 153 error_pattern = xor(tx_bits, rx_bits); 154 nerrors = sum(error_pattern); 155 bervec(ne) = nerrors/nbits; 156 if verbose 157 printf("EbNodB: % 3.1f nbits: %5d nerrors: %5d ber: %4.3f\n", EbNodB, nbits, nerrors, bervec(ne)); 158 if verbose == 2 159 figure(2); clf; 160 plot(rx_symb*exp(j*pi/4),'+','markersize', 10); 161 mx = max(abs(rx_symb)); 162 axis([-mx mx -mx mx]); 163 if sim_in.diversity && sim_in.hf_en 164 figure(3); 165 plot(1:nsymb, abs(hf_model1), 1:nsymb, abs(hf_model2), 'linewidth', 2); 166 end 167 end 168 end 169 end 170 171 sim_out.bervec = bervec; 172endfunction 173 174 175% ------------------------------------------------------------- 176 177 178function run_single 179 sim_in.verbose = 2; 180 sim_in.nbits = 1000; 181 sim_in.EbNovec = 4; 182 sim_in.dqpsk = 0; 183 sim_in.hf_en = 0; 184 sim_in.diversity = 0; 185 186 sim_qpsk = ber_test(sim_in); 187endfunction 188 189 190function run_curves 191 max_nbits = 1E5; 192 sim_in.verbose = 1; 193 sim_in.EbNovec = 0:10; 194 sim_in.dqpsk = 0; 195 sim_in.hf_en = 0; 196 sim_in.diversity = 0; 197 198 % AWGN ----------------------------- 199 200 ber_awgn_theory = 0.5*erfc(sqrt(10.^(sim_in.EbNovec/10))); 201 sim_in.nbits = min(max_nbits, floor(500 ./ ber_awgn_theory)); 202 203 sim_qpsk = ber_test(sim_in); 204 sim_in.dqpsk = 1; 205 sim_dqpsk = ber_test(sim_in); 206 207 % HF ----------------------------- 208 209 hf_sim_in = sim_in; hf_sim_in.dqpsk = 0; hf_sim_in.hf_en = 1; 210 hf_sim_in.EbNovec = 0:16; 211 212 EbNoLin = 10.^(hf_sim_in.EbNovec/10); 213 ber_hf_theory = 0.5.*(1-sqrt(EbNoLin./(EbNoLin+1))); 214 215 hf_sim_in.nbits = min(max_nbits, floor(500 ./ ber_hf_theory)); 216 sim_qpsk_hf = ber_test(hf_sim_in); 217 218 hf_sim_in.dqpsk = 1; 219 sim_dqpsk_hf = ber_test(hf_sim_in); 220 221 hf_sim_in.dqpsk = 0; 222 hf_sim_in.diversity = 1; 223 sim_qpsk_hf_div = ber_test(hf_sim_in); 224 225 % Plot results -------------------- 226 227 close all; 228 figure (1, 'position', [100, 10, 600, 400]); clf; 229 230 semilogy(sim_in.EbNovec, ber_awgn_theory,'r+-;QPSK AWGN theory;', 'linewidth', 2) 231 xlabel('Eb/No (dB)') 232 ylabel('BER') 233 grid("minor") 234 axis([min(sim_in.EbNovec) max(sim_in.EbNovec) 1E-3 1]) 235 hold on; 236 237 semilogy([0 4 4], [ber_awgn_theory(5) ber_awgn_theory(5) 1E-3],'k--', 'linewidth', 2); 238 hold off; 239 240 figure (2, 'position', [300, 10, 600, 400]); clf; 241 semilogy(sim_in.EbNovec, ber_awgn_theory,'r+-;QPSK AWGN theory;','markersize', 10, 'linewidth', 2) 242 hold on; 243 semilogy(sim_in.EbNovec, sim_qpsk.bervec,'g+-;QPSK AWGN simulated;','markersize', 10, 'linewidth', 2) 244 semilogy(sim_in.EbNovec, sim_dqpsk.bervec,'b+-;DQPSK AWGN simulated;','markersize', 10, 'linewidth', 2) 245 xlabel('Eb/No (dB)') 246 ylabel('BER') 247 grid("minor") 248 axis([min(sim_in.EbNovec) max(sim_in.EbNovec) 1E-3 1]) 249 250 figure (3, 'position', [400, 10, 600, 400]); clf; 251 semilogy(sim_in.EbNovec, ber_awgn_theory,'r+-;QPSK AWGN theory;','markersize', 10, 'linewidth', 2) 252 hold on; 253 semilogy(sim_in.EbNovec, sim_qpsk.bervec,'g+-;QPSK AWGN simulated;','markersize', 10, 'linewidth', 2) 254 semilogy(sim_in.EbNovec, sim_dqpsk.bervec,'b+-;DQPSK AWGN simulated;','markersize', 10, 'linewidth', 2) 255 semilogy(hf_sim_in.EbNovec, ber_hf_theory,'r+-;QPSK HF theory;','markersize', 10, 'linewidth', 2) 256 semilogy(hf_sim_in.EbNovec, sim_dqpsk_hf.bervec,'b+-;DQPSK HF simulated;','markersize', 10, 'linewidth', 2) 257 semilogy(hf_sim_in.EbNovec, sim_qpsk_hf.bervec,'g+-;QPSK HF simulated;','markersize', 10, 'linewidth', 2) 258 semilogy(hf_sim_in.EbNovec, sim_qpsk_hf_div.bervec,'c+-;QPSK Diversity HF simulated;','markersize', 10, 'linewidth', 2) 259 hold off; 260 xlabel('Eb/No (dB)') 261 ylabel('BER') 262 grid("minor") 263 axis([min(hf_sim_in.EbNovec) max(hf_sim_in.EbNovec) 1E-3 1]) 264 265endfunction 266 267% ------------------------------------------------------------- 268 269more off; 270rand('seed',1); randn('seed', 1); 271run_curves 272#run_single 273