1%fsk_llr_sam Test MFSK using David's sample-based mod/demod with CML LLR routines 2% using 2k rate 3/4 RA LDPC code. Bill, July 2020 3% 4%LLR conversion options: (Ltype) 5% 1 David's original M=2 SD to LLR routine 6% 2 Bill's pdf-based M=2 or 4 SD to LLR 7% 3 Some simple HD conversions 8% 4 CML approach using symbol likelihoods, then converted to bit LLRs 9% 10% Results from this sim are stored in "res" -- use fsk_llr_plot to see BER figs 11% Use "plt" to see some useful plots (for selected Ltypes) : 12% eg 1 is SD pdf histograms; 2 is Rx PSD; 3 is bit LLR histograms 13% 14% Adjust Evec and Nbits as required before running. 15 16#{ 17 Example 1 2FSK: 18 octave:4> fsk_cml_sam 19 octave:5> fsk_llr_plot 20 21 Example 2 4FSK: 22 octave:4> M=4; fsk_cml_sam 23 octave:4> fsk_llr_plot 24#} 25 26ldpc; 27 28% define Rician pdf 29% note that Valenti uses an approximation that avoids Bessel evaluation 30function y = rice(x,v,s) 31 s2 = s*s; 32 y = (x / s2) .* exp(-0.5 * (x.^2 + v.^2)/s2) .* besseli(0, x*v/s2); 33endfunction 34 35function plot_pdf(v,s) 36 x=(0:0.1:2*v); 37 y= rice(x, v, s); 38 figure(201); hold on 39 plot(x,y,'g'); 40 %title('Rician pdf: signal carrier') 41 y= rice(x, 0, s); 42 plot(x,y,'b'); 43 title('Rician pdf: signal and noise-only carriers') 44 pause(0.01); 45endfunction 46 47% single Eb/No point simulation --------------- 48function [raw_ber rx_filt rx_bits tx_symbols demapper sig_est SNRest v_est] = ... 49 run_single(tx_bits, M, EcNodB, plt) 50 % Ec/N0 is per channel bit 51 bps = log2(M); % bits per symbol 52 Ts = 16; % length of each symbol in samples 53 54 if length(tx_bits)==1 55 Nbits = tx_bits; 56 tx_bits = randi(2,1,Nbits)-1; % make random bits 57 end 58 59 Nbits = length(tx_bits); 60 Nsymbols = Nbits/log2(M); 61 tx_symbols = zeros(1,Nsymbols); 62 63 mapper = bps:-1:1; 64 % look up table demapper from symbols to bits (hard decision) 65 demapper=zeros(M,bps); 66 for m=1:M 67 for b=1:bps 68 if bitand(m-1,b) demapper(m,bps-b+1) = 1; end 69 end 70 end 71 72 % continuous phase mFSK modulator 73 74 w(1:M) = 2*pi*(1:M)/Ts; 75 tx_phase = 0; 76 tx = zeros(1,Ts*Nsymbols); 77 78 for s=1:Nsymbols 79 bits_for_this_symbol = tx_bits(bps*(s-1)+1:bps*s); 80 symbol_index = bits_for_this_symbol * mapper' + 1; 81 tx_symbols(s) = symbol_index; 82 assert(demapper(symbol_index,:) == bits_for_this_symbol); 83 for k=1:Ts 84 tx_phase = tx_phase + w(symbol_index); 85 tx((s-1)*Ts+k) = exp(j*tx_phase); 86 end 87 end 88 89 % AWGN channel noise 90 91 EsNodB = EcNodB + 10*log10(bps); 92 EsNo = 10^(EsNodB/10); 93 variance = Ts/EsNo; 94 noise = sqrt(variance/2)*(randn(1,Nsymbols*Ts) + j*randn(1,Nsymbols*Ts)); 95 rx = tx + noise; 96 97 if plt==2, % check the Spectrum 98 [psd,Fpsd] =pwelch(rx,128,0.5,128,Ts); 99 figure(110); plot(Fpsd,10*log10(psd)); 100 title('Rx Signal: PSD '); 101 xlabel('Freq/Rs'); 102 %figure(111);plot(unwrap(arg(tx))); 103 pause(0.01); 104 end 105 106 107 % integrate and dump demodulator 108 109 rx_bits = zeros(1,Nbits); 110 rx_filt = zeros(Nsymbols,M); 111 rx_pows = zeros(1,M); 112 rx_nse_pow = 0.0; rx_sig_pow =0.0; 113 for s=1:Nsymbols 114 arx_symb = rx((s-1)*Ts + (1:Ts)); 115 for m=1:M 116 r= sum(exp(-j*w(m)*(1:Ts)) .* arx_symb); 117 rx_pows(m)= r * conj(r); 118 rx_filt(s,m) = abs(r); 119 end 120 [tmp symbol_index] = max(rx_filt(s,:)); 121 rx_sig_pow = rx_sig_pow + rx_pows(symbol_index); 122 rx_pows(symbol_index)=[]; 123 rx_nse_pow = rx_nse_pow + sum(rx_pows)/(M-1); 124 rx_bits(bps*(s-1)+1:bps*s) = demapper(symbol_index,:); 125 end 126 127 % using Rxpower = v^2 + sigmal^2 128 rx_sig_pow = rx_sig_pow/Nsymbols; 129 rx_nse_pow = rx_nse_pow/Nsymbols; 130 v_est = sqrt(rx_sig_pow-rx_nse_pow); 131 SNRest = rx_sig_pow/rx_nse_pow; 132 sig_est = sqrt(rx_nse_pow/2); % for Rayleigh: 2nd raw moment = 2 .sigma^2 133 Kest = rx_sig_pow/(2.0*sig_est^2) -1.0; 134 135 Nerrors = sum(xor(tx_bits, rx_bits)); 136 raw_ber = Nerrors/Nbits; 137 printf('Ec (dB): %4.1f M: %2d Total bits: %5d HD Errs: %6d (Raw) BER: %1.3f\n', ... 138 EcNodB, M, Nbits, Nerrors, raw_ber); 139 if plt==1, plot_hist(rx_filt,tx_symbols, M); end 140 141endfunction 142 143% simulate one codeword of 2 or 4 FSK -------------------------- 144function [Nerrors raw_ber EcNodB] = run_single_ldpc(M, Ltype, Nbits,EbNodB, plt, HRA) 145 146 disp([num2str(M) 'FSK coded test ... ']) 147 if M==2 148 bps = 1; modulation = 'FSK'; mod_order=2; mapping = 'gray'; 149 elseif M==4 150 bps = 2; modulation = 'FSK'; mod_order=4; mapping = 'gray'; 151 else 152 error('sorry - bad value of M!'); 153 end 154 decoder_type = 0; max_iterations = 100; 155 156 Hsize=size(HRA); 157 Krate = (Hsize(2)-Hsize(1))/Hsize(2); % 158 EcNodB = EbNodB + 10*log10(Krate); 159 code_param = ldpc_init_user(HRA, modulation, mod_order, mapping); 160 Nframes = floor(Nbits/code_param.data_bits_per_frame); 161 Nbits = Nframes*code_param.data_bits_per_frame; 162 163 % Encoder 164 data_bits = round(rand(1,code_param.data_bits_per_frame)); 165 tx_bits = []; 166 for f=1:Nframes; 167 codeword_bits = LdpcEncode(data_bits, code_param.H_rows, code_param.P_matrix); 168 tx_bits = [tx_bits codeword_bits]; 169 end 170 171 % modem/channel simulation 172 [raw_ber rx_filt rx_bits tx_symbols demapper sig_est SNRlin v_est] = ... 173 run_single(tx_bits,M,EcNodB, plt ); 174 175 % Decoder 176 Nerrors = 0; 177 for f=1:Nframes 178 st = (f-1)*code_param.coded_bits_per_frame/bps + 1; 179 en = st + code_param.coded_bits_per_frame/bps - 1; 180 if or(Ltype==1, Ltype==3) 181 if bps==1, 182 sd = rx_filt(st:en,1) - rx_filt(st:en,2); 183 % OR ind = rx_filt(st:en,1) > rx_filt(st:en,2); 184 % llr = ind'*2 -1; % this works but use SNR scaling 185 if Ltype==3, HD=1; else, HD = 0; end 186 llr = sd_to_llr(sd, HD)'; % David's orig 2FSK routine 187 end 188 if bps==2, 189 if Ltype==3, 190 llr = mfsk_hd_to_llrs(rx_filt(st:en,:), demapper); 191 else 192 error('Ltype =1 not provided for coded 4FSK'); 193 end 194 end 195 end 196 if Ltype==2, % SDs are converted to LLRs 197 % use the SD amp estimates, try Bill's new LLR routine 198 if plt==1, plot_pdf(v_est, sig_est); end 199 llr = mfsk_sd_to_llrs(rx_filt(st:en,:), demapper, v_est, sig_est); 200 end 201 if Ltype==4, 202 % use CML demod: non-coherent; still seems to need amplitude estimate 203 symL = DemodFSK(1/v_est*rx_filt(st:en,:)', SNRlin, 1); 204 llr = -Somap(symL); % now convert to bit LLRs 205 end 206 if plt==3, figure(204); hist(llr); 207 title('Histogram LLR for decoder inputs'); pause(0.01); end 208 209 [x_hat, PCcnt] = MpDecode(llr, code_param.H_rows, code_param.H_cols, ... 210 max_iterations, decoder_type, 1, 1); 211 Niters = sum(PCcnt~=0); 212 detected_data = x_hat(Niters,:); 213 Nerrors = Nerrors + sum(xor(data_bits, detected_data(1:code_param.data_bits_per_frame))); 214 end 215 216 ber = Nerrors/Nbits; 217 printf('Eb (dB): %4.1f Data Bits: %4d Num CWs: %5d Tot Errs: %5d (Cod) BER: %1.3f\n', ... 218 EbNodB,Nbits , Nframes, Nerrors, ber); 219endfunction 220 221function plot_hist(rx_filt,tx_symbols, M) 222 % more general version of previous fn; - plots histograms for any Tx patterns 223 Smax = 36; 224 X = 0:Smax-1; 225 H = zeros(1,Smax); H2 = zeros(1,Smax); s2=0.0; 226 for m = 1:M 227 ind = tx_symbols==m; 228 ind2 = tx_symbols~=m; 229 H= H+ hist(rx_filt(ind,m),X); 230 H2= H2+ hist(rx_filt(ind2,m),X); 231 x=rx_filt(ind2,m); 232 s2 =s2 + sum(x(:).^2)/length(x); 233 end 234 disp('noise RMS is '); sqrt(s2/4) 235 figure(207); clf, plot(X,H); 236 title([num2str(M) 'FSK pdf for rx=tx symbol']) 237 hold on, plot(X,H2,'g'); 238 title([num2str(M) 'FSK pdf for rx!=tx symbol']) 239 pause(0.1); 240endfunction 241 242 243% 2FSK SD -> LLR mapping that we used for Wenet SSTV system 244function llr = sd_to_llr(sd, HD=0) % original 2FSK + HD option 245 sd = sd / mean(abs(sd)); 246 x = sd - sign(sd); 247 sumsq = sum(x.^2); 248 summ = sum(x); 249 mn = summ/length(sd); 250 estvar = sumsq/length(sd) - mn*mn; 251 estEsN0 = 1/(2* estvar + 1E-3); 252 if HD==0, 253 llr = 4 * estEsN0 * sd; 254 else 255 llr = 8 * estEsN0 * sign(sd); %%%% *4 >> *8 256 endif 257endfunction 258 259 260 261function llrs = mfsk_sd_to_llrs(SD, map, v, sig) 262 % Array of MFSK SD is converted to bit LLR vector according to mapping 'map' 263 % SD is M elements wide; map is M by log2(M); v and sig are params of Rician pdf 264 % Bill (VK5DSP) for Codec2, version 24/6/20 265 266 Llim = 1e-7; 267 XX =1.0; % check LLR scaling? 268 Smap = size(map); 269 Ssd = size(SD); 270 if Smap(1) == 2 271 llrs = zeros(1, Ssd(1)); 272 bps = 1; 273 assert(2^bps == Ssd(2), 'wrong SD length?') 274 % note that when v=0, the pdf reverts to Rayleigh 275 for kk = 1: Ssd(1) 276 Prx_1 = rice(SD(kk,2), v, sig) * rice(SD(kk,1),0,sig); 277 Prx_0 = rice(SD(kk,2), 0, sig) * rice(SD(kk,1),v,sig); 278 279 if or(isnan(Prx_1), Prx_1<Llim), Prx_1 = Llim; end 280 if or(isnan(Prx_0), Prx_0<Llim), Prx_0 = Llim; end 281 llrs(kk) = -XX * log(Prx_1/Prx_0); 282 %% note inversion of LLRs 283 endfor 284 else 285 if map==[0 0; 0 1; 1 0; 1 1] % a specific case for 4FSK 286 llrs = zeros(2, Ssd(1)); 287 for kk = 1: Ssd(1) 288 % pn_m is prob of sub car n given bit m is transmitted 289 p1_0 = rice(SD(kk,1), 0, sig); p1_1 = rice(SD(kk,1), v, sig); 290 p2_0 = rice(SD(kk,2), 0, sig); p2_1 = rice(SD(kk,2), v, sig); 291 p3_0 = rice(SD(kk,3), 0, sig); p3_1 = rice(SD(kk,3), v, sig); 292 p4_0 = rice(SD(kk,4), 0, sig); p4_1 = rice(SD(kk,4), v, sig); 293 294 % second bit 295 Prx_1 = p1_0*p3_0*(p2_1*p4_0 + p2_0*p4_1); 296 Prx_0 = p2_0*p4_0*(p1_1*p3_0 + p1_0*p3_1); 297 298 if or(isnan(Prx_1), Prx_1<Llim), Prx_1 = Llim; end 299 if or(isnan(Prx_0), Prx_0<Llim), Prx_0 = Llim; end 300 llrs(2,kk) = -XX*log(Prx_1/Prx_0); 301 302 % 1st bit (LHS) 303 Prx_1 = p1_0*p2_0*(p3_1*p4_0 + p3_0*p4_1); 304 Prx_0 = p3_0*p4_0*(p1_1*p2_0 + p1_0*p2_1); 305 306 if or(isnan(Prx_1), Prx_1<Llim), Prx_1 = Llim; end 307 if or(isnan(Prx_0), Prx_0<Llim), Prx_0 = Llim; end 308 llrs(1,kk) = -XX* log(Prx_1/Prx_0); 309 endfor 310 llrs= llrs(:); % convert to single vector of LLRs 311 llrs = llrs'; 312 else 313 disp('the mapping is '); map 314 error('not sure how that mapping works!!') 315 endif 316 endif 317endfunction 318 319function llrs = mfsk_hd_to_llrs(sd, map, SNRlin=4) 320% for M=4, compare these results to SD method of mfsk_sd_to_llrs 321 [x, ind] = max(sd'); % find biggest SD 322 b2 = map(ind', :)'; % get bit pairs 323 b1 = b2(:)'; % covert to row vector 324 b1 = b1*2 -1; % convert to -1/ +1 325 llrs = -4*SNRlin*b1; % approx scaling; default is ~6dB SNR 326endfunction 327 328 329% main ------------------------------------------------------------------------------ 330 331format short 332more off 333init_cml(); 334 335if exist('Ctype')==0, Ctype=1, end 336 337if Ctype==1 338 load H_256_768_22.txt; HRA = H_256_768_22; % rate 1/3 339elseif Ctype==2 340 load H_256_512_4.mat; HRA=H; % rate 1/2 code 341elseif Ctype==3 342 load HRAa_1536_512.mat; % rate 3/4 2k code 343else 344 error('bad Ctype'); 345end 346 347Hsize=size(HRA); 348printf('using LDPC code length: %5d: Code rate: %5.2f\n',length(HRA), (Hsize(2)-Hsize(1))/Hsize(2)) 349 350% store results in array "res" and plot afterwards 351% retain prev sims? 352if exist('keep_res')==0, nrun = 0; clear res; end 353 354if exist('Nbits')==0, Nbits = 20000, end 355if exist('Evec')==0, Evec = [5: 0.5: 9], end 356if exist('plt')==0. plt=0; end 357 358if exist('M')==0, M=2, end 359 360for Ltype = [ 2 4] % << add other Ltype options here, if required 361 for Eb = Evec 362 if and(M==4, Ltype==1) 363 disp('skip original SD >> LLRs in 4FSK') 364 else 365 [Nerr raw_ber Ec] = run_single_ldpc(M, Ltype, Nbits, Eb, plt, HRA); 366 nrun = nrun+1; res(nrun,:) = [Eb Ec M Ltype Nbits Nerr raw_ber]; 367 endif 368 end 369end 370disp('results stored in array "res" - plot with fsk_cml_plot') 371disp(' Eb Ec M Ltype #Bits #Errs Raw-BER') 372for n = 1: nrun 373 printf(' %5.1f %5.1f %3d %3d %6d %6d %5.4f \n', res(n,:)) 374end 375 376 377