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