% Octave script for comparing Octave and C versions of OFDZM modem
% ------------------------------------------------------------------
-1;
+
+load rand_state_vec.mat
+%rand("state",rand_state)
function rx = ebno_awgn_channel(tx, Fs, Rb, EbNodB)
+ %rand("seed",2)
+ %randn("seed",3);
EbNo = 10^(EbNodB/10);
variance = Fs / (Rb*EbNo);
noise = sqrt(variance)*randn(1,length(tx));
%rx = rx * .10;
end
+function rx = snr_awgn_channel(tx, Fs, BwS, snrdb)
+ snr = 10^(snrdb/10);
+ %signal power
+ psig = sum(abs(tx))^2 / length(tx);
+ %noise power in interesting band
+ pnoiseb = (1/snr) * psig;
+ %noise power in 1hz band
+ pnoise0 = pnoiseb/BwS;
+ %noise power over sample rate bandwidth
+ pnoiset = pnoise0 * (Fs/2);
+
+ noise = sqrt(pnoiset)*randn(1,length(tx));
+
+ rx = noise + tx;
+ %rx = rx * .10;
+end
+
function nums = im_re_interleave(nim)
nums = zeros(1,length(nim)*2);
nums(1:2:length(nums)) = real(nim);
endfunction
function pass = run_ofdm_test(Nframes,sample_clock_offset_ppm,foff_hz,EbNodB = 100,nheadsamp = 0)
+ rand("seed",2)
+ randn("seed",2)
%Nframes = 30;
%sample_clock_offset_ppm = 100;
%foff_hz = 5;
states = ofdm_init(bps, Rs, Tcp, Ns, Nc);
ofdm_load_const;
- rand('seed',2);
tx_bits = round(rand(1,Nbitsperframe));
% Run tx loop
end
[rx_bits states aphase_est_pilot_log arx_np arx_amp] = ofdm_demod2(states, rxbuf_in);
-
% log some states for comparison to C
rxbuf_in_log = [rxbuf_in_log rxbuf_in];
end
-% This works best on my machine -- Brady
-graphics_toolkit('fltk')
+function [ber berc] = run_ber_test(Nframes,sample_clock_offset_ppm,foff_hz,EbNodB = 100,nheadsamp = 0)
+ rand("seed",2)
+ randn("seed",2)
+ %Nframes = 30;
+ %sample_clock_offset_ppm = 100;
+ %foff_hz = 5;
+
+ more off; format;
+ ofdm_lib;
+ autotest;
+
+ % ---------------------------------------------------------------------
+ % Run Octave version
+ % ---------------------------------------------------------------------
+
+ Ts = 0.018; Tcp = 0.002; Rs = 1/Ts; bps = 2; Nc = 16; Ns = 8;
+ Rb_real = 1400.0;
+ states = ofdm_init(bps, Rs, Tcp, Ns, Nc);
+ ofdm_load_const;
+
+ tx_bits = round(rand(1,Nbitsperframe));
+
+ % Run tx loop
+
+ tx_bits_log = []; tx_log = [];
+ for f=1:Nframes
+ tx_bits_log = [tx_bits_log tx_bits];
+ tx_log = [tx_log ofdm_mod(states, tx_bits)];
+ end
+
+ % Channel simulation ----------------------------------------------
+
+ rx_log = sample_clock_offset(tx_log, sample_clock_offset_ppm);
+ rx_log = freq_shift(rx_log, foff_hz, Fs);
+ #rx_log(1:nheadsamp) = zeros(nheadsamp,1);
+ rx_log_l = length(rx_log)
+ rx_log = [zeros(1,nheadsamp), rx_log];
+ rx_log = rx_log(1:rx_log_l);
+ %rx_log = snr_awgn_channel(rx_log, Fs, Fs, EbNodB);
+ rx_log = ebno_awgn_channel(rx_log, Fs, Rb_real, EbNodB);
+
+ rx_vec = fopen("tofdm_rx_vec","wb+");
+ fwrite(rx_vec,im_re_interleave(rx_log),"single");
+ fclose(rx_vec);
+
+ % Rx ---------------------------------------------------------------
+
+ % Init rx with ideal timing so we can test with timing estimation disabled
+
+ Nsam = length(rx_log);
+ prx = 1;
+ nin = Nsamperframe+2*(M+Ncp);
+ %states.rxbuf(Nrxbuf-nin+1:Nrxbuf) = rx_log(prx:nin);
+ %prx += nin;
+
+ rxbuf_log = []; rxbuf_in_log = []; rx_sym_log = []; foff_hz_log = [];
+ timing_est_log = []; sample_point_log = [];
+ phase_est_pilot_log = []; rx_amp_log = [];
+ rx_np_log = []; rx_bits_log = [];
+
+ states.timing_en = 1;
+ states.foff_est_en = 1;
+ states.phase_est_en = 1;
+
+ if states.timing_en == 0
+ % manually set ideal timing instant
+ states.sample_point = Ncp;
+ end
+
+ have_frame = true;
+ while have_frame
+ % insert samples at end of buffer, set to zero if no samples
+ % available to disable phase estimation on future pilots on last
+ % frame of simulation
+
+ nin = states.nin;
+ lnew = min(Nsam-prx+1,nin);
+ rxbuf_in = zeros(1,nin);
+ %printf("nin: %d prx: %d lnew: %d\n", nin, prx, lnew);
+ if lnew
+ rxbuf_in(1:lnew) = rx_log(prx:prx+lnew-1);
+ end
+ prx += nin;
+ if prx > Nsam
+ have_frame = false;
+ end
+
+ [rx_bits states aphase_est_pilot_log arx_np arx_amp] = ofdm_demod(states, rxbuf_in);
+ % log some states for comparison to C
+
+ rxbuf_in_log = [rxbuf_in_log rxbuf_in];
+ rxbuf_log = [rxbuf_log states.rxbuf];
+ rx_sym_log = [rx_sym_log; states.rx_sym];
+ phase_est_pilot_log = [phase_est_pilot_log; aphase_est_pilot_log];
+ rx_amp_log = [rx_amp_log arx_amp];
+ foff_hz_log = [foff_hz_log; states.foff_est_hz];
+ timing_est_log = [timing_est_log; states.timing_est];
+ sample_point_log = [sample_point_log; states.sample_point];
+ rx_np_log = [rx_np_log arx_np];
+ rx_bits_log = [rx_bits_log rx_bits];
+
+ end
+
+ % ---------------------------------------------------------------------
+ % Run C version and plot Octave and C states and differences
+ % ---------------------------------------------------------------------
+
+ % Override default path by setting path_to_tofdm = "/your/path/to/tofdm"
+
+ if exist("path_to_tofdm", "var") == 0
+ path_to_tofdm = "../build_linux/unittest/tofdm";
+ end
+ system(path_to_tofdm);
+
+ load tofdm_out.txt;
+ % Generated with modem probe thing
+ load ofdm_test.txt;
+
+ front_cutoff = 10000;
+ rx_np_log = rx_np_log(front_cutoff/2:length(rx_np_log));
+ rx_np_log_c = rx_np_log_c(front_cutoff/2:length(rx_np_log_c));
+
+ fg = 1;
+ figure(fg++); clf; plot(foff_hz_log_c); title('foff hz C');
+ figure(fg++); clf; plot(rx_np_log,'+'); title('Octave Scatter Diagram'); axis([-1.5 1.5 -1.5 1.5]);
+ figure(fg++); clf; plot(rx_np_log_c,'+'); title('C Scatter Diagram'); axis([-1.5 1.5 -1.5 1.5]);
+
+
+ [rx_bits_log rx_bits_log_c] = trim1(rx_bits_log,rx_bits_log_c);
+ [tx_bits_log rx_bits_log] = trim1(tx_bits_log,rx_bits_log);
+ [tx_bits_log rx_bits_log_c] = trim1(tx_bits_log,rx_bits_log_c);
+
+ bitcnt = length(rx_bits_log_c);
+ rx_bits = rx_bits_log_c(front_cutoff:length(rx_bits_log_c));
+ tx_bits = tx_bits_log(front_cutoff:length(tx_bits_log));
+
+ ber = 1;
+ for offset = (1:100)
+ nerr = sum(xor(rx_bits(offset:length(rx_bits)),tx_bits(1:length(rx_bits)+1-offset)));
+ bern = nerr/(bitcnt-offset);
+ if(bern < ber)
+ ox = offset;
+ best_nerr = nerr;
+ end
+ ber = min([ber bern]);
+ end
+ offset = ox;
+ berc = ber;
+
+ figure(fg++)
+ plot(xor(rx_bits(offset:length(rx_bits)),tx_bits(1:length(rx_bits)+1-offset)))
+ title('C bit errors')
+
+ rx_bits = rx_bits_log(front_cutoff:length(rx_bits));
+ ber = 1;
+
+ for offset = (1:100)
+ nerr = sum(xor(rx_bits(offset:length(rx_bits)),tx_bits(1:length(rx_bits)+1-offset)));
+ bern = nerr/(bitcnt-offset);
+ if(bern < ber)
+ ox = offset;
+ best_nerr = nerr;
+ end
+ ber = min([ber bern]);
+ end
+ offset = ox;
+ ber = ber;
+
+ figure(fg++)
+ plot(xor(rx_bits(offset:length(rx_bits)),tx_bits(1:length(rx_bits)+1-offset)))
+ title('Octave bit errors')
+
+
+ figure(fg++)
+ plot(20*log10(abs(fft(rx_log))), (1:length(rx_log)) .* (Fs/length(rx_log) ))
+
+endfunction
+
+% FLTK plots work better than gnuplot or the default on my machine -- Brady
+try
+ graphics_toolkit('fltk')
+catch err
+ warning(err.identifier, err.message);
+end_try_catch
-run_ofdm_test(60,100,5,8,10000)
+%run_ofdm_test(60,100,60,5,10000)
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;
}
}
+ if(ratio_out != NULL){
+ *ratio_out = (max_corr + 1e-9)/(max_mag + 1e-9);
+ }
+
/* Skip coarse frequency est if we don't have anywhere to put result */
if(foff_out == NULL){
return 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;
- if(ratio_out != NULL){
- *ratio_out = (max_corr + 1e-9)/(max_mag + 1e-9);
- }
-
return t_est;
}
complex float shift_nco = 1;
complex float shift_nco_dph = 0;
- if(sync <= 0 || 1 ) {
+ if(sync <= 0 ) {
/* Do coarse search over frequency range in 20hz slices*/
- //for(fshift = -100; fshift < 100; fshift+=40){
- fshift = 0;{
+ for(fshift = -103; fshift < 100; fshift+=40){
+ //fshift = 0;{
/* Shift block of samples by fshift */
shift_nco_dph = cexpf(2*M_PI*(fshift/(float)Fs)*I);
shift_nco = 1;
}
/* Do a coarse search for the pilot */
- frame_pos_est = coarse_sync(ofdm,rxbuf_temp,RxBufSize,&foff_out,&corr_ratio);
+ frame_pos_est = coarse_sync(ofdm,rxbuf_temp,RxBufSize,NULL,&corr_ratio);
if(corr_ratio > corr_ratio_max){
corr_ratio_max = corr_ratio;
- fshift_max = -fshift + foff_out;
+ fshift_max = fshift;
frame_pos_max = frame_pos_est;
}
}
frame_pos_est = (frame_pos_max % SampsPerFrame) - (M+Ncp) + 1;
if(corr_ratio_max > 0.5 && sync <=0){
- ofdm->sync_count = 3;
+ /* Now that we think we can sync, do a freq. domain search in the smaller range for our signal */
+ shift_nco_dph = cexpf(2*M_PI*(fshift_max/(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;
+ }
+ coarse_sync(ofdm,rxbuf_temp,RxBufSize,&foff_out,&corr_ratio);
+
+ ofdm->sync_count = 2;
ofdm->frame_point = frame_pos_est;
- ofdm->foff_est_hz = fshift_max;
- fprintf(stderr,"syncing\n");
+ ofdm->foff_est_hz = -fshift_max + foff_out;
+ fprintf(stderr,"syncing, shift %f,fine foff %f\n",fshift_max,-fshift_max + foff_out);
}
- fprintf(stderr,"Best ratio %f freq %f offset %d\n",corr_ratio_max,fshift_max,frame_pos_est);
+ fprintf(stderr,"Best ratio %f freq %f offset %d fine %f\n",corr_ratio_max,fshift_max,frame_pos_est,ofdm->foff_est_hz);
}
}
ft_est = coarse_sync(ofdm, work, (en - st), &freq_offset,&ratio_out);
fprintf(stderr,"good ratio %f\n",ratio_out);
// fprintf(stderr,"\nx>");
-
+ if(ratio_out < .7){
+ ofdm->sync_count--;
+ }
//ofdm->foff_est_hz += (freq_offset/2);
ofdm->timing_est += (ft_est - ceilf(FtWindowWidth / 2));
ft_est_other = (ft_est_other % SampsPerFrame) - (M+Ncp) + 1;
// fprintf(stderr,"\n");
//coarser_sync(ofdm, work, (en - st), &freq_offset, NULL);
- if(ratio_out<.5){
- ofdm->sync_count--;
+
+ fprintf(stderr,"ratio is %f t_offset is %d f_off is %f freq_off is %f\n",ratio_out,ft_est_other,ofdm->foff_est_hz,freq_offset);
+ if(freq_offset > 1){
+ //ofdm->foff_est_hz += freq_offset/4;
}
- fprintf(stderr,"ratio is %f t_offset is %d\n",ratio_out,ft_est_other);
}
/*
projects / freetel-svn-tracking.git / commitdiff