1% Test MFSK at symbol level in CML using RA LDPC codes, Bill, June 2020 2% Simulate in AWGM and plot BERs. Assumes that CML is installed! 3% 4% If required setup numb of codewords (Ncw) and channel bits/sym (bps, 1 to 4) 5% may also select FEC code with Ctype (1 to 3); use plt to for debug plots of LLR values 6 7% July 1 version allows M=8 and pads required vectors if required to makeup 3rd bit 8 9ldpc; 10 11%rand('seed',1); 12%randn('seed',1); 13format short g 14more off 15init_cml(); 16 17if exist('Ncw')==0, Ncw=100, end %setup defaults 18if exist('plt')==0, plt=0; end 19if exist('bps')==0, bps=4; end 20if exist('Ctype')==0, Ctype=1, end 21 22if Ctype==1 23 load H_256_768_22.txt; HRA = H_256_768_22; % rate 1/3 24elseif Ctype==2 25 load H_256_512_4.mat; HRA=H; % K=256, rate 1/2 code 26 % above code might be improved -- but still works better than rate 1/3 27elseif Ctype==3 28 load HRAa_1536_512.mat; % rate 3/4, N=2k code 29else 30 error('bad Ctype'); 31end 32 33M=2^bps; nos =0; clear res 34modulation = 'FSK'; mod_order=M; mapping = 'gray'; 35code_param = ldpc_init_user(HRA, modulation, mod_order, mapping); 36 37[Nr Nc] = size(HRA); 38Nbits = Nc - Nr; 39Krate = (Nc-Nr)/Nc 40framesize = Nc; 41%{ 42[H_rows, H_cols] = Mat2Hrows(HRA); 43code_param.H_rows = H_rows; 44code_param.H_cols = H_cols; 45code_param.P_matrix = []; 46%} 47 48 49S =CreateConstellation('FSK', M); 50if M==2, Ebvec=[7:0.2: 8.7], end 51if M==4, Ebvec=[4:0.25: 7], end 52if M==8, Ebvec =[3.5: 0.25: 5.5]; end 53if M==16, Ebvec=[2.5:0.25:4.8]; end 54 55disp(['Symbol-based ' num2str(M) 'FSK sim with K=' ... 56 num2str(Nbits) ', code rate=' num2str(Krate) ', #CWs=' num2str(Ncw)]) 57 58 59 60% if M=8, for 3 bits/symbol, may need to pad codeword with bits ... 61if floor(Nc/bps)*bps ~= Nc 62 Npad = ceil(Nc/bps) *bps-Nc 63 disp('padding codeword') 64else 65 Npad=0; 66end 67 68for Eb = Ebvec 69 70 Ec = Eb + 10*log10(Krate); 71 Es = Ec + 10*log10(bps); 72 Eslin = 10^(Es/10); %Es/N0 = 1/2k_n^2 73 74 Terrs =0; 75 for nn = 1:Ncw 76 77 txbits = randi(2,1,Nbits) -1; 78 79 codeword = LdpcEncode( txbits, code_param.H_rows, code_param.P_matrix ); 80 code_param.code_bits_per_frame = length( codeword ); 81 code_param.data_bits_per_frame = length(txbits); 82 Nsymb = (code_param.code_bits_per_frame+Npad)/bps; 83 84 if Npad; codeword = [codeword zeros(1,Npad)]; end 85 Tx = Modulate(codeword, S); 86 87 kn = sqrt(1/(2*Eslin)); 88 Rx = Tx + kn * (randn(size(Tx)) + j*randn(size(Tx))); 89 90 SNRlin = Eslin; % Valenti calls this SNR, but seems to be Es/N0 91 symL = DemodFSK(Rx, SNRlin, 2); %demod type is nonCOH, without estimate amplitudes 92 bitL = Somap(symL); 93 94 if Npad, bitL(end-Npad+1:end)=[]; end 95 if plt>0, figure(110); hist(bitL); title('bit LLRs') 96 figure(111); hist(bitL); title('Sym Ls'), pause, 97 end 98 max_it =100; decoder_type =0; 99 100 [x_hat, PCcnt] = MpDecode( -bitL, code_param.H_rows, code_param.H_cols, ... 101 max_it, decoder_type, 1, 1); 102 Niters = sum(PCcnt~=0); 103 detected_data = x_hat(Niters,:); 104 error_positions = xor( detected_data(1:code_param.data_bits_per_frame), txbits ); 105 106 Nerrs = sum( error_positions); 107 if plt>1, figure(121); plot(error_positions); Nerrs, end 108 Terrs = Terrs + Nerrs; 109 end 110 111 BER = Terrs/ (Ncw*Nbits); 112 113 %HDs = (sign(bitL)+1)/2; 114 %NerrsHD = sum(codeword~=HDs); 115 %BER_HD = Nerrs/Nbits; 116 117 nos = nos+1; 118 disp('Eb Nerrs BER') 119 res(nos, :) = [Eb, Terrs, BER] 120 121end 122figure(91) 123semilogy(res(:,1), res(:,3), '-x'); grid on; hold on; 124%semilogy(res(:,1), berawgn(res(:,1), 'fsk', M, 'noncoherent'), 'g'); 125title([num2str(M) 'FSK BER with LDPC FEC']) 126xlabel('Eb/N0'); ylabel('BER') 127 128 129 130 131 132 133