From 51546ddbe8255d79dae981150f8cf9ce68483a4e Mon Sep 17 00:00:00 2001 From: drowe67 Date: Fri, 13 Apr 2012 01:20:40 +0000 Subject: [PATCH] working two state freq est/tracking, sample clock correction logic, tested over pathsim with a few models with reas results git-svn-id: https://svn.code.sf.net/p/freetel/code@366 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/README_fdmdv.txt | 1 + codec2-dev/octave/fdmdv.m | 160 +++++++++++++++++++++-------- codec2-dev/octave/fdmdv_demod.m | 130 ++++++++++++----------- codec2-dev/octave/fdmdv_ut.m | 39 ++++--- 4 files changed, 207 insertions(+), 123 deletions(-) diff --git a/codec2-dev/octave/README_fdmdv.txt b/codec2-dev/octave/README_fdmdv.txt index 57d29c5e..a22b65a5 100644 --- a/codec2-dev/octave/README_fdmdv.txt +++ b/codec2-dev/octave/README_fdmdv.txt @@ -18,6 +18,7 @@ TODO + does modem fall over? [ ] try non-flat channel, e.g. 3dB difference between hi and low tones + make sure all estimators keep working +[ ] test rx level sensitivity, i.e. 0 to 20dB attenuation [ ] try to run from shell script [ ] run a few tests [ ] start coding in C and repeat tests diff --git a/codec2-dev/octave/fdmdv.m b/codec2-dev/octave/fdmdv.m index 444fef8e..7f65353c 100644 --- a/codec2-dev/octave/fdmdv.m +++ b/codec2-dev/octave/fdmdv.m @@ -82,8 +82,6 @@ function tx_symbols = bits_to_qpsk(prev_tx_symbols, tx_bits, modulation) if (strcmp(modulation,'dqpsk')) - - if 1 % map to (Nc,1) DQPSK symbols for c=1:Nc @@ -102,21 +100,6 @@ function tx_symbols = bits_to_qpsk(prev_tx_symbols, tx_bits, modulation) tx_symbols(c) = -j*prev_tx_symbols(c); endif end - else - % map to pi/4 DQPSK from spra341 Eq (6) & (7) - - for c=1:Nc - - msb = tx_bits_matrix(c,1); lsb = tx_bits_matrix(c,2); - a = 2*msb - 1; - b = 2*lsb - 1; - p = prev_tx_symbols(c); - inphase = (real(p)*a - imag(p)*b)*0.707; - quadrature = (imag(p)*a + real(p)*b)*0.707; - tx_symbols(c) = inphase + j*quadrature; - end - end - else % QPSK mapping tx_symbols = -1 + 2*tx_bits_matrix(:,1) - j + 2j*tx_bits_matrix(:,2); @@ -215,7 +198,7 @@ endfunction % Frequency shift each modem carrier down to Nc baseband signals -function rx_baseband = fdm_downconvert(rx_fdm) +function rx_baseband = fdm_downconvert(rx_fdm, nin) global Fs; global M; global Nc; @@ -223,12 +206,12 @@ function rx_baseband = fdm_downconvert(rx_fdm) global phase_rx; global freq; - rx_baseband = zeros(1,M); + rx_baseband = zeros(1,nin); % Nc/2 tones below centre freq for c=1:Nc/2 - for i=1:M + for i=1:nin phase_rx(c) = phase_rx(c) * freq(c); rx_baseband(c,i) = rx_fdm(i)*phase_rx(c)'; end @@ -237,7 +220,7 @@ function rx_baseband = fdm_downconvert(rx_fdm) % Nc/2 tones above centre freq for c=Nc/2+1:Nc - for i=1:M + for i=1:nin phase_rx(c) = phase_rx(c) * freq(c); rx_baseband(c,i) = rx_fdm(i)*phase_rx(c)'; end @@ -246,7 +229,7 @@ function rx_baseband = fdm_downconvert(rx_fdm) % Pilot c = Nc+1; - for i=1:M + for i=1:nin phase_rx(c) = phase_rx(c) * freq(c); rx_baseband(c,i) = rx_fdm(i)*phase_rx(c)'; end @@ -256,7 +239,7 @@ endfunction % Receive filter each baseband signal at oversample rate P -function rx_filt = rx_filter(rx_baseband) +function rx_filt = rx_filter(rx_baseband, nin) global Nc; global M; global P; @@ -265,14 +248,14 @@ function rx_filt = rx_filter(rx_baseband) global gt_alpha5_root; global Fsep; - rx_filt = zeros(Nc+1,P); + rx_filt = zeros(Nc+1,nin*P/M); % rx filter each symbol, generate P filtered output samples for each symbol. % Note we keep memory at rate M, it's just the filter output at rate P N=M/P; j=1; - for i=1:N:M + for i=1:N:nin rx_filter_memory(:,Nfilter-N+1:Nfilter) = rx_baseband(:,i:i-1+N); rx_filt(:,j) = rx_filter_memory * gt_alpha5_root'; rx_filter_memory(:,1:Nfilter-N) = rx_filter_memory(:,1+N:Nfilter); @@ -284,7 +267,7 @@ endfunction % Estimate frequency offset of FDM signal using BPSK pilot. This is quite % sensitive to pilot tone level wrt other carriers -function foff = rx_est_freq_offset(rx_fdm, pilot, pilot_prev) +function foff = rx_est_freq_offset(rx_fdm, pilot, pilot_prev, nin) global M; global Npilotbaseband; global pilot_baseband1; @@ -295,15 +278,15 @@ function foff = rx_est_freq_offset(rx_fdm, pilot, pilot_prev) % down convert latest M samples of pilot by multiplying by % ideal BPSK pilot signal we have generated locally - pilot_baseband1(1:Npilotbaseband-M) = pilot_baseband1(M+1:Npilotbaseband); - pilot_baseband2(1:Npilotbaseband-M) = pilot_baseband2(M+1:Npilotbaseband); - for i=1:M - pilot_baseband1(Npilotbaseband-M+i) = rx_fdm(i) * conj(pilot(i)); - pilot_baseband2(Npilotbaseband-M+i) = rx_fdm(i) * conj(pilot_prev(i)); + pilot_baseband1(1:Npilotbaseband-nin) = pilot_baseband1(nin+1:Npilotbaseband); + pilot_baseband2(1:Npilotbaseband-nin) = pilot_baseband2(nin+1:Npilotbaseband); + for i=1:nin + pilot_baseband1(Npilotbaseband-nin+i) = rx_fdm(i) * conj(pilot(i)); + pilot_baseband2(Npilotbaseband-nin+i) = rx_fdm(i) * conj(pilot_prev(i)); end - [foff1 max1 pilot_lpf1] = lpf_peak_pick(pilot_baseband1, pilot_lpf1); - [foff2 max2 pilot_lpf2] = lpf_peak_pick(pilot_baseband2, pilot_lpf2); + [foff1 max1 pilot_lpf1] = lpf_peak_pick(pilot_baseband1, pilot_lpf1, nin); + [foff2 max2 pilot_lpf2] = lpf_peak_pick(pilot_baseband2, pilot_lpf2, nin); if max1 > max2 foff = foff1; @@ -315,7 +298,7 @@ endfunction % LPF and peak pick part of freq est, put in a function as we call it twice -function [foff imax pilot_lpf] = lpf_peak_pick(pilot_baseband, pilot_lpf) +function [foff imax pilot_lpf] = lpf_peak_pick(pilot_baseband, pilot_lpf, nin) global M; global Npilotlpf; global Npilotcoeff; @@ -325,9 +308,9 @@ function [foff imax pilot_lpf] = lpf_peak_pick(pilot_baseband, pilot_lpf) % LPF cutoff 200Hz, so we can handle max +/- 200 Hz freq offset - pilot_lpf(1:Npilotlpf-M) = pilot_lpf(M+1:Npilotlpf); + pilot_lpf(1:Npilotlpf-nin) = pilot_lpf(nin+1:Npilotlpf); j = 1; - for i = Npilotlpf-M+1:Npilotlpf + for i = Npilotlpf-nin+1:Npilotlpf pilot_lpf(i) = pilot_baseband(j:j+Npilotcoeff) * pilot_coeff'; j++; end @@ -354,7 +337,7 @@ endfunction % Estimate optimum timing offset, and symbol receive symbols -function [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband) +function [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband, nin) global M; global Nt; global Nc; @@ -365,10 +348,22 @@ function [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband) global Nfiltertiming; global gt_alpha5_root; + % nin adjust + % -------------------------------- + % 120 -1 (one less rate P sample) + % 160 0 (nominal) + % 200 1 (one more rate P sample) + + adjust = P - nin*P/M; + % update buffer of Nt rate P filtered symbols - rx_filter_mem_timing(:,1:(Nt-1)*P) = rx_filter_mem_timing(:,P+1:Nt*P); - rx_filter_mem_timing(:,(Nt-1)*P+1:Nt*P) = rx_filt(1:Nc,:); + rx_filter_mem_timing(:,1:(Nt-1)*P+adjust) = rx_filter_mem_timing(:,P+1-adjust:Nt*P); + %size((Nt-1)*P+1+adjust:Nt*P) + %size(rx_filt) + %adjust + %nin + rx_filter_mem_timing(:,(Nt-1)*P+1+adjust:Nt*P) = rx_filt(1:Nc,:); % sum envelopes of all carriers @@ -396,8 +391,8 @@ function [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband) % baseband signal at rate M this enables us to resample the filtered % rx symbol with M sample precision once we have rx_timing - rx_baseband_mem_timing(:,1:Nfiltertiming-M) = rx_baseband_mem_timing(:,M+1:Nfiltertiming); - rx_baseband_mem_timing(:,Nfiltertiming-M+1:Nfiltertiming) = rx_baseband; + rx_baseband_mem_timing(:,1:Nfiltertiming-nin) = rx_baseband_mem_timing(:,nin+1:Nfiltertiming); + rx_baseband_mem_timing(:,Nfiltertiming-nin+1:Nfiltertiming) = rx_baseband; % sample right in the middle of the timing estimator window, by filtering % at rate M @@ -426,7 +421,7 @@ endfunction % convert symbols back to an array of bits -function [rx_bits sync_bit] = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation) +function [rx_bits sync_bit f_err] = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation) global Nc; global Nb; global Nb; @@ -462,8 +457,10 @@ function [rx_bits sync_bit] = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulati phase_difference(Nc+1) = rx_symbols(Nc+1) .* conj(prev_rx_symbols(Nc+1)); if (real(phase_difference(Nc+1)) < 0) sync_bit = 0; + f_err = imag(phase_difference(Nc+1)); else sync_bit = 1; + f_err = -imag(phase_difference(Nc+1)); end else @@ -709,3 +706,82 @@ function [buf_out t nin] = resample(buf_in, t, ratio, nout) t -= delta; endfunction + + +% freq offset state machine. Moves between acquire and track states based +% on BPSK pilot sequence. Freq offset estimator occasionally makes mistakes +% when used continuously. So we use it unit we have acquired the BPSK pilot, +% then switch to a more robust tracking algorithm. If we lose sync we switch +% back to acquire mode for fast-requisition. + +function [track state] = freq_state(sync_bit, state) + + % acquire state, look for 6 symbol 010101 sequence from sync bit + + next_state = state; + if state == 0 + if sync_bit == 0 + next_state = 1; + end + end + if state == 1 + if sync_bit == 1 + next_state = 2; + else + next_state = 0; + end + end + if state == 2 + if sync_bit == 0 + next_state = 3; + else + next_state = 0; + end + end + if state == 3 + if sync_bit == 1 + next_state = 4; + else + next_state = 0; + end + end + if state == 4 + if sync_bit == 0 + next_state = 5; + else + next_state = 0; + end + end + if state == 5 + if sync_bit == 1 + next_state = 6; + else + next_state = 0; + end + end + + % states 6 and above are track mode, make sure we keep getting 0101 sync bit sequence + + if state == 6 + if sync_bit == 0 + next_state = 7; + else + next_state = 0; + end + end + if state == 7 + if sync_bit == 1 + next_state = 6; + else + next_state = 0; + end + end + + state = next_state; + if state >= 6 + track = 1; + else + track = 0; + end +endfunction + diff --git a/codec2-dev/octave/fdmdv_demod.m b/codec2-dev/octave/fdmdv_demod.m index 0d52c5cf..a5baa49c 100644 --- a/codec2-dev/octave/fdmdv_demod.m +++ b/codec2-dev/octave/fdmdv_demod.m @@ -1,6 +1,6 @@ % fdmdv_demod.m % -% Demodulator function for FDMDV modem. Requires 48kHz sample rate raw files +% Demodulator function for FDMDV modem. Requires 8kHz sample rate raw files % as input % % Copyright David Rowe 2012 @@ -25,7 +25,7 @@ function fdmdv_demod(rawfilename, nbits) pilot_lut = generate_pilot_lut; pilot_lut_index = 1; - prev_pilot = zeros(1,M); + prev_pilot_lut_index = 3*M+1; % BER stats @@ -39,101 +39,93 @@ function fdmdv_demod(rawfilename, nbits) rx_symbols_log = []; rx_timing_log = []; foff_log = []; + rx_fdm_log = []; - % resampler states + % misc states - t = 3; - ratio = 1.000; - F=6; - MF=M*F; - nin = MF; - nin_size = MF+6; - buf_in = zeros(1,nin_size); - rx_fdm_buf = []; - ratio_log = []; + nin = M; % timing correction for sample rate differences + foff = 0; + track_log = []; + track = 1; + fest_state = 0; % Main loop ---------------------------------------------------- for f=1:frames - % update buf_in memory - - m = nin_size - nin; - for i=1:m - buf_in(i) = buf_in(i+nin); - end - % obtain nin samples of the test input signal - - for i=m+1:nin_size - buf_in(i) = fread(fin, 1, "short")/gain; + + for i=1:nin + rx_fdm(i) = fread(fin, 1, "short")/gain; end - % resample at 48kHz and decimate to 8kHz - - [rx_fdm_mf t nin] = resample(buf_in, t, 1.0, MF); - rx_fdm = rx_fdm_mf(1:F:MF); - - %rx_fdm = buf_in(m+1:m+n); - - %for i=1:M - % rx_fdm(i) = fread(fin, 1, "short")/gain; - %end - rx_fdm_buf = [rx_fdm_buf rx_fdm]; + rx_fdm_log = [rx_fdm_log rx_fdm(1:nin)]; % frequency offset estimation and correction - for i=1:M + for i=1:nin pilot(i) = pilot_lut(pilot_lut_index); pilot_lut_index++; if pilot_lut_index > 4*M pilot_lut_index = 1; end + prev_pilot(i) = pilot_lut(prev_pilot_lut_index); + prev_pilot_lut_index++; + if prev_pilot_lut_index > 4*M + prev_pilot_lut_index = 1; + end + end + foff_coarse = rx_est_freq_offset(rx_fdm, pilot, prev_pilot, nin); + if track == 0 + foff = foff_coarse; end - foff = rx_est_freq_offset(rx_fdm, pilot, prev_pilot); - prev_pilot = pilot; foff_log = [ foff_log foff ]; - foff = 0; foff_rect = exp(j*2*pi*foff/Fs); - for i=1:M + for i=1:nin foff_phase *= foff_rect'; rx_fdm(i) = rx_fdm(i)*foff_phase; end % baseband processing - rx_baseband = fdm_downconvert(rx_fdm); - rx_filt = rx_filter(rx_baseband); + rx_baseband = fdm_downconvert(rx_fdm, nin); + rx_filt = rx_filter(rx_baseband, nin); - [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband); + [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband, nin); rx_timing_log = [rx_timing_log rx_timing]; - beta = 1E-7; - ratio += beta*rx_timing; - if (ratio > 1.002) - ratio = 1.002; + nin = M; + if rx_timing > 2*M/P + nin += M/P; end - if (ratio < 0.998) - ratio = 0.998; + if rx_timing < 0; + nin -= M/P; end - ratio_log = [ratio_log ratio]; if strcmp(modulation,'dqpsk') rx_symbols_log = [rx_symbols_log rx_symbols.*conj(prev_rx_symbols)*exp(j*pi/4)]; else rx_symbols_log = [rx_symbols_log rx_symbols]; endif - [rx_bits sync] = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation); + [rx_bits sync f_err] = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation); + foff -= 0.5*f_err; prev_rx_symbols = rx_symbols; sync_log = [sync_log sync]; + % freq est state machine + + [track fest_state] = freq_state(sync, fest_state); + track_log = [track_log track]; + % count bit errors if we find a test frame % Allow 15 frames for filter memories to fill and time est to settle [test_frame_sync bit_errors] = put_test_bits(rx_bits); - if (test_frame_sync == 1 && f > 15) + if (test_frame_sync == 1) total_bit_errors = total_bit_errors + bit_errors; total_bits = total_bits + Ntest_bits; - bit_errors_log = [bit_errors_log bit_errors]; + bit_errors_log = [bit_errors_log bit_errors/Ntest_bits]; + else + bit_errors_log = [bit_errors_log 0]; end % test frame sync state machine, just for more informative plots @@ -173,50 +165,56 @@ function fdmdv_demod(rawfilename, nbits) % Plots % --------------------------------------------------------------------- + xt = (1:frames)/Rs; + secs = frames/Rs; + figure(1) clf; [n m] = size(rx_symbols_log); - plot(real(rx_symbols_log(1:Nc+1,10:m)),imag(rx_symbols_log(1:Nc+1,10:m)),'+') + plot(real(rx_symbols_log(1:Nc+1,15:m)),imag(rx_symbols_log(1:Nc+1,15:m)),'+') axis([-2 2 -2 2]); title('Scatter Diagram'); figure(2) clf; subplot(211) - plot(rx_timing_log(15:m)) + plot(xt, rx_timing_log) title('timing offset (samples)'); subplot(212) - plot(foff_log) + plot(xt, foff_log) title('Freq offset (Hz)'); grid figure(3) clf; subplot(211) - %plot(rx_fdm_buf); - %title('FDM Rx Signal'); - plot(ratio_log-1); - title('Sampling Clock error (ppm)'); + [a b] = size(rx_fdm_log); + xt1 = (1:b)/Fs; + plot(xt1, rx_fdm_log); + title('Rx FDM Signal'); subplot(212) Nfft=Fs; - S=fft(rx_fdm_buf,Nfft); + S=fft(rx_fdm_log,Nfft); SdB=20*log10(abs(S)); plot(SdB(1:Fs/4)) - title('FDM Tx Spectrum'); + title('FDM Rx Spectrum'); figure(4) clf; subplot(311) - stem(sync_log) - axis([0 frames 0 1.5]); + stem(xt, sync_log) + axis([0 secs 0 1.5]); title('BPSK Sync') subplot(312) - stem(bit_errors_log); + stem(xt, bit_errors_log); title('Bit Errors for test frames') subplot(313) - plot(test_frame_sync_log); - axis([0 frames 0 1.5]); + plot(xt, test_frame_sync_log); + axis([0 secs 0 1.5]); title('Test Frame Sync') - mean(ratio_log) + figure(5) + clf; + plot(xt, track_log); + axis([0 secs 0 1.5]); endfunction diff --git a/codec2-dev/octave/fdmdv_ut.m b/codec2-dev/octave/fdmdv_ut.m index 65cb0d54..ad6df497 100644 --- a/codec2-dev/octave/fdmdv_ut.m +++ b/codec2-dev/octave/fdmdv_ut.m @@ -11,9 +11,9 @@ fdmdv; % load modem code % Simulation Parameters -------------------------------------- -frames = 50; +frames = 100; EbNo_dB = 7.3; -Foff_hz = 0; +Foff_hz = 10; modulation = 'dqpsk'; hpa_clip = 150; @@ -30,6 +30,8 @@ rx_baseband_log = []; rx_bits_offset = zeros(Nc*Nb*2); prev_tx_symbols = ones(Nc+1,1); prev_rx_symbols = ones(Nc+1,1); +f_err = 0; +foff = 0; foff_log = []; tx_baseband_log = []; tx_fdm_log = []; @@ -47,6 +49,7 @@ test_frame_sync_state = 0; pilot_lut = generate_pilot_lut; pilot_lut_index = 1; +prev_pilot_lut_index = 3*M+1; % fixed delay simuation @@ -111,11 +114,12 @@ for f=1:frames % frequency offset + %Foff_hz += 1/Rs; + Foff = Foff_hz; for i=1:M % Time varying freq offset - % Foff = Foff_hz + 100*sin(t*2*pi/(300*Fs)); - % t++; - Foff = Foff_hz; + %Foff = Foff_hz + 100*sin(t*2*pi/(300*Fs)); + %t++; freq_offset = exp(j*2*pi*Foff/Fs); phase_offset *= freq_offset; tx_fdm(i) = phase_offset*tx_fdm(i); @@ -155,9 +159,13 @@ for f=1:frames if pilot_lut_index > 4*M pilot_lut_index = 1; end + prev_pilot(i) = pilot_lut(prev_pilot_lut_index); + prev_pilot_lut_index++; + if prev_pilot_lut_index > 4*M + prev_pilot_lut_index = 1; + end end - foff = rx_est_freq_offset(rx_fdm_delay, pilot, prev_pilot); - prev_pilot = pilot; + %foff = rx_est_freq_offset(rx_fdm_delay, pilot, prev_pilot, M); foff_log = [ foff_log foff ]; %foff = 0; foff_rect = exp(j*2*pi*foff/Fs); @@ -169,11 +177,11 @@ for f=1:frames % baseband processing - rx_baseband = fdm_downconvert(rx_fdm_delay(1:M)); + rx_baseband = fdm_downconvert(rx_fdm_delay(1:M), M); rx_baseband_log = [rx_baseband_log rx_baseband]; - rx_filt = rx_filter(rx_baseband); + rx_filt = rx_filter(rx_baseband, M); - [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband); + [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband, M); rx_timing_log = [rx_timing_log rx_timing]; %rx_phase = rx_est_phase(rx_symbols); @@ -185,10 +193,11 @@ for f=1:frames else rx_symbols_log = [rx_symbols_log rx_symbols]; endif - [rx_bits sync] = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation); + [rx_bits sync f_err] = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation); + foff -= 0.5*f_err; prev_rx_symbols = rx_symbols; sync_log = [sync_log sync]; - + % count bit errors if we find a test frame % Allow 15 frames for filter memories to fill and time est to settle @@ -246,16 +255,16 @@ figure(1) clf; [n m] = size(rx_symbols_log); plot(real(rx_symbols_log(1:Nc+1,15:m)),imag(rx_symbols_log(1:Nc+1,15:m)),'+') -axis([-2 2 -2 2]); +axis([-3 3 -3 3]); title('Scatter Diagram'); figure(2) clf; subplot(211) -plot(rx_timing_log) +plot(rx_timing_log(15:m)) title('timing offset (samples)'); subplot(212) -plot(foff_log) +plot(foff_log(15:m)) title('Freq offset (Hz)'); figure(3) -- 2.25.1