added LDPC curve for comparison, and imported trellis results

Author: drowe67

octave/vq_compare.m

@@ -27,7 +27,8 @@
 function vq_compare(action="run_curves", vq_fn, dec=1, EbNodB=3, in_fn, out_fn)
   more off;
   randn('state',1);
-
+  graphics_toolkit("gnuplot");
+  
   if strcmp(action, "run_curves")
     run_curves(30*100);
   end
@@ -123,7 +124,7 @@ function vq_compare(action="run_curves", vq_fn, dec=1, EbNodB=3, in_fn, out_fn)
   target_ = vq(rx_indexes+1,:);
   diff = target - target_;
   mse = mean(diff(:).^2);
-  printf("Eb/No: %3.2f dB nframes: %2d nerrors: %3d BER: %4.3f PER: %3.2f mse: %3.2f %3.2f\n", 
+  printf("Eb/No: %3.2f dB nframes: %3d nerrors: %4d BER: %4.3f PER: %3.2f mse: %3.2f %3.2f\n", 
          EbNodB, nframes, nerrors, nerrors/tbits, nper/nframes, mse_noerrors, mse);
   results.ber = nerrors/tbits;	 
   results.per = nper/nframes;	 
@@ -133,10 +134,113 @@ function vq_compare(action="run_curves", vq_fn, dec=1, EbNodB=3, in_fn, out_fn)
   results.rx_indexes = rx_indexes;
 endfunction
 
+% VQ a target sequence of frames then run a test using a LDPC code
+function results = run_test_ldpc(target, vq, EbNo, verbose)
+  [frames tmp]      = size(target);
+  [vq_length tmp]   = size(vq);
+  nbits             = log2(vq_length);
+  nerrors           = 0;
+  tbits             = 0;
+  nframes           = 0;
+  nper              = 0;
+  
+  % init LDPC code
+  mod_order = 4; bps = 2;
+  modulation = 'QPSK';
+  mapping = 'gray';
+  max_iterations = 100; demod_type = 0; decoder_type = 0;
+  ldpc; init_cml('~/cml/');
+  tempStruct = load("HRA_56_56.txt");
+  b = fieldnames(tempStruct);
+  ldpcArrayName = b{1,1};
+  % extract the array from the struct
+  HRA = tempStruct.(ldpcArrayName);
+  [code_param framesize rate] = ldpc_init_user(HRA, modulation, mod_order, mapping);
+
+  % set up noise
+  EbNodB = 10*log10(EbNo);
+  EsNodB = EbNodB + 10*log10(rate) + 10*log10(bps);
+  EsNo = 10^(EsNodB/10);
+  variance = 1/EsNo;
+
+  % Vector Quantise target vectors sequence
+  [tx_indexes target_ ] = vector_quantiser_fast(vq, target, verbose);
+  % use convention of indexes starting from 0
+  tx_indexes -= 1; 
+  %  mean SD of VQ with no errors
+  diff = target - target_;
+  mse_noerrors = mean(diff(:).^2);
+
+  % construct tx frames x nbit matrix using VQ indexes
+  tx_bits = zeros(frames, nbits);
+  for f=1:frames
+    tx_bits(f,:) = dec2sd(tx_indexes(f), nbits) > 0;
+  end
+
+  % find a superframe size, that has an integer number of nbits and data_bits_per_frame frames
+  bits_per_superframe = nbits;
+  while mod(bits_per_superframe,nbits) || mod(bits_per_superframe,code_param.data_bits_per_frame)
+    bits_per_superframe += nbits;
+  end
+
+  Nsuperframes = floor(frames*nbits/bits_per_superframe);
+  Nldpc_codewords = Nsuperframes*bits_per_superframe/code_param.data_bits_per_frame;
+  frames = Nsuperframes*bits_per_superframe/nbits;
+  %printf("bits_per_superframe: %d Nldpc_codewords: %d frames: %d\n", bits_per_superframe, Nldpc_codewords, frames);
+  
+  % reshape tx_bits matrix into Nldpc_codewords x data_bits_per_frame
+  tx_bits = tx_bits(1:frames,:);
+  tx_bits_ldpc = reshape(tx_bits',code_param.data_bits_per_frame, Nldpc_codewords)';
+
+  % modulate tx symbols
+  tx_symbols = [];
+  for nn=1:Nldpc_codewords
+    [tx_codeword atx_symbols] = ldpc_enc(tx_bits_ldpc(nn,:), code_param);
+    tx_symbols = [tx_symbols atx_symbols];
+  end
+  
+  noise = sqrt(variance*0.5)*(randn(1,length(tx_symbols)) + j*randn(1,length(tx_symbols)));
+  rx_symbols = tx_symbols+noise;
+  
+  % LDPC decode  
+  for nn = 1:Nldpc_codewords
+    st = (nn-1)*code_param.coded_syms_per_frame + 1;
+    en = (nn)*code_param.coded_syms_per_frame;
+
+    arx_codeword = ldpc_dec(code_param, max_iterations, demod_type, decoder_type, rx_symbols(st:en), EsNo, ones(1,code_param.coded_syms_per_frame));
+    rx_bits_ldpc(nn,:) = arx_codeword(1:code_param.data_bits_per_frame);
+  end
+  
+  % reshape rx_bits_ldpc matrix into frames x nbits
+  rx_bits = reshape(rx_bits_ldpc',nbits,frames)';
+
+  rx_indexes = tx_indexes;
+  for f=1:frames
+    rx_indexes(f)  = sum(rx_bits(f,:) .* 2.^(nbits-1:-1:0));
+    errors         = sum(xor(tx_bits(f,:), rx_bits(f,:)));
+    nerrors       += errors;
+    if errors nper++;, end
+    tbits += nbits;
+    nframes++;
+  end
+
+  EbNodB = 10*log10(EbNo);
+  target_ = vq(rx_indexes+1,:);
+  diff = target - target_;
+  mse = mean(diff(:).^2);
+  printf("Eb/No: %3.2f dB nframes: %4d nerrors: %4d BER: %4.3f PER: %3.2f mse: %3.2f %3.2f\n", 
+         EbNodB, nframes, nerrors, nerrors/tbits, nper/nframes, mse_noerrors, mse);
+  results.ber = nerrors/tbits;	 
+  results.per = nper/nframes;	 
+  results.mse = mse;
+  results.tx_indexes = tx_indexes;
+  results.rx_indexes = rx_indexes;
+endfunction
+
 % Simulations ---------------------------------------------------------------------
 
 % top level function to set up and run a test with a specific vq
-function [results target_] = run_test_vq(vq_fn, target_fn, nframes=100, dec=1, EbNodB=3, verbose=0)
+function [results target_] = run_test_vq(vq_fn, target_fn, nframes=100, dec=1, EbNodB=3, ldpc_en=0, verbose=0)
   K = 20; K_st=2+1; K_en=16+1;
 
   % load VQ
@@ -157,7 +261,11 @@ function vq_compare(action="run_curves", vq_fn, dec=1, EbNodB=3, in_fn, out_fn)
 
   % run a test
   EbNo=10^(EbNodB/10);
-  results = run_test(target, vqsub, EbNo, verbose);
+  if ldpc_en
+    results = run_test_ldpc(target, vqsub, EbNo, verbose);
+  else
+    results = run_test(target, vqsub, EbNo, verbose);
+  end
   if verbose
     for f=2:nframes-1
       printf("f: %03d tx_index: %04d rx_index: %04d\n", f,  results.tx_indexes(f), results.rx_indexes(f));
@@ -182,43 +290,53 @@ function vq_compare(action="run_curves", vq_fn, dec=1, EbNodB=3, in_fn, out_fn)
 % generate sets of curves
 function run_curves(frames=100, dec=1)
   target_fn = "../build_linux/all_speech_8k_lim.f32";
-  results1_log = [];
   EbNodB = 0:5;
+
+  results1_ldpc_log = [];
   for i=1:length(EbNodB)
-    results = run_test_vq("../build_linux/vq_stage1.f32", target_fn, frames, dec, EbNodB(i), verbose=0);
-    results1_log = [results1_log results];
+    results = run_test_vq("../build_linux/vq_stage1.f32", target_fn, frames, dec, EbNodB(i), ldpc_en=1, verbose=0);
+    results1_ldpc_log = [results1_ldpc_log results];
   end
-  results2_log = [];
+  results4_ldpc_log = [];
   for i=1:length(EbNodB)
-    results = run_test_vq("../build_linux/vq_stage1_bs001.f32", target_fn, frames, dec, EbNodB(i), verbose=0);
-    results2_log = [results2_log results];
+    results = run_test_vq("../build_linux/vq_stage1_bs004.f32", target_fn, frames, dec, EbNodB(i), ldpc_en=1, verbose=0);
+    results4_ldpc_log = [results4_ldpc_log results];
+  end
+
+  results1_log = [];
+  for i=1:length(EbNodB)
+    results = run_test_vq("../build_linux/vq_stage1.f32", target_fn, frames, dec, EbNodB(i), ldpc_en=0, verbose=0);
+    results1_log = [results1_log results];
   end
   results4_log = [];
   for i=1:length(EbNodB)
-    results = run_test_vq("../build_linux/vq_stage1_bs004.f32", target_fn, frames, dec, EbNodB(i), verbose=0);
+    results = run_test_vq("../build_linux/vq_stage1_bs004.f32", target_fn, frames, dec, EbNodB(i), ldpc_en=0, verbose=0);
     results4_log = [results4_log results];
   end
   for i=1:length(results1_log)
     ber(i) = results1_log(i).ber;
     per(i) = results1_log(i).per;
     mse_noerrors(i) = sqrt(results1_log(i).mse_noerrors);
     mse_vq1(i) = sqrt(results1_log(i).mse);
-    mse_vq2(i) = sqrt(results2_log(i).mse);
     mse_vq4(i) = sqrt(results4_log(i).mse);
+    mse_vq1_ldpc(i) = sqrt(results1_ldpc_log(i).mse);
+    mse_vq4_ldpc(i) = sqrt(results4_ldpc_log(i).mse);
   end
 
   figure(1); clf;
   semilogy(EbNodB, ber, 'g+-;ber;','linewidth', 2); hold on;
   semilogy(EbNodB, per, 'b+-;per;','linewidth', 2);
   grid('minor'); xlabel('Eb/No(dB)');
-
-hold off;
+  hold off;
 
   figure(2); clf;
   plot(EbNodB, mse_noerrors, "b+-;no errors;"); hold on;
-  plot(EbNodB, mse_vq1, "g+-;vanilla;");
-  plot(EbNodB, mse_vq2, "r+-;binary switch;");
-  plot(EbNodB, mse_vq4, "b+-;binary switch with prob;");
+  plot(EbNodB, mse_vq1, "g+-;vanilla AWGN;");
+  plot(EbNodB, mse_vq4, "b+-;binary switch;");
+  plot(EbNodB, mse_vq1_ldpc, "r+-;ldpc (112,56);");
+  plot(EbNodB, mse_vq4_ldpc, "k+-;binary switch ldpc (112,56);");
+  load trellis_dec3_nstage3.txt
+  plot(EbNodB, rms_sd, "c+-;binary switch trellis dec3;");  
   hold off; grid; title("RMS SD (dB)"); xlabel('Eb/No(dB)');
 endfunction