sim_in.Nsymbrow = Nsymbrow = Nsymb/Nc;
sim_in.Npilotsframe = Npilotsframe = 2;
sim_in.Nsymbrowpilot = Nsymbrowpilot = Nsymbrow + Npilotsframe;
-
+
if verbose == 2
printf("Each frame contains %d data bits or %d data symbols, transmitted as %d symbols by %d carriers.", framesize, Nsymb, Nsymbrow, Nc);
printf(" There are %d pilot symbols in each carrier together at the start of each frame, then %d data symbols.", Npilotsframe, Ns);
endfunction
-% Frequency offset estimation --------------------------------------------------
-
-% This function was used in initial Nov 2014 experiments
-
-function [f_max s_max] = freq_off_est(rx_fdm, tx_pilot, offset, n)
-
- Fs = 8000;
- nc = 1800; % portion we wish to correlate over (first 2 rows on pilots)
-
- % downconvert to complex baseband to remove images
-
- f = 1500;
- 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
-
-
-% Set of functions to implement latest and greatest freq offset
-% estimation, March 2015 ----------------------
-
-% 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;
- 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
- h = 0.5 - 0.5*cos(2*pi*(0:Ncm-1)/(Ncm-1));
- T = abs(fft(coarse_mem .* h, Ndft)).^2;
- sc = Ndft/Fs;
-
- 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);
-
- 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
-
-
function [ch_symb rx_timing rx_filt rx_baseband afdmdv f_est] = rate_Fs_rx_processing(afdmdv, ch_fdm_frame, f_est, nsymb, nin, freq_track)
M = afdmdv.M;
-
+
rx_baseband = [];
rx_filt = [];
rx_timing = [];
+ ch_fdm_frame_index = 1;
+
for r=1:nsymb
% shift signal to nominal baseband, this will put Nc/2 carriers either side of 0 Hz
- [rx_fdm_frame_bb afdmdv.fbb_phase_rx] = freq_shift(ch_fdm_frame(1+(r-1)*M:r*M), -f_est, afdmdv.Fs, afdmdv.fbb_phase_rx);
+ [rx_fdm_frame_bb afdmdv.fbb_phase_rx] = freq_shift(ch_fdm_frame(ch_fdm_frame_index:ch_fdm_frame_index + nin - 1), -f_est, afdmdv.Fs, afdmdv.fbb_phase_rx);
+ ch_fdm_frame_index += nin;
% downconvert each FDM carrier to Nc separate baseband signals
ch_symb(r,:) = rx_onesym;
+ % we only allow a timing shift on one symbol per frame
+
+ if nin != M
+ nin = M;
+ end
+
% freq tracking, see test_ftrack.m for unit test. Placed in
% this function as it needs to work on a symbol by symbol basis
% rather than frame by frame. This means the control loop
% tcohpsk.m
% David Rowe Oct 2014
%
-% Octave script that:
+% Octave coherent PSK modem script that hs two modes:
%
% 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.
+% the reference versions of the same modem written in Octave.
%
-% or (ii) can optionally be used to run an Octave version of the cohpsk
-% modem to tune and develop it.
+
+% (ii) Runs the Octave version of the cohpsk modem to tune and develop
+% it, including extensive channel simulations such as AWGN noise,
+% fading/HF, frequency offset, frequency drift, and tx/rx sample
+% rate differences.
% TODO:
%
-% [ ] Test
+% [X] Test
% [X] AWGN channel
% [X] freq offset
% [X] fading channel
% [X] freq drift
-% [ ] timing drift
+% [X] timing drift
% [X] tune perf/impl loss to get closer to ideal
% [X] linear interp of phase for better fading perf
% [X] freq offset/drift feedback loop
% [X] ability to "unsync" when signal dissapears
% [ ] some calibrated tests against FreeDV 1600
% + compare sound quality at various Es/Nos
+% [ ] sync
+% + set some req & implement
% [ ] way to handle eom w/o nasties
% + like mute ouput when signal has gone or v low snr
% + instantaneous snr
% [ ] nasty rig filter passband
+% [ ] pilot based EsNo estimation
graphics_toolkit ("gnuplot");
more off;
fading_en = 0;
hf_delay_ms = 2;
compare_with_c = 1;
+ sample_rate_ppm = 0;
end
% should be BER around 0.015 to 0.02
if strcmp(test, 'awgn')
frames = 100;
- foff = 0;
- dfoff = 0;
+ foff = 58.8;
+ dfoff = -0.5/Fs;
EsNodB = 8;
fading_en = 0;
hf_delay_ms = 2;
compare_with_c = 0;
+ sample_rate_ppm = 1000;
end
% Similar to AWGN - should be BER around 0.015 to 0.02
if strcmp(test, 'fading');
frames = 100;
foff = -53.1;
- dfoff = 0.5/Fs;
+ dfoff = 0.0/Fs;
EsNodB = 12;
fading_en = 1;
hf_delay_ms = 2;
compare_with_c = 0;
+ sample_rate_ppm = -100;
end
EsNo = 10^(EsNodB/10);
% main loop --------------------------------------------------------------------
+% run mod and channel as aseparate loop so we can resample to simulate sample rate differences
+
for f=1:frames
- acohpsk.frame = f;
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);
- %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);
ch_fdm_frame += noise;
ch_fdm_frame_log = [ch_fdm_frame_log ch_fdm_frame];
+end
+
+% simulate difference in sample clocks
+
+ch_fdm_frame_log = resample(ch_fdm_frame_log, (1E6 + sample_rate_ppm), 1E6);
+
+% Now run demod ----------------------------------------------------------------
+
+ch_fdm_frame_log_index = 1;
+f = 0;
+
+while (ch_fdm_frame_log_index + acohpsk.Nsymbrowpilot*M+M/P) < length(ch_fdm_frame_log)
+ acohpsk.frame = f++;
+
+ if sync == 0
+ nin = M;
+ else
+ nin = M;
+ if rx_timing(1) > M/P
+ nin = M + M/P;
+ end
+ if rx_timing(1) < -M/P
+ nin = M - M/P;
+ end
+ end
+
+ nin_frame = (acohpsk.Nsymbrowpilot-1)*M + nin;
+ ch_fdm_frame = ch_fdm_frame_log(ch_fdm_frame_log_index:ch_fdm_frame_log_index + nin_frame - 1);
+ ch_fdm_frame_log_index += nin_frame;
%
% Demod ----------------------------------------------------------------------
% store two frames of received samples so we can rewind if we get a good candidate
- ch_fdm_frame_buf(1:(Nsw-1)*acohpsk.Nsymbrowpilot*M) = ch_fdm_frame_buf(acohpsk.Nsymbrowpilot*M+1:Nsw*acohpsk.Nsymbrowpilot*M);
- ch_fdm_frame_buf((Nsw-1)*acohpsk.Nsymbrowpilot*M+1:Nsw*acohpsk.Nsymbrowpilot*M) = ch_fdm_frame;
+ ch_fdm_frame_buf(1:Nsw*acohpsk.Nsymbrowpilot*M-nin_frame) = ch_fdm_frame_buf(nin_frame+1:Nsw*acohpsk.Nsymbrowpilot*M);
+ ch_fdm_frame_buf(Nsw*acohpsk.Nsymbrowpilot*M-nin_frame+1:Nsw*acohpsk.Nsymbrowpilot*M) = ch_fdm_frame;
next_sync = sync;
- nin = M;
% if out of sync do Initial Freq offset estimation over NSW frames to flush out memories
% BER stats
if f > 2
- error_positions = xor(prev_tx_bits2, rx_bits);
- %error_positions = xor(prev_tx_bits, rx_bits);
+ error_positions = xor(tx_bits_log((f-3)*framesize+1:(f-2)*framesize), rx_bits);
Nerrs += sum(error_positions);
nerr_log = [nerr_log sum(error_positions)];
Tbits += length(error_positions);
printf("\r [%d]", f);
end
-
- %rx_fdm_frame_bb_log = [rx_fdm_frame_bb_log rx_fdm_frame_bb];
- %rx_baseband_log = [rx_baseband_log rx_baseband];
- %rx_filt_log = [rx_filt_log rx_filt];
- %rx_timing_log = [rx_timing_log rx_timing];
- %ch_symb_log = [ch_symb_log; ch_symb];
- % ct_symb_ff_log = [ct_symb_ff_log; acohpsk.ct_symb_ff_buf(1:acohpsk.Nsymbrowpilot,:)];
+ % reset BER stats if we lose sync
-
- if sync == 0
- Nerrs = 0;
- Tbits = 0;
- nerr_log = [];
+ if sync == 1
+ %Nerrs = 0;
+ %Tbits = 0;
+ %nerr_log = [];
end
- % printf("f: %d sync: %d next_sync: %d\n", f, sync, next_sync);
[sync acohpsk] = sync_state_machine(acohpsk, sync, next_sync);
prev_tx_bits2 = prev_tx_bits;
plot(foff_log(1:length(f_est_log)) - f_est_log + Fcentre)
title('freq offset estimation error');
+ figure(8)
+ clf
+ subplot(211)
+ stem(real(rx_filt_log(1,52:4:200)))
+ subplot(212)
+ stem(imag(rx_filt_log(1,52:4:200)))
end
projects / freetel-svn-tracking.git / commitdiff