From f48021b1657848f1fb8b2909bd666f9dd5c25ed7 Mon Sep 17 00:00:00 2001 From: drowe67 Date: Tue, 17 Jun 2014 03:11:42 +0000 Subject: [PATCH] first pass at freq offset est unit test git-svn-id: https://svn.code.sf.net/p/freetel/code@1662 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/fdmdv_demod.m | 1 + codec2-dev/octave/fdmdv_demod_c.m | 2 +- codec2-dev/octave/fdmdv_ut_freq_off.m | 157 ++++++++++++++++++++++++++ codec2-dev/octave/test_dqpsk2.m | 6 +- codec2-dev/octave/test_qpsk.m | 14 +-- codec2-dev/octave/test_qpsk3.m | 10 +- 6 files changed, 176 insertions(+), 14 deletions(-) create mode 100644 codec2-dev/octave/fdmdv_ut_freq_off.m diff --git a/codec2-dev/octave/fdmdv_demod.m b/codec2-dev/octave/fdmdv_demod.m index 2e403cc6..53a303ac 100644 --- a/codec2-dev/octave/fdmdv_demod.m +++ b/codec2-dev/octave/fdmdv_demod.m @@ -207,6 +207,7 @@ function fdmdv_demod(rawfilename, nbits, NumCarriers, errorpatternfilename, symb end end end + test_frame_sync_state = next_test_frame_sync_state; test_frame_sync_log = [test_frame_sync_log test_frame_sync_state]; end diff --git a/codec2-dev/octave/fdmdv_demod_c.m b/codec2-dev/octave/fdmdv_demod_c.m index c2813e62..d63f8510 100644 --- a/codec2-dev/octave/fdmdv_demod_c.m +++ b/codec2-dev/octave/fdmdv_demod_c.m @@ -10,7 +10,7 @@ % function fdmdv_demod_c(dumpfilename, bits) - + NumCarriers = 16; fdmdv; % include modem code frames = bits/(Nc*Nb); diff --git a/codec2-dev/octave/fdmdv_ut_freq_off.m b/codec2-dev/octave/fdmdv_ut_freq_off.m new file mode 100644 index 00000000..bfc60f82 --- /dev/null +++ b/codec2-dev/octave/fdmdv_ut_freq_off.m @@ -0,0 +1,157 @@ +% fdmdv_ut_freq_off.m +% David Rowe 17 June 2014 +% +% Unit Test program for freq offset estimation in FDMDV modem. +% +% Copyright David Rowe 2012 This program is +% distributed under the terms of the GNU General Public License +% Version 2 + +% [ ] sweep of different delays +% [ ] sweep of Eb/No +% [ ] sweep of freq offsets +% [ ] plot/print pass fail/relevant stats +% + variance +% + histogram of freq ests? + +fdmdv; % load modem code + +% Simulation Parameters -------------------------------------- + +frames = 100; +EbNovec = 10; +Foff_hz = 0; + +% --------------------------------------------------------------------- +% Eb/No calculations. We need to work out Eb/No for each FDM carrier. +% Total power is sum of power in all FDM carriers +% --------------------------------------------------------------------- + +function Nsd = calc_Nsd_from_EbNo(EbNo_dB) + global Rs; + global Nb; + global Fs; + + C = 1; % power of each FDM carrier (energy/sample). Total Carrier power should = Nc*C = Nc + N = 1; % total noise power (energy/sample) of noise source across entire bandwidth + + % Eb = Carrier power * symbol time / (bits/symbol) + % = C *(1/Rs) / Nb + Eb_dB = 10*log10(C) - 10*log10(Rs) - 10*log10(Nb); + + No_dBHz = Eb_dB - EbNo_dB; + + % Noise power = Noise spectral density * bandwidth + % Noise power = Noise spectral density * Fs/2 for real signals + N_dB = No_dBHz + 10*log10(Fs/2); + Ngain_dB = N_dB - 10*log10(N); + Nsd = 10^(Ngain_dB/20); +end + +% ------------------------------------------------------------ + +modulation = 'dqpsk'; +tx_filt = zeros(Nc,M); +tx_fdm_log = []; +rx_fdm_log = []; +prev_tx_symbols = ones(Nc+1,1); +ferr = 0; +foff = 0; +foff_log = []; + +% fixed delay simuation + +Ndelay = M+20; +rx_fdm_delay = zeros(Ndelay,1); + +% --------------------------------------------------------------------- +% Eb/No calculations. We need to work out Eb/No for each FDM carrier. +% Total power is sum of power in all FDM carriers +% --------------------------------------------------------------------- + +% freq offset simulation states + +phase_offset = 1; +freq_offset = exp(j*2*pi*Foff_hz/Fs); +foff_phase = 1; +t = 0; +foff = 0; + +% --------------------------------------------------------------------- +% Main loop +% --------------------------------------------------------------------- + +for ne = 1:length(EbNovec) + EbNo_dB = EbNovec(ne); + Nsd = calc_Nsd_from_EbNo(EbNo_dB); + + for f=1:frames + + % ------------------- Modulator ------------------- + + tx_bits = get_test_bits(Nc*Nb); + tx_symbols = bits_to_psk(prev_tx_symbols, tx_bits, modulation); + prev_tx_symbols = tx_symbols; + tx_baseband = tx_filter(tx_symbols); + tx_fdm = fdm_upconvert(tx_baseband); + tx_fdm_log = [tx_fdm_log real(tx_fdm)]; + + % ------------------- Channel simulation ------------------- + + % frequency offset + + Foff = Foff_hz; for i=1:M freq_offset = exp(j*2*pi*Foff/Fs); + phase_offset *= freq_offset; tx_fdm(i) = phase_offset*tx_fdm(i); end + + rx_fdm = real(tx_fdm); + + % AWGN noise + + noise = Nsd*randn(1,M); + rx_fdm += noise; + rx_fdm_log = [rx_fdm_log rx_fdm]; + + % Delay + + rx_fdm_delay(1:Ndelay-M) = rx_fdm_delay(M+1:Ndelay); + rx_fdm_delay(Ndelay-M+1:Ndelay) = rx_fdm; + %rx_fdm_delay = rx_fdm; + + % ------------------- Freq Offset Est ------------------- + + % frequency offset estimation and correction, need to call + % rx_est_freq_offset even in track mode to keep states updated + + [pilot prev_pilot pilot_lut_index prev_pilot_lut_index] = ... + get_pilot(pilot_lut_index, prev_pilot_lut_index, M); + [foff_coarse S1 S2] = rx_est_freq_offset(rx_fdm_delay, pilot, prev_pilot, M); + foff_log(ne,f) = foff_coarse; + end +end + +% --------------------------------------------------------------------- +% Print Stats +% --------------------------------------------------------------------- + +% --------------------------------------------------------------------- +% Plots +% --------------------------------------------------------------------- + +figure(1) +clf +for ne = 1:length(EbNovec) + foff_std(ne) = std(foff_log(ne,:)); +end +plot(EbNovec,foff_std) +xlabel("Eb/No (dB)") +ylabel("Std Dev") + +figure(2) +clf; +plot(foff_log(1,:)) +xlabel("Frames") +ylabel("Freq offset estimate") + +figure(3) +clf; +hist(foff_log(1,:)) diff --git a/codec2-dev/octave/test_dqpsk2.m b/codec2-dev/octave/test_dqpsk2.m index 01f45677..08787510 100644 --- a/codec2-dev/octave/test_dqpsk2.m +++ b/codec2-dev/octave/test_dqpsk2.m @@ -407,8 +407,10 @@ pdB = pdB(4:4:length(pdB)); % Use linear mapping function in dB domain to map to symbol power -power_map_x = [ 0 20 24 40 50 ]; -power_map_y = [-6 -6 0 6 6]; +%power_map_x = [ 0 20 24 40 50 ]; +%power_map_y = [--6 -6 0 6 6]; +power_map_x = [ 0 50 ]; +power_map_y = [ -15 12]; mapped_pdB = interp1(power_map_x, power_map_y, pdB); sim_in.symbol_amp = 10 .^ (mapped_pdB/20); diff --git a/codec2-dev/octave/test_qpsk.m b/codec2-dev/octave/test_qpsk.m index 987841cf..6d8796e0 100644 --- a/codec2-dev/octave/test_qpsk.m +++ b/codec2-dev/octave/test_qpsk.m @@ -6,6 +6,7 @@ % Design to test coherent demodulation ideas on HF channels without % building a full blown modem. Uses 'genie provided' estimates for % timing estimation, frame sync. +% 1; @@ -219,7 +220,6 @@ function sim_out = ber_test(sim_in, modulation) comb = conj(spread(sc:sc+M-1))' + conj(spread_2ms(sc:sc+M-1))'; if hf_phase_only tx_filt = tx_filt.*exp(j*angle(comb)); - hf_angle_log = [hf_angle_log angle(comb)]; else if hf_mag_only comb = conj(spread(sc:sc+M-1))' + conj(spread_2ms(sc:sc+M-1))'; @@ -239,6 +239,7 @@ function sim_out = ber_test(sim_in, modulation) C_log = [C_log abs(comb)*hf_gain]; end tx_filt_log = [tx_filt_log tx_filt]; + hf_angle_log = [hf_angle_log angle(comb)]; % AWGN noise and phase/freq offset channel simulation % 0.5 factor ensures var(noise) == variance , i.e. splits power between Re & Im @@ -413,12 +414,11 @@ function sim_out = ber_test(sim_in, modulation) figure(3); clf; - if hf_phase_only - plot(hf_angle_log); - axis([1 10000 min(hf_angle_log) max(hf_angle_log)]) - else - plot(C_log); - end + subplot(211) + plot(C_log); + subplot(212) + plot(hf_angle_log); + axis([1 10000 min(hf_angle_log) max(hf_angle_log)]) end endfunction diff --git a/codec2-dev/octave/test_qpsk3.m b/codec2-dev/octave/test_qpsk3.m index 37efad87..eaf69995 100644 --- a/codec2-dev/octave/test_qpsk3.m +++ b/codec2-dev/octave/test_qpsk3.m @@ -573,12 +573,14 @@ function ideal sim_in = standard_init(); - sim_in.verbose = 1; + sim_in.sim_coh_dpsk = 0; + sim_in.newldpc = 1; + sim_in.verbose = 2; sim_in.plot_scatter = 1; sim_in.Esvec = 5; sim_in.hf_sim = 1; - sim_in.Ntrials = 100; + sim_in.Ntrials = 30; sim_qpsk_hf = ber_test(sim_in, 'qpsk'); @@ -948,8 +950,8 @@ endfunction more off; -%ideal(); +ideal(); %phase_est_hf(); %phase_est_awgn(); %test_dpsk(); -gen_error_pattern_qpsk +%gen_error_pattern_qpsk -- 2.25.1