fading model working OK for different symbol and sampling rates

git-svn-id: https://svn.code.sf.net/p/freetel/code@2144 01035d8c-6547-0410-b346-abe4f91aad63

diff --git a/codec2-dev/octave/cohpsk.m b/codec2-dev/octave/cohpsk.m
index b6a83148c53ff8cd1dab066b5fbda33f52e4888d..868fa0d93822cfde32f3f1d2c4db5fd0027f8620 100644 (file)
--- a/codec2-dev/octave/cohpsk.m
+++ b/codec2-dev/octave/cohpsk.m
@@ -245,7 +245,7 @@ function [rx_symb rx_bits rx_symb_linear amp_ phi_ EsNo_ cohpsk] = qpsk_symbols_
     % better on fading channels, but increases BER slighlty for AWGN
     % channels.
 
-    sampling_points = [1 2 7 8];
+    sampling_points = [1 2 cohpsk.Nsymbrow+3 cohpsk.Nsymbrow+4];
     pilot2 = [cohpsk.pilot(1,:); cohpsk.pilot(2,:); cohpsk.pilot(1,:); cohpsk.pilot(2,:);];
     phi_ = zeros(Nsymbrow, Nc*Nd);
     amp_ = zeros(Nsymbrow, Nc*Nd);
@@ -369,13 +369,13 @@ endfunction
 
 % Init HF channel model from stored sample files of spreading signal ----------------------------------
 
-function [spread spread_2ms hf_gain] = init_hf_model(Fs, Rs, nsam)
+function [spread spread_2ms hf_gain] = init_hf_model(Fs, nsam)
 
-    % convert "spreading" samples from 1kHz carrier at Fs to complex
+    % convert "spreading" samples from 1kHz carrier at Fss 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;
+    Fc = 1000; Fss = 8000;
     fspread = fopen("../raw/sine1k_2Hz_spread.raw","rb");
     spread1k = fread(fspread, "int16")/10000;
     fclose(fspread);
@@ -384,11 +384,11 @@ function [spread spread_2ms hf_gain] = init_hf_model(Fs, Rs, nsam)
     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))');
+    spreadbb = spread1k.*exp(-j*(2*pi*Fc/Fss)*(1:length(spread1k))');
+    spreadbb_2ms = spread1k_2ms.*exp(-j*(2*pi*Fc/Fss)*(1:length(spread1k_2ms))');
 
     % remove -2000 Hz image
-    b = fir1(50, 5/Fs);
+    b = fir1(50, 5/Fss);
     spread = filter(b,1,spreadbb);
     spread_2ms = filter(b,1,spreadbb_2ms);
    
@@ -398,10 +398,10 @@ function [spread spread_2ms hf_gain] = init_hf_model(Fs, Rs, nsam)
     spread    = spread(1000:length(spread));
     spread_2ms = spread_2ms(1000:length(spread_2ms));
 
-    % decimate down to Rs
+    % change output samples so they are at rate Fs reqd by caller
     
-    spread = spread(1:M:length(spread));
-    spread_2ms = spread_2ms(1:M:length(spread_2ms));
+    spread = resample(spread, Fs, Fss);
+    spread_2ms = resample(spread_2ms, Fs, Fss);
 
     % Determine "gain" of HF channel model, so we can normalise
     % carrier power during HF channel sim to calibrate SNR.  I imagine
@@ -682,7 +682,7 @@ function [next_sync cohpsk] = frame_sync_fine_freq_est(cohpsk, ch_symb, sync, ne
  
   % sample pilots at start of this frame and start of next frame 
 
-  sampling_points = [1 2 7 8];
+  sampling_points = [1 2 cohpsk.Nsymbrow+3 cohpsk.Nsymbrow+4];
   pilot2 = [ cohpsk.pilot(1,:); cohpsk.pilot(2,:); cohpsk.pilot(1,:); cohpsk.pilot(2,:);];
 
   if sync == 0
@@ -812,7 +812,7 @@ function sim_out = ber_test(sim_in)
     f_off            = sim_in.f_off;
     div_time_shift   = sim_in.div_timeshift;
 
-    [spread spread_2ms hf_gain] = init_hf_model(Fs, Rs, Nsymbrowpilot*(Ntrials+2));
+    [spread spread_2ms hf_gain] = init_hf_model(Rs, Nsymbrowpilot*(Ntrials+2));
 
     if strcmp(modulation,'dqpsk')
       Nsymbrowpilot = Nsymbrow;

diff --git a/codec2-dev/octave/fdmdv.m b/codec2-dev/octave/fdmdv.m
index b872697b6733a234030ada5333776785c8820b91..68c137bb9619986c23a49a8fcd7800fd497ae657 100644 (file)
--- a/codec2-dev/octave/fdmdv.m
+++ b/codec2-dev/octave/fdmdv.m
@@ -890,7 +890,8 @@ function [pilot prev_pilot pilot_lut_index prev_pilot_lut_index] = get_pilot(pil
   end
 endfunction
 
-
+if 0
+% want to use Octave resample function!
 
 % Change the sample rate by a small amount, for example 1000ppm (ratio
 % = 1.001).  Always returns nout samples in buf_out, but uses a
@@ -928,7 +929,7 @@ function [buf_out t nin] = resample(buf_in, t, ratio, nout)
   t -= delta;
 
 endfunction
-
+end
 
 % freq offset state machine.  Moves between acquire and track states based
 % on BPSK pilot sequence.  Freq offset estimator occasionally makes mistakes

diff --git a/codec2-dev/octave/tcohpsk.m b/codec2-dev/octave/tcohpsk.m
index ffeebef37ccac9dc9efffdae762e52a5768463c7..3abc785b5008068c9167def97a89980cbc9489d7 100644 (file)
--- a/codec2-dev/octave/tcohpsk.m
+++ b/codec2-dev/octave/tcohpsk.m
@@ -46,8 +46,8 @@ randn('state',1);
 % select which test  ----------------------------------------------------------
 
 %test = 'compare to c';
-test = 'awgn';
-%test = 'fading';
+%test = 'awgn';
+test = 'fading';
 
 % some parameters that can be over ridden, e.g. to disable parts of modem
 
@@ -69,10 +69,10 @@ end
 % should be BER around 0.015 to 0.02
 
 if strcmp(test, 'awgn')
-  frames = 10;
+  frames = 100;
   foff =  0;
   dfoff = 0;
-  EsNodB = 80;
+  EsNodB = 8;
   fading_en = 0;
   hf_delay_ms = 2;
   compare_with_c = 0;
@@ -82,7 +82,7 @@ end
 
 if strcmp(test, 'fading');
   frames = 100;
-  foff = -52.3;
+  foff = -53.1;
   dfoff = 0.5/Fs;
   EsNodB = 12;
   fading_en = 1;
@@ -95,19 +95,29 @@ EsNo = 10^(EsNodB/10);
 % modem constants ----------------------------------------------------------
 
 Rs = 75;               % symbol rate in Hz
-Nc = 4;                % number of carriers
+Nc = 7;                % number of carriers
 Nd = 2;                % diveristy factor
-framesize = 32;        % number of payoad data bits in the frame
+framesize = 56;        % number of payload data bits in the frame
 
 Nsw = 3;               % frames we demod for initial sync window
 afdmdv.Nsym = 6;       % size of tx/tx root nyquist filter in symbols
 afdmdv.Nt = 5;         % number of symbols we estimate timing over
 
+clip = 6.5;            % Clipping of tx signal to reduce PAPR. Adjust by 
+                       % experiment as Nc and Nd change.  Check out no noise 
+                       % scatter diagram and AWGN/fading BER perf
+                       % at operating points
+
 % FDMDV init ---------------------------------------------------------------
 
 Fs = afdmdv.Fs = 7500;
 afdmdv.Nc = Nd*Nc-1;
 afdmdv.Rs = Rs;
+if Fs/afdmdv.Rs != floor(Fs/afdmdv.Rs)
+  printf("\n  Oops, Fs/Rs must be an integer!\n\n");
+  return
+end
+
 M = afdmdv.M  = afdmdv.Fs/afdmdv.Rs;
 afdmdv.Nfilter = afdmdv.Nsym*M;
 afdmdv.tx_filter_memory = zeros(afdmdv.Nc+1, afdmdv.Nfilter);
@@ -135,6 +145,9 @@ afdmdv.phase_rx = ones(afdmdv.Nc+1,1);
 afdmdv.Nfilter = afdmdv.Nsym*afdmdv.M;
 afdmdv.rx_fdm_mem = zeros(1,afdmdv.Nfilter + afdmdv.M);
 Q = afdmdv.Q = afdmdv.M/4;
+if Q != floor(Q)
+  printf("\n  Yeah .... if (Fs/Rs)/4 = M/4 isn't an integer we will just go and break things.\n\n");
+end
 
 afdmdv.rx_filter_mem_timing = zeros(afdmdv.Nc+1, afdmdv.Nt*afdmdv.P);
 afdmdv.Nfiltertiming = afdmdv.M + afdmdv.Nfilter + afdmdv.M;
@@ -196,7 +209,7 @@ tx_bits_coh = round(rand(1,framesize*10));
 ptx_bits_coh = 1;
 
 Nerrs = Tbits = 0;
-prev_tx_bits = [];
+prev_tx_bits = prev_tx_bits2 = [];
 
 phase_ch = 1;
 sync = initial_sync;
@@ -205,7 +218,7 @@ acohpsk.f_fine_est = 0;
 acohpsk.ct = 4;
 acohpsk.ftrack_en = ftrack_en;
 
-[spread spread_2ms hf_gain] = init_hf_model(Fs, Fs, frames*acohpsk.Nsymbrowpilot*afdmdv.M);
+[spread spread_2ms hf_gain] = init_hf_model(Fs, frames*acohpsk.Nsymbrowpilot*afdmdv.M);
 hf_n = 1;
 nhfdelay = floor(hf_delay_ms*Fs/1000);
 ch_fdm_delay = zeros(1, acohpsk.Nsymbrowpilot*M + nhfdelay);
@@ -234,7 +247,6 @@ for f=1:frames
     tx_fdm_frame = [tx_fdm_frame tx_fdm];
   end
 
-  clip = 5;
   ind = find(abs(tx_fdm_frame) > clip);
   tx_fdm_frame(ind) = clip*exp(j*angle(tx_fdm_frame(ind)));
 
@@ -488,7 +500,7 @@ else
   
   papr = max(tx_fdm_frame_log.*conj(tx_fdm_frame_log)) / mean(tx_fdm_frame_log.*conj(tx_fdm_frame_log));
   papr_dB = 10*log10(papr);
-  printf("av tx pwr: %f PAPR: %4.2f dB av rx pwr: %f\n", var(tx_fdm_frame_log), papr_dB, var(ch_fdm_frame_log));
+  printf("av tx pwr: %4.2f PAPR: %4.2f av rx pwr: %4.2f\n", var(tx_fdm_frame_log), papr_dB, var(ch_fdm_frame_log));
 
   % some other useful plots
 
@@ -502,26 +514,43 @@ else
   title('tx fdm imag');
 
   figure(2)
-  clf;
-  plot(rx_symb_log*exp(j*pi/4),'+')
-  title('Scatter');
+  clf
+  spec = 20*log10(abs(fft(tx_fdm_frame_log)));
+  l = length(spec);
+  plot((Fs/l)*(1:l), spec)
+  axis([1 Fs/2 0 max(spec)]);
+  title('tx spectrum');
+  ylabel('Amplitude (dB)')
+  xlabel('Frequency (Hz)')
+  grid;
 
   figure(3)
   clf;
-  plot(rx_timing_log)
-  title('rx timing');
+  plot(rx_symb_log*exp(j*pi/4),'+')
+  title('Scatter');
 
   figure(4)
   clf;
-  stem(nerr_log)
-  title('Bit Errors');
+  subplot(211)
+  plot(rx_phi_log)
+  subplot(212)
+  plot(rx_amp_log)
 
   figure(5)
   clf;
+  subplot(211)
+  plot(rx_timing_log)
+  title('rx timing');
+  subplot(212)
   stem(ratio_log)
   title('Sync ratio');
 
-  figure(6);
+  figure(6)
+  clf;
+  stem(nerr_log)
+  title('Bit Errors');
+
+  figure(7);
   clf;
   subplot(211)
   plot(foff_log,';freq offset;');
@@ -534,15 +563,6 @@ else
   plot(foff_log(1:length(f_est_log)) - f_est_log + Fcentre)
   title('freq offset estimation error');
 
-  figure(7)
-  clf
-  spec = 20*log10(abs(fft(tx_fdm_frame_log)));
-  l = length(spec);
-  plot((Fs/l)*(1:l), spec)
-  axis([1 Fs/2 0 max(spec)]);
-  title('tx spectrum');
-  ylabel('Amplitude (dB)')
-  xlabel('Frequency (Hz)')
 
 end
 

diff --git a/codec2-dev/octave/test_cohpsk.m b/codec2-dev/octave/test_cohpsk.m
index 3b0abed7c633a023b23d76005f08a85ea967402f..6864488ebaf19cba25254ff3e8f6bca5d6600dfc 100644 (file)
--- a/codec2-dev/octave/test_cohpsk.m
+++ b/codec2-dev/octave/test_cohpsk.m
@@ -554,7 +554,7 @@ endfunction
 more off;
 %close all;
 %test_curves();
-%test_single();
+test_single();
 %rate_Fs_tx("tx_zero.raw");
 %rate_Fs_tx("tx.raw");
 %rate_Fs_rx("tx_-4dB.wav")

diff --git a/codec2-dev/octave/test_foff.m b/codec2-dev/octave/test_foff.m
index ee5df9ddda915ff0acd78f5e5d900995c4487a9b..3d26e89e0bf1ee092bdb3dc3a73d2268ad3b9518 100644 (file)
--- a/codec2-dev/octave/test_foff.m
+++ b/codec2-dev/octave/test_foff.m
@@ -91,7 +91,7 @@ function sim_out = freq_off_est_test(sim_in)
   rand('state',1); \r
   tx_bits_coh = round(rand(1,framesize*10));\r
   ptx_bits_coh = 1;\r
-  [spread spread_2ms hf_gain] = init_hf_model(Fs, Fs, frames*acohpsk.Nsymbrowpilot*afdmdv.M);\r
+  [spread spread_2ms hf_gain] = init_hf_model(Fs, frames*acohpsk.Nsymbrowpilot*afdmdv.M);\r
 \r
   hf_n = 1;\r
   nhfdelay = floor(hf_delay_ms*Fs/1000);\r