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