moved hf channel sim into a function
authordrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Wed, 15 Oct 2014 01:32:04 +0000 (01:32 +0000)
committerdrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Wed, 15 Oct 2014 01:32:04 +0000 (01:32 +0000)
git-svn-id: https://svn.code.sf.net/p/freetel/code@1889 01035d8c-6547-0410-b346-abe4f91aad63

codec2-dev/octave/test_ldpc.m

index 392039d78742d876f32261597047d65500a1a52a..fcdd4cad868abe737482b838d651717489706d9c 100644 (file)
@@ -21,11 +21,124 @@ randn('state',1);
 
 1;
 
+% Symbol rate processing for tx side (modulator)
+
+function [tx_symb tx_bits prev_sym_tx] = tx_symbol_rate(sim_in, tx_bits, code_param, prev_sym_tx)
+    ldpc_code     = sim_in.ldpc_code;
+    framesize     = sim_in.framesize;
+    rate          = sim_in.rate;
+    Nsymbrow      = sim_in.Nsymbrow;
+    Nsymbrowpilot = sim_in.Nsymbrowpilot;
+    Nc            = sim_in.Nc;
+    Npilotsframe  = sim_in.Npilotsframe;
+    Ns            = sim_in.Ns;
+    Nchip         = sim_in.Nchip;
+    modulation    = sim_in.modulation;
+    pilot         = sim_in.pilot;
+
+    if ldpc_code
+        [tx_bits, tmp] = ldpc_enc(tx_bits, code_param);
+    end
+
+    % modulate --------------------------------------------
+
+    % organise symbols into a Nsymbrow rows by Nc cols
+    % data and parity bits are on separate carriers
+
+    tx_symb = zeros(Nsymbrow,Nc);
+
+    for c=1:Nc
+      for r=1:Nsymbrow
+        i = (c-1)*Nsymbrow + r;
+        tx_symb(r,c) = qpsk_mod(tx_bits(2*(i-1)+1:2*i));
+      end
+    end
+
+    % Optionally insert pilots, one every Ns data symbols
+
+    tx_symb_pilot = zeros(Nsymbrowpilot, Nc);
+            
+    for p=1:Npilotsframe
+      tx_symb_pilot((p-1)*(Ns+1)+1,:)          = pilot(p,:);                 % row of pilots
+      %printf("%d %d %d %d\n", (p-1)*(Ns+1)+2, p*(Ns+1), (p-1)*Ns+1, p*Ns);
+      tx_symb_pilot((p-1)*(Ns+1)+2:p*(Ns+1),:) = tx_symb((p-1)*Ns+1:p*Ns,:); % payload symbols
+    end
+    tx_symb = tx_symb_pilot;
+
+    % Optionally copy to other carriers (spreading)
+
+    for c=Nc+1:Nc:Nc*Nchip
+      tx_symb(:,c:c+Nc-1) = tx_symb(:,1:Nc);
+    end
+            
+    % Optionally DQPSK encode
+    if strcmp(modulation,'dqpsk')
+      for c=1:Nc*Nchip
+        for r=1:Nsymbrowpilot
+          tx_symb(r,c) *= prev_sym_tx(c);
+          prev_sym_tx(c) = tx_symb(r,c);
+        end
+      end               
+    end
+
+    % ensures energy/symbol is normalised when spreading
+
+    tx_symb = tx_symb/sqrt(Nchip);
+end
+
+
+% Init HF channel model from stored sample files of spreading signal ----------------------------------
+
+function [spread spread_2ms hf_gain] = init_hf_model(Fs, Rs, nsam)
+
+    % convert "spreading" samples from 1kHz carrier at Fs to complex
+    % baseband, generated by passing a 1kHz sine wave through PathSim
+    % with the ccir-poor model, enabling one path at a time.
+    
+    Fc = 1000; M = Fs/Rs;
+    fspread = fopen("../raw/sine1k_2Hz_spread.raw","rb");
+    spread1k = fread(fspread, "int16")/10000;
+    fclose(fspread);
+    fspread = fopen("../raw/sine1k_2ms_delay_2Hz_spread.raw","rb");
+    spread1k_2ms = fread(fspread, "int16")/10000;
+    fclose(fspread);
+
+    % down convert to complex baseband
+    spreadbb = spread1k.*exp(-j*(2*pi*Fc/Fs)*(1:length(spread1k))');
+    spreadbb_2ms = spread1k_2ms.*exp(-j*(2*pi*Fc/Fs)*(1:length(spread1k_2ms))');
+
+    % remove -2000 Hz image
+    b = fir1(50, 5/Fs);
+    spread = filter(b,1,spreadbb);
+    spread_2ms = filter(b,1,spreadbb_2ms);
+   
+    % discard first 1000 samples as these were near 0, probably as
+    % PathSim states were ramping up
+
+    spread    = spread(1000:length(spread));
+    spread_2ms = spread_2ms(1000:length(spread_2ms));
+
+    % decimate down to Rs
+    
+    spread = spread(1:M:length(spread));
+    spread_2ms = spread_2ms(1:M:length(spread_2ms));
+
+    % Determine "gain" of HF channel model, so we can normalise
+    % carrier power during HF channel sim to calibrate SNR.  I imagine
+    % different implementations of ccir-poor would do this in
+    % different ways, leading to different BER results.  Oh Well!
+
+    hf_gain = 1.0/sqrt(var(spread(1:nsam))+var(spread_2ms(1:nsam)));
+endfunction
+
+
 % main test function 
 
-function sim_out = ber_test(sim_in, modulation)
+function sim_out = ber_test(sim_in)
     Fs = 8000;
 
+    modulation       = sim_in.modulation;
     verbose          = sim_in.verbose;
     framesize        = sim_in.framesize;
     Ntrials          = sim_in.Ntrials;
@@ -45,12 +158,14 @@ function sim_out = ber_test(sim_in, modulation)
     Np               = sim_in.Np;     % number of pilots to use
     Ns               = sim_in.Ns;     % step size between pilots
     ldpc_code        = sim_in.ldpc_code;
+    sim_in.rate = rate = sim_in.ldpc_code_rate; 
 
     bps              = 2;
-    Nsymb            = framesize/bps;
-    Nsymbrow         = Nsymb/Nc;
-    Npilotsframe     = Nsymbrow/Ns;
-    Nsymbrowpilot    = Nsymbrow + Npilotsframe;
+
+    sim_in.Nsymb         = Nsymb            = framesize/bps;
+    sim_in.Nsymbrow      = Nsymbrow         = Nsymb/Nc;
+    sim_in.Npilotsframe  = Npilotsframe     = Nsymbrow/Ns;
+    sim_in.Nsymbrowpilot = Nsymbrowpilot    = Nsymbrow + Npilotsframe;
 
     printf("Each frame is %d bits or %d symbols, transmitted as %d symbols by %d carriers.",
            framesize, Nsymb, Nsymbrow, Nc);
@@ -74,6 +189,7 @@ function sim_out = ber_test(sim_in, modulation)
     for c=1:Nc
       pilot(:,c) = [ones(1,floor(Npilotsframe/2)) -ones(1,ceil(Npilotsframe/2))]';
     end
+    sim_in.pilot = pilot;
     tx_pilot_buf = [pilot; pilot; pilot];
    
     % Init LDPC --------------------------------------------------------------------
@@ -91,7 +207,6 @@ function sim_out = ber_test(sim_in, modulation)
 
         ldpc;
 
-        rate = sim_in.ldpc_code_rate; 
         mod_order = 4; 
         modulation2 = 'QPSK';
         mapping = 'gray';
@@ -109,46 +224,7 @@ function sim_out = ber_test(sim_in, modulation)
         rate = 1;
     end
 
-    % Init HF channel model from stored sample files of spreading signal ----------------------------------
-
-    % convert "spreading" samples from 1kHz carrier at Fs to complex
-    % baseband, generated by passing a 1kHz sine wave through PathSim
-    % with the ccir-poor model, enabling one path at a time.
-    
-    Fc = 1000; M = Fs/Rs;
-    fspread = fopen("../raw/sine1k_2Hz_spread.raw","rb");
-    spread1k = fread(fspread, "int16")/10000;
-    fclose(fspread);
-    fspread = fopen("../raw/sine1k_2ms_delay_2Hz_spread.raw","rb");
-    spread1k_2ms = fread(fspread, "int16")/10000;
-    fclose(fspread);
-
-    % down convert to complex baseband
-    spreadbb = spread1k.*exp(-j*(2*pi*Fc/Fs)*(1:length(spread1k))');
-    spreadbb_2ms = spread1k_2ms.*exp(-j*(2*pi*Fc/Fs)*(1:length(spread1k_2ms))');
-
-    % remove -2000 Hz image
-    b = fir1(50, 5/Fs);
-    spread = filter(b,1,spreadbb);
-    spread_2ms = filter(b,1,spreadbb_2ms);
-   
-    % discard first 1000 samples as these were near 0, probably as
-    % PathSim states were ramping up
-
-    spread    = spread(1000:length(spread));
-    spread_2ms = spread_2ms(1000:length(spread_2ms));
-
-    % decimate down to Rs
-    
-    spread = spread(1:M:length(spread));
-    spread_2ms = spread_2ms(1:M:length(spread_2ms));
-
-    % Determine "gain" of HF channel model, so we can normalise
-    % carrier power during HF channel sim to calibrate SNR.  I imagine
-    % different implementations of ccir-poor would do this in
-    % different ways, leading to different BER results.  Oh Well!
-
-    hf_gain = 1.0/sqrt(var(spread(1:Nsymbrow*Ntrials))+var(spread_2ms(1:Nsymbrow*Ntrials)));
+    [spread spread_2ms hf_gain] = init_hf_model(Fs, Rs, Nsymbrowpilot*Ntrials);
 
     % Start Simulation ----------------------------------------------------------------
 
@@ -157,7 +233,7 @@ function sim_out = ber_test(sim_in, modulation)
         EsNo = 10^(EsNodB/10);
     
         variance = 1/EsNo;
-         if verbose > 1
+        if verbose > 1
             printf("EsNo (dB): %f EsNo: %f variance: %f\n", EsNodB, EsNo, variance);
         end
         
@@ -185,60 +261,16 @@ function sim_out = ber_test(sim_in, modulation)
         for nn = 1:Ntrials+2
                   
             if ldpc_code
-                tx_bits = round(rand(1,framesize*rate));                       
-                [tx_bits, tmp] = ldpc_enc(tx_bits, code_param);
+              tx_bits = round(rand(1,framesize*rate));                       
             else
-                tx_bits = round(rand(1,framesize));                       
+              tx_bits = round(rand(1,framesize));                       
             end
+
+            [s_ch tx_bits prev_sym_tx] = tx_symbol_rate(sim_in, tx_bits, code_param, prev_sym_tx);
+   
             tx_bits_buf(1:framesize) = tx_bits_buf(framesize+1:2*framesize);
             tx_bits_buf(framesize+1:2*framesize) = tx_bits;
 
-            % modulate --------------------------------------------
-
-            % organise symbols into a Nsymbrow rows by Nc cols
-            % data and parity bits are on separate carriers
-
-            tx_symb = zeros(Nsymbrow,Nc);
-
-            for c=1:Nc
-              for r=1:Nsymbrow
-                i = (c-1)*Nsymbrow + r;
-                tx_symb(r,c) = qpsk_mod(tx_bits(2*(i-1)+1:2*i));
-              end
-            end
-
-            % Optionally insert pilots, one every Ns data symbols
-
-            tx_symb_pilot = zeros(Nsymbrowpilot, Nc);
-            
-            for p=1:Npilotsframe
-              tx_symb_pilot((p-1)*(Ns+1)+1,:)          = pilot(p,:);                 % row of pilots
-              %printf("%d %d %d %d\n", (p-1)*(Ns+1)+2, p*(Ns+1), (p-1)*Ns+1, p*Ns);
-              tx_symb_pilot((p-1)*(Ns+1)+2:p*(Ns+1),:) = tx_symb((p-1)*Ns+1:p*Ns,:); % payload symbols
-            end
-            tx_symb = tx_symb_pilot;
-
-            % Optionally copy to other carriers (spreading)
-
-            for c=Nc+1:Nc:Nc*Nchip
-              tx_symb(:,c:c+Nc-1) = tx_symb(:,1:Nc);
-            end
-            
-            % Optionally DQPSK encode
-            if strcmp(modulation,'dqpsk')
-              for c=1:Nc*Nchip
-                for r=1:Nsymbrowpilot
-                  tx_symb(r,c) *= prev_sym_tx(c);
-                  prev_sym_tx(c) = tx_symb(r,c);
-                end
-              end               
-            end
-
-            % ensures energy/symbol is normalised when spreading
-
-            s_ch = tx_symb/sqrt(Nchip);
-
             % HF channel simulation  ------------------------------------
             
             hf_fading = ones(1,Nsymb);
@@ -628,12 +660,13 @@ function test_single
   sim_in.ldpc_code_rate   = 0.5;
   sim_in.ldpc_code        = 1;
 
-  sim_in.Ntrials          = 20;
+  sim_in.Ntrials          = 10;
   sim_in.Esvec            = 7; 
   sim_in.hf_sim           = 1;
   sim_in.hf_mag_only      = 0;
-  
-  sim_qpsk_hf             = ber_test(sim_in, 'qpsk');
+  sim_in.modulation       = 'qpsk';
+
+  sim_qpsk_hf             = ber_test(sim_in);
 
   fep=fopen("errors_450.bin","wb"); fwrite(fep, sim_qpsk_hf.ldpc_errors_log, "short"); fclose(fep);
 endfunction