% returns an estimate of frequency offset, advances to next sync state
function [next_sync cohpsk] = coarse_freq_offset_est(cohpsk, fdmdv, ch_fdm_frame, sync, next_sync)
- Fcentre = fdmdv.Fcentre;
- Nc = fdmdv.Nc;
- Fsep = fdmdv.Fsep;
- M = fdmdv.M;
- Fs = fdmdv.Fs;
- Ndft = cohpsk.Ndft;
+ Fcentre = fdmdv.Fcentre;
+ Nc = fdmdv.Nc;
+ Fsep = fdmdv.Fsep;
+ M = fdmdv.M;
+ Fs = fdmdv.Fs;
+ Ndft = cohpsk.Ndft;
+ coarse_mem = cohpsk.coarse_mem;
+ Ncm = cohpsk.Ncm;
+
+ % ll
+ % |--------|-----|
+
+ ll = length(ch_fdm_frame);
+ sz_mem = Ncm-ll;
+ for i=1:sz_mem
+ coarse_mem(i) = coarse_mem(i+ll);
+ end
+ coarse_mem(Ncm-ll+1:Ncm) = ch_fdm_frame;
if sync == 0
- f_start = Fcentre - ((Nc/2)+2)*Fsep;
- f_stop = Fcentre + ((Nc/2)+2)*Fsep;
- ll = length(ch_fdm_frame);
- h = 0.5 - 0.5*cos(2*pi*(0:ll-1)/(ll-1));
- T = abs(fft(ch_fdm_frame .* h, Ndft)).^2;
+ h = 0.5 - 0.5*cos(2*pi*(0:Ncm-1)/(Ncm-1));
+ T = abs(fft(coarse_mem .* h, Ndft)).^2;
sc = Ndft/Fs;
- bin_start = floor(f_start*sc+0.5)+1;
- bin_stop = floor(f_stop*sc+0.5)+1;
+ for i=1:5
+ f_start = Fcentre - ((Nc/2)+2)*Fsep;
+ f_stop = Fcentre + ((Nc/2)+2)*Fsep;
+ bin_start = floor(f_start*sc+0.5)+1;
+ bin_stop = floor(f_stop*sc+0.5)+1;
+ x = bin_start-1:bin_stop-1;
+ bin_est = x*T(bin_start:bin_stop)'/sum(T(bin_start:bin_stop));
+ f_est = floor(bin_est/sc+0.5);
+ Fcentre = f_est;
+ end
%printf("f_start: %f f_stop: %f sc: %f bin_start: %d bin_stop: %d\n", f_start, f_stop, sc, bin_start, bin_stop);
- x = bin_start-1:bin_stop-1;
- bin_est = x*T(bin_start:bin_stop)'/sum(T(bin_start:bin_stop));
- cohpsk.f_est = floor(bin_est/sc+0.5);
+ cohpsk.f_est = f_est;
+
printf(" coarse freq est: %f\n", cohpsk.f_est);
next_sync = 1;
-
+ figure(5)
+ clf
+ subplot(211)
+ plot(T)
+ hold on
+ plot([bin_est bin_est],[0 max(T)],'g')
+ hold off
+ axis([bin_start bin_stop 0 max(T)])
+
end
+ cohpsk.coarse_mem = coarse_mem;
endfunction
% sample correlation over 2D grid of time and fine freq points
max_corr = 0;
- for f_fine=-20:1:20
+ for f_fine=-20:0.25:20
f_fine_rect = exp(-j*f_fine*2*pi*sampling_points/Rs)';
for t=0:cohpsk.Nsymbrowpilot-1
corr = 0; mag = 0;
end
printf(" fine freq f: %f max_corr: %f max_mag: %f ct: %d\n", cohpsk.f_fine_est, abs(max_corr), max_mag, cohpsk.ct);
- if max_corr/max_mag > 0.9
- printf(" in sync!\n");
+ if abs(max_corr/max_mag) > 0.7
+ printf(" [%d] in sync!\n", cohpsk.frame);
+ cohpsk.sync_timer = 0;
+ %cohpsk.f_est -= cohpsk.f_fine_est;
+ %cohpsk.f_fine_est = 0;
+ %cohpsk.ff_rect = 1;
+ printf(" .... adjusting to %f\n", cohpsk.f_est);
next_sync = 4;
else
next_sync = 0;
printf(" back to coarse freq offset est...\n");
end
+ cohpsk.ratio = abs(max_corr/max_mag);
+ end
+
+
+ if sync == 4
+
+ % we are in sync so just sample correlation over 1D grid of fine freq points
+
+ max_corr = 0;
+ st = cohpsk.f_fine_est - 1;
+ en = cohpsk.f_fine_est + 1;
+ for f_fine = st:0.25*en
+ f_fine_rect = exp(-j*f_fine*2*pi*sampling_points/Rs)';
+ corr = 0; mag = 0;
+ for c=1:Nc*Nd
+ f_corr_vec = f_fine_rect .* ct_symb_buf(cohpsk.ct+sampling_points,c);
+ for p=1:length(sampling_points)
+ corr += pilot2(p,c-Nc*floor((c-1)/Nc)) * f_corr_vec(p);
+ mag += abs(f_corr_vec(p));
+ end
+ end
+ if corr > max_corr
+ max_corr = corr;
+ max_mag = mag;
+ f_fine_est = f_fine;
+ end
+ end
+
+ %cohpsk.f_est -= 0.5*f_fine_est;
+ %printf(" coarse: %f fine: %f\n", cohpsk.f_est, f_fine_est);
+ cohpsk.ratio = abs(max_corr/max_mag);
end
+
+
endfunction
% havent fine freq corrected ct_symb_buf if next_sync == 4 correct
% all 8 if sync == 2 correct latest 6
- if (next_sync == 4) || (sync == 4)
- if (next_sync == 4) && (sync == 2)
+ if (next_sync == 4) && (sync == 2)
% first frame, we've just gotten sync so fine freq correct all Nsymbrowpilot+2 samples
acohpsk.ff_phase *= acohpsk.ff_rect';
ct_symb_ff_buf(r,:) *= acohpsk.ff_phase;
end
+ end
- else
-
+ if sync == 4
% second and subsequent frames, just fine freq correct the latest Nsymbrowpilot
ct_symb_ff_buf(1:2,:) = ct_symb_ff_buf(acohpsk.Nsymbrowpilot+1:acohpsk.Nsymbrowpilot+2,:);
acohpsk.ff_phase *= acohpsk.ff_rect';
ct_symb_ff_buf(r,:) *= acohpsk.ff_phase;
end
+ end
- end
+ mag = abs(acohpsk.ff_phase);
+ acohpsk.ff_phase /= mag;
- mag = abs(acohpsk.ff_phase);
- acohpsk.ff_phase /= mag;
+ acohpsk.ct_symb_ff_buf = ct_symb_ff_buf;
- acohpsk.ct_symb_ff_buf = ct_symb_ff_buf;
- end
endfunction
% misc sync state machine code, just wanted it in a function
-function sync = sync_state_machine(sync, next_sync)
+function [sync cohpsk] = sync_state_machine(cohpsk, sync, next_sync)
+
if sync == 1
next_sync = 2;
end
+ if sync == 5
+ next_sync = 4;
+ end
+
+ if sync == 4
+
+ % check that sync is still good, fall out of sync on consecutive bad frames */
+
+ if cohpsk.ratio < 0.5
+ cohpsk.sync_timer++;
+ else
+ cohpsk.sync_timer = 0;
+ end
+ %printf(" ratio: %f sync timer: %d\n", cohpsk.ratio, cohpsk.sync_timer);
+
+ if cohpsk.sync_timer == 5
+ printf(" lost sync ....\n");
+ next_sync = 0;
+ end
+ end
+
sync = next_sync;
endfunction
% tcohpsk.m
% David Rowe Oct 2014
%
-% Octave script that tests the C port of the coherent PSK modem. This
-% script loads the output of unittest/tcohpsk.c and compares it to the
-% output of the reference versions of the same functions written in
-% Octave.
+% Octave script that:
%
-
+% i) tests the C port of the coherent PSK modem. This script loads
+% the output of unittest/tcohpsk.c and compares it to the output of
+% the reference versions of the same functions written in Octave.
+%
+% or (ii) can optionally be used to run an Octave version of the cohpsk
+% modem to tune and develop it.
graphics_toolkit ("gnuplot");
+more off;
cohpsk;
fdmdv;
% test parameters ----------------------------------------------------------
-n = 840;
-frames = 35;
-foff = 10.5;
-EsNodB = 8;
+frames = 100;
+foff = 100;
+dfoff = 0;
+EsNodB = 10;
+fading_en = 0;
+hf_delay_ms = 2;
+compare_with_c = 0;
EsNo = 10^(EsNodB/10);
rx_baseband_log = [];
rx_fdm_frame_log = [];
ct_symb_ff_log = [];
+rx_timing_log = [];
+ratio_log = [];
-% BER measurement -------------------------------------------------------------
+% Channel modeling and BER measurement ----------------------------------------
rand('state',1);
tx_bits_coh = round(rand(1,framesize*10));
phase_ch = 1;
sync = 0;
-% Output vectors from C port ---------------------------------------------------
-
-load ../build_linux/unittest/tcohpsk_out.txt
+[spread spread_2ms hf_gain] = init_hf_model(Fs, 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);
% main loop --------------------------------------------------------------------
for i=1:frames
+ acohpsk.frame = i;
tx_bits = tx_bits_coh(ptx_bits_coh:ptx_bits_coh+framesize-1);
ptx_bits_coh += framesize;
if ptx_bits_coh > length(tx_bits_coh)
% Channel --------------------------------------------------------------------
%
- [ch_fdm_frame phase_ch] = freq_shift(tx_fdm_frame, foff, Fs, phase_ch);
-
+ %[ch_fdm_frame phase_ch] = freq_shift(tx_fdm_frame, foff, Fs, phase_ch);
+ %foff += dfoff*acohpsk.Nsymbrowpilot*M;
+
+ ch_fdm_frame = zeros(1,acohpsk.Nsymbrowpilot*M);
+ for i=1:acohpsk.Nsymbrowpilot*M
+ foff_rect = exp(j*2*pi*foff/Fs);
+ foff += dfoff;
+ phase_ch *= foff_rect;
+ ch_fdm_frame(i) = tx_fdm_frame(i) * phase_ch;
+ end
+ phase_ch /= abs(phase_ch);
+
+ if fading_en
+ ch_fdm_delay(1:nhfdelay) = ch_fdm_delay(acohpsk.Nsymbrowpilot*M+1:nhfdelay+acohpsk.Nsymbrowpilot*M);
+ ch_fdm_delay(nhfdelay+1:nhfdelay+acohpsk.Nsymbrowpilot*M) = ch_fdm_frame;
+
+ for i=1:acohpsk.Nsymbrowpilot*M
+ ahf_model = hf_gain*(spread(hf_n)*ch_fdm_frame(i) + spread_2ms(hf_n)*ch_fdm_delay(i));
+ ch_fdm_frame(i) = ahf_model;
+ hf_n++;
+ end
+ end
+
% each carrier has power = 2, total power 2Nc, total symbol rate NcRs, noise BW B=Fs
% Es/No = (C/Rs)/(N/B), N = var = 2NcFs/NcRs(Es/No) = 2Fs/Rs(Es/No)
next_sync = sync;
[next_sync acohpsk] = coarse_freq_offset_est(acohpsk, afdmdv, ch_fdm_frame, sync, next_sync);
-
nin = M;
% shift entire FDM signal to 0 Hz
rx_filt_log = [rx_filt_log rx_filt];
[rx_onesym rx_timing env afdmdv] = rx_est_timing(afdmdv, rx_filt, nin);
+ rx_timing_log = [rx_timing_log rx_timing];
+
ch_symb(r,:) = rx_onesym;
end
rx_phi_log = [rx_phi_log; phi_];
rx_bits_log = [rx_bits_log rx_bits];
tx_bits_prev_log = [tx_bits_prev_log prev_tx_bits2];
+ ratio_log = [ratio_log acohpsk.ratio];
% BER stats
Tbits += length(error_positions);
end
- %printf("f: %d sync: %d next_sync: %d\n", i, sync, next_sync);
- sync = sync_state_machine(sync, next_sync);
+ if sync == 0
+ Nerrs = 0;
+ Tbits = 0;
+ nerr_log = [];
+ end
+
+ % printf("f: %d sync: %d next_sync: %d\n", i, sync, next_sync);
+ [sync acohpsk] = sync_state_machine(acohpsk, sync, next_sync);
prev_tx_bits2 = prev_tx_bits;
prev_tx_bits = tx_bits;
end
-stem_sig_and_error(1, 111, tx_bits_log_c(1:n), tx_bits_log(1:n) - tx_bits_log_c(1:n), 'tx bits', [1 n -1.5 1.5])
-stem_sig_and_error(2, 211, real(tx_symb_log_c(1:n)), real(tx_symb_log(1:n) - tx_symb_log_c(1:n)), 'tx symb re', [1 n -1.5 1.5])
-stem_sig_and_error(2, 212, imag(tx_symb_log_c(1:n)), imag(tx_symb_log(1:n) - tx_symb_log_c(1:n)), 'tx symb im', [1 n -1.5 1.5])
-
-stem_sig_and_error(3, 211, real(tx_fdm_frame_log_c(1:n)), real(tx_fdm_frame_log(1:n) - tx_fdm_frame_log_c(1:n)), 'tx fdm frame re', [1 n -10 10])
-stem_sig_and_error(3, 212, imag(tx_fdm_frame_log_c(1:n)), imag(tx_fdm_frame_log(1:n) - tx_fdm_frame_log_c(1:n)), 'tx fdm frame im', [1 n -10 10])
-stem_sig_and_error(4, 211, real(ch_fdm_frame_log_c(1:n)), real(ch_fdm_frame_log(1:n) - ch_fdm_frame_log_c(1:n)), 'ch fdm frame re', [1 n -10 10])
-stem_sig_and_error(4, 212, imag(ch_fdm_frame_log_c(1:n)), imag(ch_fdm_frame_log(1:n) - ch_fdm_frame_log_c(1:n)), 'ch fdm frame im', [1 n -10 10])
-stem_sig_and_error(5, 211, real(rx_fdm_frame_bb_log_c(1:n)), real(rx_fdm_frame_bb_log(1:n) - rx_fdm_frame_bb_log_c(1:n)), 'rx fdm frame bb re', [1 n -10 10])
-stem_sig_and_error(5, 212, imag(rx_fdm_frame_bb_log_c(1:n)), imag(rx_fdm_frame_bb_log(1:n) - rx_fdm_frame_bb_log_c(1:n)), 'rx fdm frame bb im', [1 n -10 10])
-
-[n m] = size(ch_symb_log);
-stem_sig_and_error(6, 211, real(ch_symb_log_c), real(ch_symb_log - ch_symb_log_c), 'ch symb re', [1 n -1.5 1.5])
-stem_sig_and_error(6, 212, imag(ch_symb_log_c), imag(ch_symb_log - ch_symb_log_c), 'ch symb im', [1 n -1.5 1.5])
-stem_sig_and_error(7, 211, real(ct_symb_ff_log_c), real(ct_symb_ff_log - ct_symb_ff_log_c), 'ct symb ff re', [1 n -1.5 1.5])
-stem_sig_and_error(7, 212, imag(ct_symb_ff_log_c), imag(ct_symb_ff_log - ct_symb_ff_log_c), 'ct symb ff im', [1 n -1.5 1.5])
-stem_sig_and_error(8, 211, rx_amp_log_c, rx_amp_log - rx_amp_log_c, 'Amp Est', [1 n -1.5 1.5])
-stem_sig_and_error(8, 212, rx_phi_log_c, rx_phi_log - rx_phi_log_c, 'Phase Est', [1 n -4 4])
-stem_sig_and_error(9, 211, real(rx_symb_log_c), real(rx_symb_log - rx_symb_log_c), 'rx symb re', [1 n -1.5 1.5])
-stem_sig_and_error(9, 212, imag(rx_symb_log_c), imag(rx_symb_log - rx_symb_log_c), 'rx symb im', [1 n -1.5 1.5])
-stem_sig_and_error(10, 111, rx_bits_log_c, rx_bits_log - rx_bits_log_c, 'rx bits', [1 n -1.5 1.5])
-
-check(tx_bits_log, tx_bits_log_c, 'tx_bits');
-check(tx_symb_log, tx_symb_log_c, 'tx_symb');
-check(tx_fdm_frame_log, tx_fdm_frame_log_c, 'tx_fdm_frame');
-check(ch_fdm_frame_log, ch_fdm_frame_log_c, 'ch_fdm_frame');
-check(rx_fdm_frame_bb_log, rx_fdm_frame_bb_log_c, 'rx_fdm_frame_bb', 0.01);
-check(ch_symb_log, ch_symb_log_c, 'ch_symb',0.01);
-check(ct_symb_ff_log, ct_symb_ff_log_c, 'ct_symb_ff',0.01);
-check(rx_amp_log, rx_amp_log_c, 'rx_amp_log',0.01);
-check(rx_phi_log, rx_phi_log_c, 'rx_phi_log',0.01);
-check(rx_symb_log, rx_symb_log_c, 'rx_symb',0.01);
-check(rx_bits_log, rx_bits_log_c, 'rx_bits');
-
-% Determine bit error rate
-
-sz = length(tx_bits_log_c);
-Nerrs_c = sum(xor(tx_bits_prev_log, rx_bits_log_c));
-Tbits_c = length(tx_bits_prev_log);
-ber_c = Nerrs_c/Tbits_c;
ber = Nerrs/Tbits;
+printf("EsNodB: %4.1f ber..: %4.3f Nerrs..: %d Tbits..: %d\n", EsNodB, ber, Nerrs, Tbits);
+
+if compare_with_c
+
+ % Output vectors from C port ---------------------------------------------------
+
+ load ../build_linux/unittest/tcohpsk_out.txt
+
+ stem_sig_and_error(1, 111, tx_bits_log_c(1:n), tx_bits_log(1:n) - tx_bits_log_c(1:n), 'tx bits', [1 n -1.5 1.5])
+ stem_sig_and_error(2, 211, real(tx_symb_log_c(1:n)), real(tx_symb_log(1:n) - tx_symb_log_c(1:n)), 'tx symb re', [1 n -1.5 1.5])
+ stem_sig_and_error(2, 212, imag(tx_symb_log_c(1:n)), imag(tx_symb_log(1:n) - tx_symb_log_c(1:n)), 'tx symb im', [1 n -1.5 1.5])
+
+ stem_sig_and_error(3, 211, real(tx_fdm_frame_log_c(1:n)), real(tx_fdm_frame_log(1:n) - tx_fdm_frame_log_c(1:n)), 'tx fdm frame re', [1 n -10 10])
+ stem_sig_and_error(3, 212, imag(tx_fdm_frame_log_c(1:n)), imag(tx_fdm_frame_log(1:n) - tx_fdm_frame_log_c(1:n)), 'tx fdm frame im', [1 n -10 10])
+ stem_sig_and_error(4, 211, real(ch_fdm_frame_log_c(1:n)), real(ch_fdm_frame_log(1:n) - ch_fdm_frame_log_c(1:n)), 'ch fdm frame re', [1 n -10 10])
+ stem_sig_and_error(4, 212, imag(ch_fdm_frame_log_c(1:n)), imag(ch_fdm_frame_log(1:n) - ch_fdm_frame_log_c(1:n)), 'ch fdm frame im', [1 n -10 10])
+ stem_sig_and_error(5, 211, real(rx_fdm_frame_bb_log_c(1:n)), real(rx_fdm_frame_bb_log(1:n) - rx_fdm_frame_bb_log_c(1:n)), 'rx fdm frame bb re', [1 n -10 10])
+ stem_sig_and_error(5, 212, imag(rx_fdm_frame_bb_log_c(1:n)), imag(rx_fdm_frame_bb_log(1:n) - rx_fdm_frame_bb_log_c(1:n)), 'rx fdm frame bb im', [1 n -10 10])
+
+ [n m] = size(ch_symb_log);
+ stem_sig_and_error(6, 211, real(ch_symb_log_c), real(ch_symb_log - ch_symb_log_c), 'ch symb re', [1 n -1.5 1.5])
+ stem_sig_and_error(6, 212, imag(ch_symb_log_c), imag(ch_symb_log - ch_symb_log_c), 'ch symb im', [1 n -1.5 1.5])
+ stem_sig_and_error(7, 211, real(ct_symb_ff_log_c), real(ct_symb_ff_log - ct_symb_ff_log_c), 'ct symb ff re', [1 n -1.5 1.5])
+ stem_sig_and_error(7, 212, imag(ct_symb_ff_log_c), imag(ct_symb_ff_log - ct_symb_ff_log_c), 'ct symb ff im', [1 n -1.5 1.5])
+ stem_sig_and_error(8, 211, rx_amp_log_c, rx_amp_log - rx_amp_log_c, 'Amp Est', [1 n -1.5 1.5])
+ stem_sig_and_error(8, 212, rx_phi_log_c, rx_phi_log - rx_phi_log_c, 'Phase Est', [1 n -4 4])
+ stem_sig_and_error(9, 211, real(rx_symb_log_c), real(rx_symb_log - rx_symb_log_c), 'rx symb re', [1 n -1.5 1.5])
+ stem_sig_and_error(9, 212, imag(rx_symb_log_c), imag(rx_symb_log - rx_symb_log_c), 'rx symb im', [1 n -1.5 1.5])
+ stem_sig_and_error(10, 111, rx_bits_log_c, rx_bits_log - rx_bits_log_c, 'rx bits', [1 n -1.5 1.5])
+
+ check(tx_bits_log, tx_bits_log_c, 'tx_bits');
+ check(tx_symb_log, tx_symb_log_c, 'tx_symb');
+ check(tx_fdm_frame_log, tx_fdm_frame_log_c, 'tx_fdm_frame');
+ check(ch_fdm_frame_log, ch_fdm_frame_log_c, 'ch_fdm_frame');
+ check(rx_fdm_frame_bb_log, rx_fdm_frame_bb_log_c, 'rx_fdm_frame_bb', 0.01);
+
+ check(ch_symb_log, ch_symb_log_c, 'ch_symb',0.01);
+ check(ct_symb_ff_log, ct_symb_ff_log_c, 'ct_symb_ff',0.01);
+ check(rx_amp_log, rx_amp_log_c, 'rx_amp_log',0.01);
+ check(rx_phi_log, rx_phi_log_c, 'rx_phi_log',0.01);
+ check(rx_symb_log, rx_symb_log_c, 'rx_symb',0.01);
+ check(rx_bits_log, rx_bits_log_c, 'rx_bits');
+
+ % Determine bit error rate
+
+ sz = length(tx_bits_log_c);
+ Nerrs_c = sum(xor(tx_bits_prev_log, rx_bits_log_c));
+ Tbits_c = length(tx_bits_prev_log);
+ ber_c = Nerrs_c/Tbits_c;
+ printf("EsNodB: %4.1f ber_c: %4.3f Nerrs_c: %d Tbits_c: %d\n", EsNodB, ber_c, Nerrs_c, Tbits_c);
+
+else
+
+ % some other useful plots
+
+ figure(1)
+ plot(rx_symb_log*exp(j*pi/4),'+')
+ figure(2)
+ plot(rx_timing_log)
+ figure(3)
+ stem(nerr_log)
+ figure(4)
+ stem(ratio_log)
-printf("EsNodB: %4.1f ber_c: %3.2f Nerrs_c: %d Tbits_c: %d\n", EsNodB, ber_c, Nerrs_c, Tbits_c);
-
-% some other useful plots
-
-if 0
-figure
-plot(rx_symb_log,'+')
-figure
-plot(xor(tx_bits_prev_log, rx_bits_log_c),'+')
end
-% C header file of noise samples so C version gives extacly the same results
+
+% function to write C header file of noise samples so C version gives
+% extactly the same results
function write_noise_file(uvnoise_log)
sim_in.ldpc_code = 0;
sim_in.Ntrials = 400;
- sim_in.Esvec = 13;
+ sim_in.Esvec = 12;
sim_in.hf_sim = 1;
sim_in.hf_mag_only = 0;
sim_in.modulation = 'qpsk';
endfunction
-% ideas: cld estimate timing with freq offset and decimate to save cpu load
-% fft to do cross correlation
-function [f_max s_max] = test_freq_off_est(rx_filename, offset, n)
- fpilot = fopen("tx_zero.raw","rb"); tx_pilot = fread(fpilot, "short"); fclose(fpilot);
- frx=fopen(rx_filename,"rb"); rx_fdm = fread(frx, "short"); fclose(frx);
- Fs = 8000;
- nc = 1800; % portion we wish to correlate over (first 2 rows on pilots)
-
- % downconvert to complex baseband to remove images
-
- f = 1000;
- foff_rect = exp(j*2*pi*f*(1:2*n)/Fs);
- tx_pilot_bb = tx_pilot(1:n) .* foff_rect(1:n)';
- rx_fdm_bb = rx_fdm(offset:offset+2*n-1) .* foff_rect';
-
- % remove -2000 Hz image
-
- b = fir1(50, 1000/Fs);
- tx_pilot_bb_lpf = filter(b,1,tx_pilot_bb);
- rx_fdm_bb_lpf = filter(b,1,rx_fdm_bb);
-
- % decimate by M
-
- M = 4;
- tx_pilot_bb_lpf = tx_pilot_bb_lpf(1:M:length(tx_pilot_bb_lpf));
- rx_fdm_bb_lpf = rx_fdm_bb_lpf(1:M:length(rx_fdm_bb_lpf));
- n /= M;
- nc /= M;
-
- % correlate over a range of frequency offsets and delays
-
- c_max = 0;
- f_n = 1;
- f_range = -75:2.5:75;
- c_log=zeros(n, length(f_range));
-
- for f=f_range
- foff_rect = exp(j*2*pi*(f*M)*(1:nc)/Fs);
- for s=1:n
-
- c = abs(tx_pilot_bb_lpf(1:nc)' * (rx_fdm_bb_lpf(s:s+nc-1) .* foff_rect'));
- c_log(s,f_n) = c;
- if c > c_max
- c_max = c;
- f_max = f;
- s_max = s;
- end
- end
- f_n++;
- %printf("f: %f c_max: %f f_max: %f s_max: %d\n", f, c_max, f_max, s_max);
- end
-
- figure(1);
- y = f_range;
- x = max(s_max-25,1):min(s_max+25, n);
- mesh(y,x, c_log(x,:));
- grid
-
- s_max *= M;
- s_max -= floor(s_max/6400)*6400;
- printf("f_max: %f s_max: %d\n", f_max, s_max);
-
- % decimated position at sample rate. need to relate this to symbol
- % rate position.
-
-endfunction
function snr = EsNo_to_SNR(EsNo)
snr = interp1([20 11.8 9.0 6.7 4.9 3.3 -10], [ 3 3 0 -3 -6 -9 -9], EsNo);
endfunction
+
function gen_test_bits()
sim_in = standard_init();
framesize = 32*10;
test_bits_coh_file(tx_bits);
endfunction
+
+
+
% Start simulations ---------------------------------------
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")
%rate_Fs_rx("tx.raw")
-%test_freq_off_est("tx.raw",40,6400)
%gen_test_bits();
--- /dev/null
+% test_foff.m
+% David Rowe April 2015
+%
+% Octave script for testing the cohpsk freq offset estimator
+
+graphics_toolkit ("gnuplot");
+more off;
+
+cohpsk;
+fdmdv;
+autotest;
+
+rand('state',1);
+randn('state',1);
+
+% Core function for testing frequency offset estimator. Performs one test
+
+function sim_out = freq_off_est_test(sim_in)
+ global Nfilter;
+ global M;
+
+ Rs = 50;
+ Nc = 4;
+ Nd = 2;
+ framesize = 32;
+ Fs = 8000;
+ Fcentre = 1500;
+
+ afdmdv.Nsym = 2;
+ afdmdv.Nt = 3;
+
+ afdmdv.Fs = 8000;
+ afdmdv.Nc = Nd*Nc-1;
+ afdmdv.Rs = Rs;
+ afdmdv.M = afdmdv.Fs/afdmdv.Rs;
+ afdmdv.tx_filter_memory = zeros(afdmdv.Nc+1, Nfilter);
+ afdmdv.Nfilter = Nfilter;
+ afdmdv.gt_alpha5_root = gen_rn_coeffs(0.5, 1/Fs, Rs, afdmdv.Nsym, afdmdv.M);
+ afdmdv.Fsep = 75;
+ afdmdv.phase_tx = ones(afdmdv.Nc+1,1);
+ freq_hz = afdmdv.Fsep*( -Nc*Nd/2 - 0.5 + (1:Nc*Nd) );
+ afdmdv.freq_pol = 2*pi*freq_hz/Fs;
+ afdmdv.freq = exp(j*afdmdv.freq_pol);
+ afdmdv.Fcentre = 1500;
+
+ afdmdv.fbb_rect = exp(j*2*pi*Fcentre/Fs);
+ afdmdv.fbb_phase_tx = 1;
+ afdmdv.fbb_phase_rx = 1;
+ afdmdv.phase_rx = ones(afdmdv.Nc+1,1);
+
+ nin = M;
+
+ P = afdmdv.P = 4;
+ afdmdv.Nfilter = afdmdv.Nsym*afdmdv.M;
+ afdmdv.rx_filter_mem_timing = zeros(afdmdv.Nc+1, afdmdv.Nt*afdmdv.P);
+ afdmdv.Nfiltertiming = afdmdv.M + afdmdv.Nfilter + afdmdv.M;
+ afdmdv.rx_filter_memory = zeros(afdmdv.Nc+1, afdmdv.Nfilter);
+
+ acohpsk = standard_init();
+ acohpsk.framesize = framesize;
+ acohpsk.ldpc_code = 0;
+ acohpsk.ldpc_code_rate = 1;
+ acohpsk.Nc = Nc;
+ acohpsk.Rs = Rs;
+ acohpsk.Ns = 4;
+ acohpsk.Nd = Nd;
+ acohpsk.modulation = 'qpsk';
+ acohpsk.do_write_pilot_file = 0;
+ acohpsk = symbol_rate_init(acohpsk);
+ acohpsk.Ncm = 10*acohpsk.Nsymbrowpilot*M;
+ acohpsk.coarse_mem = zeros(1,acohpsk.Ncm);
+ acohpsk.Ndft = 2^(ceil(log2(acohpsk.Ncm)));
+
+ frames = sim_in.frames;
+ EsNodB = sim_in.EsNodB;
+ foff = sim_in.foff;
+ dfoff = sim_in.dfoff;
+ fading_en = sim_in.fading_en;
+
+ EsNo = 10^(EsNodB/10);
+ hf_delay_ms = 2;
+ phase_ch = 1;
+
+ rand('state',1);
+ tx_bits_coh = round(rand(1,framesize*10));
+ ptx_bits_coh = 1;
+ [spread spread_2ms hf_gain] = init_hf_model(Fs, 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);
+
+ sync = 0; next_sync = 1;
+ sync_start = 1;
+ freq_offset_log = [];
+ sync_time_log = [];
+ ch_fdm_frame_log = [];
+ ch_symb_log = [];
+ tx_fdm_frame_log = [];
+
+ for f=1:frames
+
+ acohpsk.frame = f;
+
+ %
+ % Mod --------------------------------------------------------------------
+ %
+
+ tx_bits = tx_bits_coh(ptx_bits_coh:ptx_bits_coh+framesize-1);
+ ptx_bits_coh += framesize;
+ if ptx_bits_coh > length(tx_bits_coh)
+ ptx_bits_coh = 1;
+ end
+
+ [tx_symb tx_bits] = bits_to_qpsk_symbols(acohpsk, tx_bits, [], []);
+
+ tx_fdm_frame = [];
+ for r=1:acohpsk.Nsymbrowpilot
+ tx_onesymb = tx_symb(r,:);
+ [tx_baseband afdmdv] = tx_filter(afdmdv, tx_onesymb);
+ [tx_fdm afdmdv] = fdm_upconvert(afdmdv, tx_baseband);
+ tx_fdm_frame = [tx_fdm_frame tx_fdm];
+ end
+ tx_fdm_frame_log = [tx_fdm_frame_log tx_fdm_frame];
+
+ %
+ % Channel --------------------------------------------------------------------
+ %
+
+ ch_fdm_frame = zeros(1,acohpsk.Nsymbrowpilot*M);
+ for i=1:acohpsk.Nsymbrowpilot*M
+ foff_rect = exp(j*2*pi*foff/Fs);
+ foff += dfoff;
+ phase_ch *= foff_rect;
+ ch_fdm_frame(i) = tx_fdm_frame(i) * phase_ch;
+ end
+ phase_ch /= abs(phase_ch);
+
+ if fading_en
+ ch_fdm_delay(1:nhfdelay) = ch_fdm_delay(acohpsk.Nsymbrowpilot*M+1:nhfdelay+acohpsk.Nsymbrowpilot*M);
+ ch_fdm_delay(nhfdelay+1:nhfdelay+acohpsk.Nsymbrowpilot*M) = ch_fdm_frame;
+
+ for i=1:acohpsk.Nsymbrowpilot*M
+ ahf_model = hf_gain*(spread(hf_n)*ch_fdm_frame(i) + spread_2ms(hf_n)*ch_fdm_delay(i));
+ ch_fdm_frame(i) = ahf_model;
+ hf_n++;
+ end
+ end
+
+ % each carrier has power = 2, total power 2Nc, total symbol rate NcRs, noise BW B=Fs
+ % Es/No = (C/Rs)/(N/B), N = var = 2NcFs/NcRs(Es/No) = 2Fs/Rs(Es/No)
+
+ variance = 2*Fs/(acohpsk.Rs*EsNo);
+ uvnoise = sqrt(0.5)*(randn(1,acohpsk.Nsymbrowpilot*M) + j*randn(1,acohpsk.Nsymbrowpilot*M));
+ noise = sqrt(variance)*uvnoise;
+
+ ch_fdm_frame += noise;
+ ch_fdm_frame_log = [ch_fdm_frame_log ch_fdm_frame];
+
+ %
+ % Try to achieve sync --------------------------------------------------------------------
+ %
+
+ next_sync = sync;
+
+ if sync == 0
+ next_sync = 2;
+ acohpsk.f_est = Fcentre;
+ end
+
+ [rx_fdm_frame_bb afdmdv.fbb_phase_rx] = freq_shift(ch_fdm_frame, -acohpsk.f_est, Fs, afdmdv.fbb_phase_rx);
+
+ for r=1:acohpsk.Nsymbrowpilot
+
+ % downconvert each FDM carrier to Nc separate baseband signals
+
+ [rx_baseband afdmdv] = fdm_downconvert(afdmdv, rx_fdm_frame_bb(1+(r-1)*M:r*M), nin);
+ [rx_filt afdmdv] = rx_filter(afdmdv, rx_baseband, nin);
+ [rx_onesym rx_timing env afdmdv] = rx_est_timing(afdmdv, rx_filt, nin);
+
+ ch_symb(r,:) = rx_onesym;
+ end
+ ch_symb_log = [ch_symb_log; ch_symb];
+
+ % coarse timing (frame sync) and initial fine freq est ---------------------------------------------
+
+ [next_sync acohpsk] = frame_sync_fine_freq_est(acohpsk, ch_symb, sync, next_sync);
+
+ % if we've acheived sync gather stats
+
+ if (next_sync == 4)
+ freq_offset_log = [freq_offset_log acohpsk.f_fine_est+foff];
+ sync_time_log = [sync_time_log f-sync_start];
+ sync = 0; next_sync = 2; sync_start = f;
+ end
+
+ %printf("f: %d sync: %d next_sync: %d\n", f, sync, next_sync);
+ [sync acohpsk] = sync_state_machine(acohpsk, sync, next_sync);
+
+ end
+
+ % ftx=fopen("coarse_tx.raw","wb"); fwrite(ftx, 1000*ch_fdm_frame_log, "short"); fclose(ftx);
+
+ sim_out.freq_offset_log = freq_offset_log;
+ sim_out.sync_time_log = sync_time_log;
+ sim_out.ch_fdm_frame_log = ch_fdm_frame_log;
+ sim_out.ch_symb_log = ch_symb_log;
+ sim_out.tx_fdm_frame_log = tx_fdm_frame_log;
+endfunction
+
+
+function freq_off_est_test_single
+ sim_in.frames = 100;
+ sim_in.EsNodB = 20;
+ sim_in.foff = -15;
+ sim_in.dfoff = 0;
+ sim_in.fading_en = 1;
+
+ sim_out = freq_off_est_test(sim_in);
+
+ figure(1)
+ subplot(211)
+ plot(sim_out.freq_offset_log)
+ subplot(212)
+ hist(sim_out.freq_offset_log)
+
+ figure(2)
+ subplot(211)
+ plot(sim_out.sync_time_log)
+ subplot(212)
+ hist(sim_out.sync_time_log)
+
+ figure(3)
+ subplot(211)
+ plot(real(sim_out.tx_fdm_frame_log(1:2*960)))
+ subplot(212)
+ plot(real(sim_out.ch_symb_log(1:24,:)),'+')
+endfunction
+
+
+function [freq_off_log EsNodBSet] = freq_off_est_test_curves
+ EsNodBSet = [20 12 8];
+
+ sim_in.frames = 100;
+ sim_in.foff = -20;
+ sim_in.dfoff = 0;
+ sim_in.fading_en = 1;
+ freq_off_log = 1E6*ones(sim_in.frames, length(EsNodBSet) );
+ sync_time_log = 1E6*ones(sim_in.frames, length(EsNodBSet) );
+
+ for i=1:length(EsNodBSet)
+ sim_in.EsNodB = EsNodBSet(i);
+ printf("%f\n", sim_in.EsNodB);
+
+ sim_out = freq_off_est_test(sim_in);
+ freq_off_log(1:length(sim_out.freq_offset_log),i) = sim_out.freq_offset_log;
+ sync_time_log(1:length(sim_out.sync_time_log),i) = sim_out.sync_time_log;
+ end
+
+ figure(1)
+ clf
+ hold on;
+ for i=1:length(EsNodBSet)
+ data = freq_off_log(find(freq_off_log(:,i) < 1E6),i);
+ s = std(data);
+ m = mean(data);
+ stdbar = [m-s; m+s];
+ plot(EsNodBSet(i), data, '+')
+ plot([EsNodBSet(i) EsNodBSet(i)]+0.5, stdbar,'+-')
+ end
+ hold off
+
+ axis([6 22 -25 25])
+ if sim_in.fading_en
+ title_str = sprintf('foff = %d Hz Fading', sim_in.foff);
+ else
+ title_str = sprintf('foff = %d Hz AWGN', sim_in.foff);
+ end
+ title(title_str);
+ xlabel('Es/No (dB)')
+ ylabel('freq offset error (Hz)');
+
+ figure(2)
+ clf
+ hold on;
+ for i=1:length(EsNodBSet)
+ leg = sprintf("%d;%d dB;", i, EsNodBSet(i));
+ plot(freq_off_log(find(freq_off_log(:,i) < 1E6),i),leg)
+ end
+ hold off;
+ title(title_str);
+ xlabel('test');
+ ylabel('freq offset error (Hz)');
+ legend('boxoff');
+
+ figure(3)
+ clf
+ hold on;
+ for i=1:length(EsNodBSet)
+ data = sync_time_log(find(sync_time_log(:,i) < 1E6),i);
+ if length(data)
+ s = std(data);
+ m = mean(data);
+ stdbar = [m-s; m+s];
+ plot(EsNodBSet(i), data, '+')
+ plot([EsNodBSet(i) EsNodBSet(i)]+0.5, stdbar,'+-')
+ end
+ end
+ hold off;
+ axis([6 22 0 10])
+ ylabel('sync time (frames)')
+ xlabel('Es/No (dB)');
+ title(title_str);
+
+ figure(4)
+ clf
+ hold on;
+ for i=1:length(EsNodBSet)
+ leg = sprintf("%d;%d dB;", i, EsNodBSet(i));
+ plot(sync_time_log(find(sync_time_log(:,i) < 1E6),i),leg)
+ end
+ hold off;
+ title(title_str);
+ xlabel('test');
+ ylabel('sync time (frames)');
+ legend('boxoff');
+
+endfunction
+
+
+freq_off_est_test_single;
+%freq_off_est_test_curves;
+
+% 1. start with +/- 20Hz offset
+% 2. Measure frames to sync. How to define sync? Foff to withn 1 Hz. Sync state
+% Need to see if we get false sync
+% 3. Try shortened filter
+% 4. Extend to parallel demods at +/-
projects / freetel-svn-tracking.git / commitdiff