Finally have my head on straight about how to go about ofdm timing/freq est
authorbaobrien <baobrien@01035d8c-6547-0410-b346-abe4f91aad63>
Tue, 23 Jan 2018 19:39:19 +0000 (19:39 +0000)
committerbaobrien <baobrien@01035d8c-6547-0410-b346-abe4f91aad63>
Tue, 23 Jan 2018 19:39:19 +0000 (19:39 +0000)
git-svn-id: https://svn.code.sf.net/p/freetel/code@3392 01035d8c-6547-0410-b346-abe4f91aad63

codec2-dev/octave/ofdm_lib.m
codec2-dev/octave/tofdm.m
codec2-dev/src/fmfsk.c
codec2-dev/src/ofdm.c

index a618bda35e4181ad56d917d7c9ea4044faea6cb2..b04c9a3065a56fba27384d8fc58547a5f1f5d76b 100644 (file)
@@ -168,6 +168,7 @@ function states = ofdm_init(bps, Rs, Tcp, Ns, Nc)
   states.sample_point = states.timing_est = 1;
   states.nin = states.Nsamperframe;
 
+  states.sync = 0;
   % generate OFDM pilot symbol, used for timing and freq offset est
 
   rate_fs_pilot_samples = states.pilots * W/states.M;
@@ -430,162 +431,19 @@ endfunction
 
 function [rx_bits states aphase_est_pilot_log rx_np rx_amp] = ofdm_demod2(states, rxbuf_in)
   ofdm_load_const;
-
-  % insert latest input samples into rxbuf
-
-  rxbuf(1:Nrxbuf-states.nin) = rxbuf(states.nin+1:Nrxbuf);
-  rxbuf(Nrxbuf-states.nin+1:Nrxbuf) = rxbuf_in;
-
-  % get latest freq offset estimate
-
-  woff_est = 2*pi*foff_est_hz/Fs;
-
-  % update timing estimate --------------------------------------------------
-
-  delta_t = 0;
-  if timing_en
-    % update timing at start of every frame
-
-    st = M+Ncp + Nsamperframe + 1 - floor(ftwindow_width/2) + (timing_est-1);
-    en = st + Nsamperframe-1 + M+Ncp + ftwindow_width-1;
-          
-    ft_est = coarse_sync(states, rxbuf(st:en) .* exp(-j*woff_est*(st:en)), rate_fs_pilot_samples);
-    timing_est = timing_est + ft_est - ceil(ftwindow_width/2);
-
-    if verbose > 1
-      printf("  ft_est: %2d timing_est: %2d sample_point: %2d\n", ft_est, timing_est, sample_point);
+  [rx_bits states aphase_est_pilot_log rx_np rx_amp] =ofdm_demod(states,rxbuf_in);
+
+  %If out of sync, try a coarse sync on RX buffer
+  if states.sync <= 0
+    [t_est foff_est] = coarse_sync(states,states.rxbuf,rate_fs_pilot_samples);
+    %Correct to be within one frame and at head of frame instead of after first symbol
+    t_est = mod(t_est,Nsamperframe) - (M+Ncp);
+    if(t_est > 30)
+      states.nin = t_est;
     end
 
-    % Black magic to keep sample_point inside cyclic prefix.  Or something like that.
-
-    delta_t = ft_est - ceil(ftwindow_width/2);
-    sample_point = max(timing_est+Ncp/4, sample_point);
-    sample_point = min(timing_est+Ncp, sample_point);
+    printf("t_est is %d, foff_est is %d\n",t_est,foff_est)
   end
-
-  % down convert at current timing instant----------------------------------
-
-  % todo: this cld be more efficent, as pilot r becomes r-Ns on next frame
-
-  rx_sym = zeros(1+Ns+1+1, Nc+2);
-
-  % previous pilot
-  
-  st = M+Ncp + Nsamperframe + (-Ns)*(M+Ncp) + 1 + sample_point; en = st + M - 1;
-
-  for c=1:Nc+2
-    acarrier = rxbuf(st:en) .* exp(-j*woff_est*(st:en)) .* conj(W(c,:));
-    rx_sym(1,c) = sum(acarrier);
-  end
-
-  % pilot - this frame - pilot
-
-  for rr=1:Ns+1 
-    st = M+Ncp + Nsamperframe + (rr-1)*(M+Ncp) + 1 + sample_point; en = st + M - 1;
-    for c=1:Nc+2
-      acarrier = rxbuf(st:en) .* exp(-j*woff_est*(st:en)) .* conj(W(c,:));
-      rx_sym(rr+1,c) = sum(acarrier);
-    end
-  end
-
-  % next pilot
-
-  st = M+Ncp + Nsamperframe + (2*Ns)*(M+Ncp) + 1 + sample_point; en = st + M - 1;
-  for c=1:Nc+2
-    acarrier = rxbuf(st:en) .* exp(-j*woff_est*(st:en)) .* conj(W(c,:));
-    rx_sym(Ns+3,c) = sum(acarrier);
-  end
-      
-  % est freq err based on all carriers ------------------------------------
-      
-  if foff_est_en
-    freq_err_rect = sum(rx_sym(2,:))' * sum(rx_sym(2+Ns,:));
-
-    % prevent instability in atan(im/re) when real part near 0 
-
-    freq_err_rect += 1E-6;
-
-    %printf("freq_err_rect: %f %f angle: %f\n", real(freq_err_rect), imag(freq_err_rect), angle(freq_err_rect));
-    freq_err_hz = angle(freq_err_rect)*Rs/(2*pi*Ns);
-    foff_est_hz = foff_est_hz + foff_est_gain*freq_err_hz;
-  end
-
-  % OK - now estimate and correct phase  ----------------------------------
-
-  aphase_est_pilot = 10*ones(1,Nc+2);
-  aamp_est_pilot = zeros(1,Nc+2);
-  for c=2:Nc+1
-
-    % estimate phase using average of 6 pilots in a rect 2D window centred
-    % on this carrier
-    % PPP
-    % DDD
-    % DDD
-    % PPP
-          
-    cr = c-1:c+1;
-    aphase_est_pilot_rect = sum(rx_sym(2,cr)*pilots(cr)') + sum(rx_sym(2+Ns,cr)*pilots(cr)');
-
-    % use next step of pilots in past and future
-
-    aphase_est_pilot_rect += sum(rx_sym(1,cr)*pilots(cr)');
-    aphase_est_pilot_rect += sum(rx_sym(2+Ns+1,cr)*pilots(cr)');
-
-    aphase_est_pilot(c) = angle(aphase_est_pilot_rect);
-    aamp_est_pilot(c) = abs(aphase_est_pilot_rect/12);
-  end
-  % correct phase offset using phase estimate, and demodulate
-  % bits, separate loop as it runs across cols (carriers) to get
-  % frame bit ordering correct
-
-  aphase_est_pilot_log = [];
-  rx_bits = []; rx_np = []; rx_amp = [];
-  for rr=1:Ns-1
-    for c=2:Nc+1
-      if phase_est_en
-        rx_corr = rx_sym(rr+2,c) * exp(-j*aphase_est_pilot(c));
-      else
-        rx_corr = rx_sym(rr+2,c);
-      end
-      rx_np = [rx_np rx_corr];
-      rx_amp = [rx_amp aamp_est_pilot(c)];
-      if bps == 1
-        abit = real(rx_corr) > 0;
-      end
-      if bps == 2
-        abit = qpsk_demod(rx_corr);
-      end
-      rx_bits = [rx_bits abit];
-    end % c=2:Nc+1
-    aphase_est_pilot_log = [aphase_est_pilot_log; aphase_est_pilot(2:Nc+1)];
-  end 
-
-  % Adjust nin to take care of sample clock offset
-
-  nin = Nsamperframe;
-  if timing_en
-    thresh = (M+Ncp)/8;
-    tshift = (M+Ncp)/4;
-    if timing_est > thresh
-      nin = Nsamperframe+tshift;
-      timing_est -= tshift;
-      sample_point -= tshift;
-    end
-    if timing_est < -thresh
-      nin = Nsamperframe-tshift;
-      timing_est += tshift;
-      sample_point += tshift;
-    end
-  end
-
-  states.rx_sym = rx_sym;
-  states.rxbuf = rxbuf;
-  states.nin = nin;
-  states.timing_est = timing_est;
-  states.sample_point = sample_point;
-  states.delta_t = delta_t;
-  states.foff_est_hz = foff_est_hz;
 endfunction
 
 
index 88a6dd6676c31d496077fd757e744e1274c5acfd..c0177a061b40fad5f31f24da4b4dbf0f9a027eb2 100644 (file)
@@ -10,8 +10,9 @@ function rx = ebno_awgn_channel(tx, Fs, Rb, EbNodB)
   EbNo = 10^(EbNodB/10);
   variance = Fs / (Rb*EbNo);
   noise = sqrt(variance)*randn(1,length(tx));
-  avg = sum(abs(tx))/length(tx)
+  avg = sum(abs(tx))/length(tx);
   rx = noise.*avg + tx;
+  %rx = rx * .10;
 end
 
 function nums = im_re_interleave(nim)
@@ -197,4 +198,4 @@ end
 % This works best on my machine -- Brady
 graphics_toolkit('fltk')
 
-run_ofdm_test(60,100,.5,60 ,0)
\ No newline at end of file
+run_ofdm_test(60,20,0.5,40,100)
\ No newline at end of file
index ab56945fbba0a446cc50752b617d95d2d48adac1..2b52335241fd31b8a40613a4fb31ace891768350 100644 (file)
@@ -198,8 +198,8 @@ void fmfsk_demod(struct FMFSK *fmfsk, uint8_t rx_bits[],float fmfsk_in[]){
     #endif
     
     /* Shift in nin samples */
-    memcpy(&oldsamps[0]   , &oldsamps[nmem-nold], sizeof(float)*nold);
-    memcpy(&oldsamps[nold], &fmfsk_in[0]        , sizeof(float)*nin );
+    memmove(&oldsamps[0]   , &oldsamps[nmem-nold], sizeof(float)*nold);
+    memcpy (&oldsamps[nold], &fmfsk_in[0]        , sizeof(float)*nin );
     
     /* Allocate memory for filtering */
     float *rx_filt = alloca(sizeof(float)*(nsym+1)*Ts);
index bcca1effc164eef99d21db41d59b492419796739..99734d27c31b233b186da456fcc93c5f27bda4d2 100644 (file)
@@ -192,7 +192,7 @@ static int coarser_sync(struct OFDM *ofdm, complex float *rx, int length, float
     float max_mag = 0;
     float mag = 0.0f;
     int t_est = 0;
-    int fmax = 60;
+    int fmax = 20;
     int pmax_i,nmax_i;
     float pmax,nmax;
     float foff_est;
@@ -278,7 +278,7 @@ static int coarser_sync(struct OFDM *ofdm, complex float *rx, int length, float
     }
     foff_est = (pmax > nmax) ? pmax_i : (nmax_i-Fs+1); 
 
-    fprintf(stderr,"foff_est is %f, tpos is %d, ratio is %f\n",foff_est,t_est,max_mag/max_corr);
+    //fprintf(stderr,"foff_est is %f, tpos is %d, ratio is %f\n",foff_est,t_est,max_corr/max_mag);
 
     *foff_out = foff_est;
 
@@ -388,7 +388,7 @@ static int coarse_sync(struct OFDM *ofdm, complex float *rx, int length, float *
     }
     foff_est = (pmax > nmax) ? pmax_i : (nmax_i-Fs+1); 
     //foff_est = pmax_i;
-    //`fprintf(stderr,"foff_est is %f, tpos is %d\n",foff_est,t_est);
+    fprintf(stderr,"foff_est is %f, tpos is %d, ratio is %f\n",foff_est,t_est%1280,max_corr/max_mag);
 
     *foff_out = foff_est;
 
@@ -778,6 +778,7 @@ void ofdm_mod(struct OFDM *ofdm, COMP result[], const int *tx_bits) {
 /*
  * Like ofdm_demod, but handles sync and large frame offset
  */
+//static int coarser_sync(struct OFDM *ofdm, complex float *rx, int length, float *foff_out, float *ratio_out)
 void ofdm_demod_coarse(struct OFDM *ofdm, int *rx_bits, COMP *rxbuf_in){
     const int SampsPerFrame       = ofdm->config.SampsPerFrame;
     const int RxBufSize       = ofdm->config.RxBufSize;
@@ -789,16 +790,42 @@ void ofdm_demod_coarse(struct OFDM *ofdm, int *rx_bits, COMP *rxbuf_in){
     complex float * coarse_rxbuf = ofdm->coarse_rxbuf;
     complex float rxbuf_temp[RxBufSize];
 
+
+    /* Copy in nin samples to end of buffer */
+    memmove(&coarse_rxbuf[0],&coarse_rxbuf[nin], (RxBufSize - nin) * sizeof(coarse_rxbuf[0]));
+    memcpy (&coarse_rxbuf[RxBufSize - nin],rxbuf_in,    nin * sizeof(coarse_rxbuf[0]));
+
+    int frame_pos_est;
+    float fshift;
+    float corr_ratio_max = 0;
+    float corr_ratio;
+    float foff_out;
+    float fshift_max;
+    int frame_pos_max = 0;
     /* NCO for freq. shifting */
     complex float shift_nco = 1;
     complex float shift_nco_dph = 0;
 
-    /* Copy in nin samples to end of buffer */
-    memmove(&coarse_rxbuf[0],&coarse_rxbuf[nin],nin);
-    memcpy (&coarse_rxbuf[RxBufSize - nin],rxbuf_in,nin);
-
-    if(!sync) {
+    if(sync <= 0 ) {
+        /* Do coarse search over frequency range in 20hz slices*/
+        for(fshift = -100; fshift < 100; fshift+=40){
+            /* Shift block of samples by fshift */
+            shift_nco_dph = cexpf(2*M_PI*(fshift/(float)Fs)*I);
+            shift_nco = 1;
+            for(int i=0; i<RxBufSize; i++){
+                rxbuf_temp[i] = coarse_rxbuf[i] * shift_nco;
+                shift_nco = shift_nco * shift_nco_dph;
+            }
 
+            /* Do a coarse search for the pilot */
+            frame_pos_est = coarser_sync(ofdm,rxbuf_temp,RxBufSize,&foff_out,&corr_ratio);
+            if(corr_ratio > corr_ratio_max){
+                corr_ratio_max = corr_ratio;
+                fshift_max = -fshift + foff_out;
+                frame_pos_max = frame_pos_est;          
+            }
+        }
+        fprintf(stderr,"Best ratio %f freq %f offset %d\n",corr_ratio_max,fshift_max,frame_pos_max);
     }
 }
 
@@ -832,21 +859,21 @@ void ofdm_demod(struct OFDM *ofdm, int *rx_bits, COMP *rxbuf_in) {
     complex float * rxbuf = ofdm->rxbuf;
 
 
-    //memmove(&rxbuf[0],&rxbuf[nin],nin);
-    //memcpy (&rxbuf[RxBufSize - nin],rxbuf_in,nin);
+    memmove(&rxbuf[0],&rxbuf[nin], (RxBufSize - nin) * sizeof(rxbuf[0]));
+    memcpy (&rxbuf[RxBufSize - nin],rxbuf_in,    nin * sizeof(rxbuf[0]));
 
     /* shift the buffer left based on nin */
 
-    for (i = 0, j = ofdm->nin; i < (RxBufSize - ofdm->nin); i++, j++) {
-        ofdm->rxbuf[i] = ofdm->rxbuf[j];
-    }
+    // for (i = 0, j = ofdm->nin; i < (RxBufSize - ofdm->nin); i++, j++) {
+    //     ofdm->rxbuf[i] = ofdm->rxbuf[j];
+    // }
 
     /* insert latest input samples onto tail of rxbuf */
     /* Note: COMP and complex float have the same memory layout */
     /* see iso/iec 9899:1999 sec 6.2.5.13 */
-    for (i = (RxBufSize - ofdm->nin), j = 0; i < RxBufSize; i++, j++) {
-        ofdm->rxbuf[i] = rxbuf_in[j].real + rxbuf_in[j].imag * I;
-    }
+    // for (i = (RxBufSize - ofdm->nin), j = 0; i < RxBufSize; i++, j++) {
+    //     ofdm->rxbuf[i] = rxbuf_in[j].real + rxbuf_in[j].imag * I;
+    // }
 
 
     /*
@@ -876,9 +903,12 @@ void ofdm_demod(struct OFDM *ofdm, int *rx_bits, COMP *rxbuf_in) {
             work[j] = ofdm->rxbuf[i] * cexpf(-I * woff_est * i);
         }
 
+        fprintf(stderr,"n>");
         ft_est = coarse_sync(ofdm, work, (en - st), &freq_offset);
-
-        coarser_sync(ofdm, work, (en - st), &freq_offset, NULL);
+        fprintf(stderr,"\nx>");
+        coarse_sync(ofdm,ofdm->rxbuf,RxBufSize,&freq_offset);
+        fprintf(stderr,"\n");
+        //coarser_sync(ofdm, work, (en - st), &freq_offset, NULL);
 
         //ofdm->foff_est_hz += (freq_offset/2);
         ofdm->timing_est += (ft_est - ceilf(FtWindowWidth / 2));