1% fsk_llr_test.m
2%
3% 2/4FSK simulation to develop LLR estimation algorithms for 4FSK/LDPC modems
4% Modified version of David's fsk_llr.m;   Bill
5
6#{
7  TODO
8  The 'v' param of the Ricean pdf is the signal-only amplitude: genie value=16
9  In practice, given varying input levels, this value needs to be estimated.
10
11  A small scaling factor seems to improve 2FSK performance -- probably the 'sig'
12  estimate can be improved.
13
14  Only tested with short code -- try a longer one!
15
16  Simulation should be updated to exit Eb after given Nerr reached
17
18#}
19
20ldpc;
21
22% define Rician pdf
23function y = rice(x,v,s)
24  s2 = s*s;
25  y = (x / s2) .* exp(-0.5 * (x.^2 + v.^2)/s2) .* besseli(0, x*v/s2);
26endfunction
27
28function plot_pdf(v,s)
29  x=(0:0.1:2*v);
30  y= rice(x, v, s);
31  figure(201);   hold on
32  plot(x,y,'g');
33  %title('Rician pdf: signal carrier')
34  y= rice(x, 0, s);
35  plot(x,y,'b');
36  title('Rician pdf: signal and noise-only carriers')
37  pause(0.01);
38endfunction
39
40% single Eb/No point simulation    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
41
42function [raw_ber rx_filt rx_bits tx_symbols demapper sig_est ] = run_single(tx_bits, M, EcNodB,  plt=0)
43  % Ec/N0 is per channel bit
44  bps = log2(M);  % bits per symbol
45  Ts = 16;        % length of each symbol in samples
46
47  if length(tx_bits)==1
48     Nbits = tx_bits;
49     tx_bits = randi(2,1,Nbits)-1;    % make random bits
50  endif
51
52  Nbits = length(tx_bits);
53  Nsymbols = Nbits/log2(M);
54  tx_symbols = zeros(1,Nsymbols);
55
56  mapper = bps:-1:1;
57  % look up table demapper from symbols to bits (hard decision)
58  demapper=zeros(M,bps);
59  for m=1:M
60    for b=1:bps
61      if  bitand(m-1,b) demapper(m,bps-b+1) = 1; end
62    end
63  end
64
65  % continuous phase mFSK modulator
66
67  w(1:M) = 2*pi*(1:M)/Ts;
68  tx_phase = 0;
69  tx = zeros(1,Ts*Nsymbols);
70
71  for s=1:Nsymbols
72    bits_for_this_symbol = tx_bits(bps*(s-1)+1:bps*s);
73    symbol_index = bits_for_this_symbol * mapper' + 1;
74    tx_symbols(s) = symbol_index;
75    assert(demapper(symbol_index,:) == bits_for_this_symbol);
76    for k=1:Ts
77      tx_phase += w(symbol_index);
78      tx((s-1)*Ts+k) = exp(j*tx_phase);
79    end
80  end
81
82  % AWGN channel noise
83
84
85  EsNodB = EcNodB + 10*log10(bps)
86  EsNo = 10^(EsNodB/10);
87  variance = Ts/EsNo;
88  noise = sqrt(variance/2)*(randn(1,Nsymbols*Ts) + j*randn(1,Nsymbols*Ts));
89  rx = tx + noise;
90
91  if plt==2,    % check the Spectrum
92    [psd,Fpsd] =pwelch(rx,128,0.5,128,Ts);
93    figure(110); plot(Fpsd,10*log10(psd));
94    title('Rx Signal:  PSD ');
95    xlabel('Freq/Rs');
96    %figure(111);plot(unwrap(arg(tx)));
97    pause(0.01);
98  endif
99
100
101  % integrate and dump demodulator
102
103  rx_bits = zeros(1,Nbits);
104  rx_filt = zeros(Nsymbols,M);
105  rx_pows = zeros(1,M);
106  rx_nse_pow  = 0.0; rx_sig_pow =0.0;
107  for s=1:Nsymbols
108    arx_symb = rx((s-1)*Ts + (1:Ts));
109    for m=1:M
110      r= sum(exp(-j*w(m)*(1:Ts)) .* arx_symb);
111      rx_pows(m)= r * conj(r);
112      rx_filt(s,m) = abs(r);
113    end
114    [tmp symbol_index] = max(rx_filt(s,:));
115    rx_sig_pow = rx_sig_pow + rx_pows(symbol_index);
116    rx_pows(symbol_index)=[];
117    rx_nse_pow = rx_nse_pow + sum(rx_pows)/(M-1);
118    rx_bits(bps*(s-1)+1:bps*s) = demapper(symbol_index,:);
119  end
120  % using Rxpower = v^2 + sigmal^2
121
122  rx_sig_pow = rx_sig_pow/Nsymbols;
123  rx_nse_pow = rx_nse_pow/Nsymbols;
124  sig_est = sqrt(rx_nse_pow/2)    % for Rayleigh: 2nd raw moment = 2 .sigma^2
125  Kest = rx_sig_pow/(2.0*sig_est^2) -1.0
126
127  Nerrors = sum(xor(tx_bits, rx_bits));
128  raw_ber = Nerrors/Nbits;
129  printf("EcNodB: %4.1f  M: %2d Uncoded Nbits: %5d Nerrors: %4d (Raw) BER: %1.3f\n", ...
130          EcNodB, M, Nbits, Nerrors, raw_ber);
131  if plt==1, plot_hist(rx_filt,tx_symbols, M); end
132
133endfunction
134
135
136% Plot histograms of Rx filter outputs
137function plot_hist(rx_filt,tx_symbols, M)
138   % more general version of previous fn; - plots histograms for any Tx patterns
139   Smax = 36;
140   X = 0:Smax-1;
141   H = zeros(1,Smax);  H2 = zeros(1,Smax); s2=0.0;
142   for m = 1:M
143     ind = tx_symbols==m;
144     ind2 =  tx_symbols~=m;
145     H= H+ hist(rx_filt(ind,m),X);
146     H2= H2+ hist(rx_filt(ind2,m),X);
147     x=rx_filt(ind2,m);
148     s2 =s2 + sum(x(:).^2)/length(x);
149   end
150   disp('noise RMS is ');  sqrt(s2/4)
151   figure(1);  clf; plot(X,H);
152   title([num2str(M) 'FSK pdf for rx=tx symbol'])
153   figure(2); clf;  plot(X,H2);
154   title([num2str(M) 'FSK pdf for rx!=tx symbol'])
155   pause(0.1);
156
157endfunction
158
159% 2FSK SD -> LLR mapping that we used for Wenet SSTV system
160function llr = sd_to_llr(sd, HD=0)     %    original 2FSK  + HD option
161    sd = sd / mean(abs(sd));
162    x = sd - sign(sd);
163    sumsq = sum(x.^2);
164    summ = sum(x);
165    mn = summ/length(sd);
166    estvar = sumsq/length(sd) - mn*mn;
167    estEsN0 = 1/(2* estvar + 1E-3);
168    if HD==0,
169      llr = 4 * estEsN0 * sd;
170    else
171       llr = 4 * estEsN0 * sign(sd);
172    endif
173endfunction
174
175
176% single point LDPC encoded frame simulation, usin 2FSK as a tractable starting point
177% Note: ~/cml/matCreateConstellation.m has some support for FSK - can it do 4FSK?
178
179%%%%%%%%%%%%%%%%%%%%%%%%%
180  function [Nerrors raw_ber EcNodB] = run_single_ldpc(M, Ltype, Nbits,EbNodB, plt=0)
181
182  disp([num2str(M) 'FSK coded test ...  '])
183  if M==2
184    bps = 1; modulation = 'FSK'; mod_order=2; mapping = 'gray';
185  elseif M==4
186    bps = 2; modulation = 'FSK'; mod_order=4; mapping = 'gray';
187  else
188    error('sorry - bad value of M!');
189  endif
190  decoder_type = 0; max_iterations = 100;
191
192  load H_256_768_22.txt
193  Krate = 1/3;
194  EcNodB = EbNodB + 10*log10(Krate);
195  code_param = ldpc_init_user(H_256_768_22, modulation, mod_order, mapping);
196  Nframes = floor(Nbits/code_param.data_bits_per_frame)
197  Nbits = Nframes*code_param.data_bits_per_frame
198
199  % Encoder
200  data_bits = round(rand(1,code_param.data_bits_per_frame));
201  tx_bits = [];
202  for f=1:Nframes;
203    codeword_bits = LdpcEncode(data_bits, code_param.H_rows, code_param.P_matrix);
204    tx_bits = [tx_bits codeword_bits];
205  end
206  %tx_bits = zeros(1,length(tx_bits));
207
208  % modem/channel simulation
209  [raw_ber rx_filt rx_bits tx_symbols demapper sig_est ] = run_single(tx_bits,M,EcNodB, 0 );
210
211  % Decoder
212  Nerrors = 0;
213  for f=1:Nframes
214    st = (f-1)*code_param.coded_bits_per_frame/bps + 1;
215    en = st + code_param.coded_bits_per_frame/bps - 1;
216
217    if or(Ltype==1, Ltype==3)
218        if bps==1,
219	  sd = rx_filt(st:en,1) - rx_filt(st:en,2);
220          %  OR ind = rx_filt(st:en,1) > rx_filt(st:en,2);
221          %     llr = ind'*2 -1;   % this works but use SNR scaling
222          if Ltype==3, HD=1; else, HD = 0;  endif
223          llr = sd_to_llr(sd, HD)';
224        endif
225        if bps==2,
226           if Ltype==3,
227	  llr = mfsk_hd_to_llrs(rx_filt(st:en,:), demapper);
228	   else
229              error('Ltype =1 not provided for coded 4FSK');
230           endif
231         endif
232    endif
233    if Ltype==2,     % SDs are converted to LLRs
234       v=16;
235       if plt==1, plot_pdf(v, sig_est);  endif
236       llr = mfsk_sd_to_llrs(rx_filt(st:en,:), demapper, v, sig_est);
237    endif
238
239    [x_hat, PCcnt] = MpDecode(llr, code_param.H_rows, code_param.H_cols, ...
240                            max_iterations, decoder_type, 1, 1);
241    Niters = sum(PCcnt!=0);
242    detected_data = x_hat(Niters,:);
243    Nerrors += sum(xor(data_bits, detected_data(1:code_param.data_bits_per_frame)));
244  endfor
245  ber = Nerrors/Nbits;
246  printf("EbNodB: %4.1f  Coded   Nbits: %5d Nerrors: %4d BER: %1.3f\n", EbNodB, Nbits, Nerrors, ber);
247endfunction
248
249%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
250
251%rand('seed',1);
252%randn('seed',1);
253format short
254more off
255init_cml();
256
257% store results in array "res" and plot afterwards
258% comment the following line if you want to retain prev sims
259nrun = 0; clear res;
260
261Nbits = 20000;  plt=0;
262
263#{
264disp(' uncoded runs')
265for M=  [2 4]
266for Eb = [6:10]
267  raw_ber = run_single(Nbits, M,  Eb, plt)    % 2fsk coded
268   nrun = nrun+1; res(nrun,:) =  [Eb Eb M -1 Nbits -1  raw_ber]
269endfor
270endfor
271#}
272
273disp(' coded runs ');
274
275M=2,
276for Ltype = [1 2 3]
277for Eb = [7: 0.5: 9]
278  [Nerr raw_ber Ec] = run_single_ldpc(M, Ltype, Nbits, Eb, plt)
279   nrun = nrun+1; res(nrun,:) =  [Eb Ec M Ltype Nbits Nerr raw_ber]
280endfor
281endfor
282
283M=4,   %v=16;
284for Ltype = [2 3]
285for Eb = [8.0 8.3 8.6 ]
286  [Nerr raw_ber Ec] = run_single_ldpc(M, Ltype, Nbits, Eb, plt)
287   nrun = nrun+1; res(nrun,:) =  [Eb Ec M Ltype Nbits Nerr raw_ber]
288endfor
289endfor
290
291
292date = datestr(now)
293save 'mfsk_test_res.mat'  res date
294
295