variable power simulation for energy modulation works quite well, some quality improv...
authordrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Sun, 20 Apr 2014 08:28:35 +0000 (08:28 +0000)
committerdrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Sun, 20 Apr 2014 08:28:35 +0000 (08:28 +0000)
git-svn-id: https://svn.code.sf.net/p/freetel/code@1513 01035d8c-6547-0410-b346-abe4f91aad63

codec2-dev/octave/fuzzy_gray.m

index abb9e12e6f7314efb11cce596dcf1cc5f646d4cd..7fcca72cb383d57f20a77d9c600d9e3c696acd8d 100644 (file)
@@ -2,8 +2,9 @@
 % David Rowe
 % 10 April 2014
 %
-% Experiments in fuzzy gray codes.  Idea is that with one bit error
-% in the codeword only changes the encoded value by at most 1.
+% Various experiments in fuzzy gray codes and quantising and
+% transmitting scalars.  Idea is that with one bit error in the
+% codeword only changes the encoded value by at most 1.
 
 1;
 
@@ -69,11 +70,12 @@ function natural = gray_to_binary(gray)
     end
 endfunction
 
-function sim_out = test_baseline_uncoded(Ebvec, Nbits, Ntrials)
+function sim_out = test_baseline_uncoded(Ebvec, Nbits, Ntrials, enable_error_log)
     Nlevels = 2.^ Nbits; powersOfTwo = 2 .^ fliplr(0:(Nbits-1));
     Nsymb   = Nbits;
 
     sim_out.qnoise_log = zeros(length(Ebvec),Ntrials);
+    sim_out.error_log  = [];
 
     for ne = 1:length(Ebvec)
         EbNodB = Ebvec(ne);
@@ -103,6 +105,65 @@ function sim_out = test_baseline_uncoded(Ebvec, Nbits, Ntrials)
             Nerrs = sum(error_positions);
             Terrs += Nerrs;
             Tbits += length(tx_bits);
+            if enable_error_log
+                sim_out.error_log  = [sim_out.error_log error_positions];
+            end
+
+            rx_index = (powersOfTwo  * rx_bits');
+            rx_value = unquantise_value(rx_index, 0, 1, Nlevels);
+
+            qsignal += tx_value*tx_value;
+            qnoise  += (tx_value - rx_value) .^ 2;
+            sim_out.qnoise_log(ne,nn) = tx_value - rx_value;
+        end
+
+        sim_out.TERvec(ne) = Terrs;
+        sim_out.BERvec(ne) = Terrs/Tbits;
+        sim_out.QSNRvec(ne) = 10*log10(qsignal/qnoise);
+        printf("EbNo (dB): %3.2f  Terrs: %6d BER %1.4f QSNR (dB): %3.2f\n", 
+        EbNodB, Terrs, Terrs/Tbits, 10*log10(qsignal/qnoise));
+    end
+
+endfunction
+
+function sim_out = test_varpower(Ebvec, Nbits, Ntrials, amps, enable_error_log)
+    Nlevels = 2.^ Nbits; powersOfTwo = 2 .^ fliplr(0:(Nbits-1));
+    Nsymb   = Nbits;
+
+    sim_out.qnoise_log = zeros(length(Ebvec), Ntrials);
+    sim_out.error_log  = [];
+
+    for ne = 1:length(Ebvec)
+        EbNodB = Ebvec(ne);
+        EbNo = 10^(EbNodB/10);
+    
+        variance = 1/EbNo;
+        
+        Terrs = 0;  Tbits = 0;
+        qsignal = qnoise = 0;
+        for nn = 1:Ntrials
+                  
+            tx_value = rand(1,1);
+            tx_index = quantise_value(tx_value, 0, 1, Nlevels);
+            tx_bits = dec2bin(tx_index, Nbits) - '0';
+            tx_symbols = (-1 + 2*tx_bits) .* amps; 
+
+            % AWGN noise and phase/freq offset channel simulation
+            % 0.5 factor ensures var(noise) == variance , i.e. splits power between Re & Im
+
+            noise = sqrt(variance*0.5)*(randn(1,Nsymb) + j*randn(1,Nsymb));
+            rx_symbols = tx_symbols + noise;
+
+            rx_bits = rx_symbols > 0;
+
+            error_positions = xor(rx_bits, tx_bits);
+            if enable_error_log
+                sim_out.error_log  = [sim_out.error_log error_positions];
+            end
+            Nerrs = sum(error_positions);
+            Terrs += Nerrs;
+            Tbits += length(tx_bits);
 
             rx_index = (powersOfTwo  * rx_bits');
             rx_value = unquantise_value(rx_index, 0, 1, Nlevels);
@@ -132,7 +193,7 @@ function valid_codewords = fuzzy_code_create(ndata,nparity)
     bad_distance = 0;
     for i=1:Nvalid
         for k=i+1:Nvalid
-            distance = sum(bitxor(valid_codewords(i,:), valid_codewords(k,:)));
+            distance = sum(bitxor(valid_codewords(i,:), valid_codewords(k,:)))
             if distance < 2
                 bad_distance++;
             end
@@ -162,33 +223,47 @@ function value = fuzzy_code_decode(codewords, rx_symbols)
     end
 endfunction
 
-function valid_codewords = fuzzy_code_create(ndata, nparity)
+function set_codewords = fuzzy_code_create(ndata, nparity)
     Nbits = ndata + nparity;
     Nvalid = 2 .^ ndata;
-    codewords = binary_to_gray(0:(2 .^ Nbits)-1);
-    valid_codewords = dec2bin(codewords(1:2:(2 .^ Nbits)), Nbits) - '0';
+    M = 2 .^ Nbits;
 
-    % check all valid codewords have a hamming distance of at least 2^nparity    
+    codewords = dec2bin(0:(M-1), Nbits) - '0';
 
-    bad_distance = 0;
-    for i=1:Nvalid
-        for k=i+1:Nvalid
-            distance = sum(bitxor(valid_codewords(i,:), valid_codewords(k,:)));
-            if distance < 2
-                bad_distance++;
+    set_codewords(1,:) = codewords(1,:);
+
+    for i=1:Nvalid-1
+
+        % measure distance from current set to all codewords
+
+        for m=1:i
+            printf("%d: ", m);
+            for k=1:length(codewords);
+                d(m,k) = sum(bitxor(set_codewords(m,:), codewords(k,:)));
+                printf("%d ", d(m,k));
             end
+            printf("\n");
         end
-    end
-    if bad_distance != 0
-       printf("Error: Nvalid: %d bad_distance: %d\n", Nvalid, bad_distance);
-       return;
-    end
 
+        % choose new codeword for set that has maximum distance from current set
+
+        dmax = 0; best_cw = 1;
+        for k=1:length(codewords);
+            dmin = min(d(:,k));
+            if dmin > dmax
+                dmax = dmin;
+                best_cw = k;
+            end
+        end
+
+        printf("[%d %d]\n", dmax, best_cw);
+        set_codewords(i+1,:) = codewords(best_cw,:)
+    end
+    
 endfunction
 
 
-function sim_out = test_fuzzy_code(Ebvec, Ndata, Ntrials)
-    Nparity = 1;
+function sim_out = test_fuzzy_code(Ebvec, Ndata, Nparity, Ntrials)
     Nbits   = Ndata + Nparity;
     Nlevels = 2 .^ Ndata;
     Nsymb   = Nbits;
@@ -232,38 +307,129 @@ function sim_out = test_fuzzy_code(Ebvec, Ndata, Ntrials)
     end
 endfunction
 
-Ebvec   = 0:7;
-Ntrials = 1000;
-Nbits   = 5;
-
-baseline = test_baseline_uncoded(Ebvec, Nbits, Ntrials);
-fuzzy = test_fuzzy_code(Ebvec, Nbits, Ntrials);
-
-figure(1);
-clf;
-semilogy(Ebvec, baseline.BERvec)
-xlabel('Eb/N0')
-ylabel('BER')
-grid("minor")
-
-figure(2);
-clf;
-plot(Ebvec, baseline.QSNRvec,'b;baseline;')
-hold on;
-plot(Ebvec, fuzzy.QSNRvec,'r;fuzzy;')
-hold off;
-xlabel('Eb/N0')
-ylabel('SNR')
-grid("minor")
-
-figure(3);
-subplot(211)
-hist(baseline.qnoise_log(4,:),50);
-subplot(212)
-hist(fuzzy.qnoise_log(4,:),50);
-
-figure(4)
-plot(baseline.qnoise_log(4,1:250),'b;baseline;')
-hold on;
-plot(fuzzy.qnoise_log(4,1:250),'r;fuzzy;')
-hold off;
+function compare_baseline_fuzzy
+    Ebvec   = 0:7;
+    Ntrials = 1000;
+    Nbits   = 3; Nparity = 1;
+
+    baseline = test_baseline_uncoded(Ebvec, Nbits, Ntrials);
+    fuzzy = test_fuzzy_code(Ebvec, Nbits, Nparity, Ntrials);
+
+    figure(1);
+    clf;
+    semilogy(Ebvec, baseline.BERvec)
+    xlabel('Eb/N0')
+    ylabel('BER')
+    grid("minor")
+
+    figure(2);
+    clf;
+    plot(Ebvec, baseline.QSNRvec,'b;baseline;')
+    hold on;
+    plot(Ebvec, fuzzy.QSNRvec,'r;fuzzy;')
+    hold off;
+    xlabel('Eb/N0')
+    ylabel('SNR')
+    grid("minor")
+
+    figure(3);
+    subplot(211)
+    hist(baseline.qnoise_log(4,:),50);
+    subplot(212)
+    hist(fuzzy.qnoise_log(4,:),50);
+
+    figure(4)
+    plot(baseline.qnoise_log(1,1:250),'b;baseline;')
+    hold on;
+    plot(fuzzy.qnoise_log(1,1:250),'r;fuzzy;')
+    hold off;
+endfunction
+
+function compare_baseline_varpower_plot
+    Ebvec   = -2:5;
+    Ntrials = 5000;
+    Nbits   = 5;
+
+    baseline = test_baseline_uncoded(Ebvec, Nbits, Ntrials, 0);
+    amps = [2 1.5 1.0 0.5 0.5];
+    av_pwr = (amps*amps')/length(amps);
+    amps_norm = amps/sqrt(av_pwr)
+    varpower = test_varpower(Ebvec, Nbits, Ntrials,  amps_norm, 0);
+
+    figure(1);
+    clf;
+    semilogy(Ebvec, baseline.BERvec)
+    xlabel('Eb/No (dB)')
+    ylabel('BER')
+    grid("minor")
+    title('BER versus Eb/No')
+
+    figure(2);
+    clf;
+    plot(Ebvec, baseline.QSNRvec,'b;baseline;')
+    hold on;
+    plot(Ebvec, varpower.QSNRvec,'r;varpower;')
+    hold off;
+    xlabel('Eb/No (dB)')
+    ylabel('SNR (dB)')
+    grid("minor")
+    title('Quantiser SNR versus Eb/No')
+
+    figure(3);
+    subplot(211)
+    hist(baseline.qnoise_log(1,:),50);
+    title('Baseline and Variable Power Error Histograms')
+    subplot(212)
+    hist(varpower.qnoise_log(1,:),50);
+
+    figure(4)
+    subplot(211)
+    plot(baseline.qnoise_log(1,1:250),'b;baseline;')
+    title('Baseline and Variable Power Error plots for Eb/No = -2dB')
+    subplot(212)
+    plot(varpower.qnoise_log(1,1:250),'r;varpower;')
+endfunction
+
+function compare_baseline_varpower_error_files
+    Ebvec   = -2;
+    Fs      = 25;         % number of samples per second
+    Nsec    = 15;         % seconds to simulate
+    Ntrials = Fs*Nsec;
+    Nbits   = 5;
+    bits_per_frame = 52;
+    bits_per_frame_rounded = ceil(bits_per_frame/8)*8; % c2enc uses integer number of bytes/frame
+    start_bit = 12;                                    % first energy bit (after 4 voicing, 7 Wo bits)
+
+    baseline = test_baseline_uncoded(Ebvec, Nbits, Ntrials, 1);
+    amps = [2 1.5 1.0 0.5 0.5];
+    av_pwr = (amps*amps')/length(amps);
+    amps_norm = amps/sqrt(av_pwr);
+    varpower = test_varpower(Ebvec, Nbits, Ntrials, amps_norm, 1);
+
+    % construct error patterns to apply to c2enc bit stream
+
+    baseline_errors = [];
+    for i=1:Ntrials
+        error_positions = baseline.error_log(Nbits*(i-1)+1:Nbits*i);
+        baseline_errors = [baseline_errors zeros(1,start_bit-1) error_positions ...
+                           zeros(1, bits_per_frame_rounded - Nbits - (start_bit-1))];
+    end
+    varpower_errors = [];
+    for i=1:Ntrials
+        error_positions = varpower.error_log(Nbits*(i-1)+1:Nbits*i);
+        varpower_errors = [varpower_errors zeros(1,start_bit-1) error_positions ...
+                           zeros(1, bits_per_frame_rounded - Nbits - (start_bit-1))];
+    end
+    % save error patterns
+
+    fep=fopen("energy_errors_baseline.bin","wb"); fwrite(fep, baseline_errors, "short"); fclose(fep);
+    fep=fopen("energy_errors_varpower.bin","wb"); fwrite(fep, varpower_errors, "short"); fclose(fep);
+endfunction
+
+%compare_baseline_varpower_plot
+%compare_baseline_varpower_error_files
+
+%compare_baseline_fuzzy
+%fuzzy_code_create(3,1)