global Nfilter = Nsym*M;
global Nfiltertiming = M+Nfilter+M;
alpha = 0.5;
-global snr_coeff = 0.9 % SNR est averaging filter coeff
+global snr_coeff;
+snr_coeff = 0.9; % SNR est averaging filter coeff
+global Nph;
+Nph = 9; % number of symbols to estimate phase over
+ % must be odd number as we take centre symbol
% root raised cosine (Root Nyquist) filter
endfunction
-% Phase estimation function - probably won't work over a HF channel
-% Tries to operate over a single symbol but uses phase information from
-% all Nc carriers which should increase the SNR of phase estimate.
-% Maybe phase is coherent over a couple of symbols in HF channel,not
-% sure but it's worth 3dB so worth experimenting or using coherent as
-% an option.
+% Experimental "feed forward" phase estimation function - estimates
+% phase over a windows of Nph (e.g. Nph = 9) symbols. May not work
+% well on HF channels but lets see. Has a phase ambiguity of m(pi/4)
+% where m=0,1,2 which needs to be corrected outside of this function
-function rx_phase = rx_est_phase(prev_rx_symbols, rx_symbols)
+function [phase_offsets ferr] = rx_est_phase(rx_symbols)
+ global rx_symbols_mem;
+ global prev_phase_offsets;
+ global phase_amb;
+ global Nph;
+ global Nc;
- % modulation strip
+ % keep record of Nph symbols
- rx_phase = angle(sum(rx_symbols .^ 4))/4;
+ rx_symbols_mem(:,1:Nph-1) = rx_symbols_mem(:,2:Nph);
+ rx_symbols_mem(:,Nph) = rx_symbols;
+ % estimate and correct phase offset based of modulation stripped samples
+
+ phase_offsets = zeros(Nc+1,1);
+ for c=1:Nc+1
+
+ % rotate QPSK constellation to a single point
+ mod_stripped = abs(rx_symbols_mem(c,:)) .* exp(j*4*angle(rx_symbols_mem(c,:)));
+
+ % find average phase offset, which will be on -pi/4 .. pi/4
+ sum_real = sum(real(mod_stripped));
+ sum_imag = sum(imag(mod_stripped));
+ phase_offsets(c) = atan2(sum_imag, sum_real)/4;
+
+ % determine if phase has jumped from - -> +
+ if (prev_phase_offsets(c) < -pi/8) && (phase_offsets(c) > pi/8)
+ phase_amb(c) -= pi/2;
+ if (phase_amb(c) < -pi)
+ phase_amb(c) += 2*pi;
+ end
+ end
+
+ % determine if phase has jumped from + -> -
+ if (prev_phase_offsets(c) > pi/8) && (phase_offsets(c) < -pi/8)
+ phase_amb(c) += pi/2;
+ if (phase_amb(c) > pi)
+ phase_amb(c) -= 2*pi;
+ end
+ end
+ end
+
+ ferr = mean(phase_offsets - prev_phase_offsets);
+ prev_phase_offsets = phase_offsets;
+
endfunction
% map (Nc,1) DQPSK symbols back into an (1,Nc*Nb) array of bits
for c=1:Nc
+ msb = lsb = 0;
d = phase_difference(c);
if ((real(d) >= 0) && (imag(d) >= 0))
msb = 0; lsb = 0;
% vector of mags, one for each carrier.
s = abs(phase_difference);
-
+
% signal mag estimate for each carrier is a smoothed version
% of instantaneous magntitude, this gives us a vector of smoothed
% mag estimates, one for each carrier.
-
+
sig_est = snr_coeff*sig_est + (1 - snr_coeff)*s;
+ %printf("s: %f sig_est: %f snr_coeff: %f\n", s(1), sig_est(1), snr_coeff);
+
% noise mag estimate is distance of current symbol from average
% location of that symbol. We reflect all symbols into the first
% quadrant for convenience.
% noise power estimates, one for each carrier.
noise_est = snr_coeff*noise_est + (1 - snr_coeff)*n;
+
endfunction
for i=1:nbits
bits(i) = test_bits(current_test_bit++);
+ %if (mod(i,2) == 0)
+ % bits(i) = 1;
+ %else
+ % bits(i) = 0;
+ %end
+
if (current_test_bit > Ntest_bits)
current_test_bit = 1;
endif
% Accepts nbits from rx and attempts to sync with test_bits sequence.
% if sync OK measures bit errors
-function [sync bit_errors] = put_test_bits(rx_bits)
+function [sync bit_errors] = put_test_bits(test_bits, rx_bits)
global Ntest_bits; % length of test sequence
- global test_bits;
global rx_test_bits_mem;
% Append to our memory
global rx_test_bits_mem;
rx_test_bits_mem = zeros(1,Ntest_bits);
+% Experimental phase estimator states ----------------------
+
+global rx_symbols_mem;
+rx_symbols_mem = zeros(Nc+1, Nph);
+global prev_phase_offsets;
+prev_phase_offsets = zeros(Nc+1, 1);
+global phase_amb;
+phase_amb = zeros(Nc+1, 1);
% count bit errors if we find a test frame
- [test_frame_sync bit_errors] = put_test_bits(rx_bits);
+ [test_frame_sync bit_errors] = put_test_bits(test_bits, rx_bits);
if (test_frame_sync == 1)
total_bit_errors = total_bit_errors + bit_errors;
total_bits = total_bits + Ntest_bits;
--- /dev/null
+% fdmdv_demod_coh.m
+%
+% Demodulator function for FDMDV modem (Octave version). Requires
+% 8kHz sample rate raw files as input. This version uses experimental
+% psuedo coherent demodulation.
+%
+% Copyright David Rowe 2013
+% This program is distributed under the terms of the GNU General Public License
+% Version 2
+%
+
+function fdmdv_demod_coh(rawfilename, nbits, pngname)
+
+ fdmdv; % include modem code
+
+ modulation = 'dqpsk';
+
+ fin = fopen(rawfilename, "rb");
+ gain = 1000;
+ frames = nbits/(Nc*Nb);
+
+ prev_rx_symbols = ones(Nc+1,1);
+ foff_phase = 1;
+
+ % BER stats
+
+ total_bit_errors = 0;
+ total_bits = 0;
+ bit_errors_log = [];
+ sync_log = [];
+ test_frame_sync_log = [];
+ test_frame_sync_state = 0;
+
+ % SNR states
+
+ sig_est = zeros(Nc+1,1);
+ noise_est = zeros(Nc+1,1);
+
+ % logs of various states for plotting
+
+ rx_symbols_log = [];
+ rx_timing_log = [];
+ foff_log = [];
+ rx_fdm_log = [];
+ snr_est_log = [];
+
+ % misc states
+
+ nin = M; % timing correction for sample rate differences
+ foff = 0;
+ track_log = [];
+ track = 0;
+ fest_state = 0;
+
+ % psuedo coherent demod states
+
+ rx_symbols_ph_log = [];
+ prev_rx_symbols_ph = ones(Nc+1,1);
+ rx_phase_offsets_log = [];
+ phase_amb_log = [];
+
+ % Main loop ----------------------------------------------------
+
+ for f=1:frames
+
+ % obtain nin samples of the test input signal
+
+ for i=1:nin
+ rx_fdm(i) = fread(fin, 1, "short")/gain;
+ end
+
+ rx_fdm_log = [rx_fdm_log rx_fdm(1:nin)];
+
+ % frequency offset estimation and correction
+
+ [pilot prev_pilot pilot_lut_index prev_pilot_lut_index] = get_pilot(pilot_lut_index, prev_pilot_lut_index, nin);
+ [foff_coarse S1 S2] = rx_est_freq_offset(rx_fdm, pilot, prev_pilot, nin);
+
+ if track == 0
+ foff = foff_coarse;
+ end
+ foff_log = [ foff_log foff ];
+ foff_rect = exp(j*2*pi*foff/Fs);
+
+ 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, nin);
+ rx_filt = rx_filter(rx_baseband, nin);
+
+ [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband, nin);
+ rx_timing_log = [rx_timing_log rx_timing];
+
+ nin = M;
+ if rx_timing > 2*M/P
+ nin += M/P;
+ end
+ if rx_timing < 0;
+ nin -= M/P;
+ end
+
+ rx_symbols_log = [rx_symbols_log rx_symbols.*(conj(prev_rx_symbols)./abs(prev_rx_symbols))*exp(j*pi/4)];
+
+ % coherent phase offset estimation ------------------------------------
+
+ [rx_phase_offsets ferr] = rx_est_phase(rx_symbols);
+ rx_phase_offsets_log = [rx_phase_offsets_log rx_phase_offsets];
+ phase_amb_log = [phase_amb_log phase_amb];
+ rx_symbols_ph = rx_symbols_mem(:,floor(Nph/2)+1) .* exp(-j*(rx_phase_offsets + phase_amb));
+ rx_symbols_ph_log = [rx_symbols_ph_log rx_symbols_ph .* exp(j*pi/4)];
+ rx_symbols_ph = -1 + 2*(real(rx_symbols_ph .* exp(j*pi/4)) > 0) + j*(-1 + 2*(imag(rx_symbols_ph .* exp(j*pi/4)) > 0));
+
+ % Std differential (used for freq offset est and BPSK sync) and psuedo coherent detection -----------------------
+
+ [rx_bits_unused sync f_err pd ] = qpsk_to_bits(prev_rx_symbols, rx_symbols, modulation);
+ [rx_bits sync_unused ferr_unused pd_unused] = qpsk_to_bits(prev_rx_symbols_ph, rx_symbols_ph, 'dqpsk');
+
+ foff -= 0.5*f_err;
+ prev_rx_symbols = rx_symbols;
+ prev_rx_symbols_ph = rx_symbols_ph;
+ sync_log = [sync_log sync];
+
+ [sig_est noise_est] = snr_update(sig_est, noise_est, pd);
+ snr_est = calc_snr(sig_est, noise_est);
+ snr_est_log = [snr_est_log snr_est];
+
+ % 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
+
+ [test_frame_sync bit_errors] = put_test_bits(test_bits, rx_bits);
+ 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/Ntest_bits];
+ else
+ bit_errors_log = [bit_errors_log 0];
+ end
+
+ % test frame sync state machine, just for more informative plots
+
+ next_test_frame_sync_state = test_frame_sync_state;
+ if (test_frame_sync_state == 0)
+ if (test_frame_sync == 1)
+ next_test_frame_sync_state = 1;
+ test_frame_count = 0;
+ end
+ end
+
+ if (test_frame_sync_state == 1)
+ % we only expect another test_frame_sync pulse every 4 symbols
+ test_frame_count++;
+ if (test_frame_count == 4)
+ test_frame_count = 0;
+ if ((test_frame_sync == 0))
+ next_test_frame_sync_state = 0;
+ 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
+
+ % ---------------------------------------------------------------------
+ % Print Stats
+ % ---------------------------------------------------------------------
+
+ ber = total_bit_errors / total_bits;
+
+ printf("%d bits %d errors BER: %1.4f\n",total_bits, total_bit_errors, ber);
+
+ % ---------------------------------------------------------------------
+ % 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,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(xt, rx_timing_log)
+ title('timing offset (samples)');
+ subplot(212)
+ plot(xt, foff_log, '-;freq offset;')
+ hold on;
+ plot(xt, track_log*75, 'r;course-fine;');
+ hold off;
+ title('Freq offset (Hz)');
+ grid
+
+ figure(3)
+ clf;
+ subplot(211)
+ [a b] = size(rx_fdm_log);
+ xt1 = (1:b)/Fs;
+ plot(xt1, rx_fdm_log);
+ title('Rx FDM Signal');
+ subplot(212)
+ spec(rx_fdm_log,8000);
+ title('FDM Rx Spectrogram');
+
+ figure(4)
+ clf;
+ subplot(311)
+ stem(xt, sync_log)
+ axis([0 secs 0 1.5]);
+ title('BPSK Sync')
+ subplot(312)
+ stem(xt, bit_errors_log);
+ title('Bit Errors for test frames')
+ subplot(313)
+ plot(xt, test_frame_sync_log);
+ axis([0 secs 0 1.5]);
+ title('Test Frame Sync')
+
+ figure(5)
+ clf;
+ plot(xt, snr_est_log);
+ title('SNR Estimates')
+
+ figure(6)
+ clf;
+ [n m] = size(rx_symbols_ph_log);
+ plot(real(rx_symbols_ph_log(1:Nc+1,15:m)),imag(rx_symbols_ph_log(1:Nc+1,15:m)),'+')
+ %plot(real(rx_symbols_ph_log(2,15:m)),imag(rx_symbols_ph_log(2,15:m)),'+')
+ axis([-2 2 -2 2]);
+ title('Scatter Diagram - after phase correction');
+
+ figure(7)
+ clf;
+ subplot(211)
+ plot(rx_phase_offsets_log(1,:))
+ subplot(212)
+ plot(phase_amb_log(1,:))
+ title('Rx Phase Offset Est')
+
+endfunction
% Simulation Parameters --------------------------------------
-frames = 25;
+frames = 100;
EbNo_dB = 7.3;
Foff_hz = 0;
modulation = 'dqpsk';
track = 0;
track_log = [];
+snr_log = [];
% ---------------------------------------------------------------------
% Main loop
% Update SNR est
[sig_est noise_est] = snr_update(sig_est, noise_est, pd);
+ snr_log = [snr_log calc_snr(sig_est, noise_est)];
% count bit errors if we find a test frame
% Allow 15 frames for filter memories to fill and time est to settle
axis([0 frames 0 1.5]);
title('Test Frame Sync')
+figure(5)
+clf
+plot(snr_log)
% Simulation Parameters --------------------------------------
-frames = 100;
-EbNo_dB = 73;
-Foff_hz = 0;
+frames = 200;
+EbNo_dB = 7;
+Foff_hz = -100;
hpa_clip = 150;
% ------------------------------------------------------------
snr_log = [];
-% PSK demod ----
-
-% Set up matrix of tx symbols, one row for each carrier.
-
-% The differential PSK encoder is driven by a sequence of test bits.
-% This test sequence is 4 frames (122 bits long), which means it
-% repeats every 4 frames. However the DPSK encoder also has 4 states,
-% so it turns out that the transmitted sequence of PSK symbols repeats
-% itself every 16 frames. Depending on the test data some carriers
-% repeat more often than that, but the sequence as a whole repeats
-% every 16 frames.
-
-% Each row in the matrix is the sequence of symbols for one carrier.
-% This repeats itself every 16 symbols. We would like to demodulate
-% and measure the BER of this block of 16 symbols. To do this we need
-% to correct any phase offset.
-
-prev_psk_test_symbols = ones(Nc+1,1);
-tx_symbols_mem = [];
-for i=1:6
- psk_test_symbols = generate_psk_test_symbols(prev_psk_test_symbols);
- prev_psk_test_symbols = psk_test_symbols(:,4);
- tx_symbols_mem = [tx_symbols_mem psk_test_symbols .* exp(j*pi/4)];
-end
-
-% PSK BER counters
-
-psk_total_bits = 0;
-psk_total_bit_errors = 0;
-
rx_symbols_ph_log = [];
+prev_rx_symbols_ph = ones(Nc+1,1);
+rx_phase_offsets_log = [];
+phase_amb_log = [];
% ---------------------------------------------------------------------
% Main loop
foff = foff_course;
end
- foff = 0; % disable for now
+ %foff = 0; % disable for now
foff_log = [ foff_log foff ];
foff_rect = exp(j*2*pi*foff/Fs);
rx_filt = rx_filter(rx_baseband, M);
[rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband, M);
+ rx_symbols_log = [rx_symbols_log rx_symbols.*(conj(prev_rx_symbols)./abs(prev_rx_symbols))*exp(j*pi/4)];
rx_timing_log = [rx_timing_log rx_timing];
- [rx_bits sync foff_fine pd] = qpsk_to_bits(prev_rx_symbols, rx_symbols, 'dqpsk');
- rx_symbols_log = [rx_symbols_log rx_symbols .* exp(j*pi/4)];
+ % coherent phase offset estimation ------------------------------------
+
+ [rx_phase_offsets ferr] = rx_est_phase(rx_symbols);
+ rx_phase_offsets_log = [rx_phase_offsets_log rx_phase_offsets];
+ phase_amb_log = [phase_amb_log phase_amb];
+ rx_symbols_ph = rx_symbols_mem(:,floor(Nph/2)+1) .* exp(-j*(rx_phase_offsets + phase_amb));
+ rx_symbols_ph_log = [rx_symbols_ph_log rx_symbols_ph .* exp(j*pi/4)];
+ rx_symbols_ph = -1 + 2*(real(rx_symbols_ph .* exp(j*pi/4)) > 0) + j*(-1 + 2*(imag(rx_symbols_ph .* exp(j*pi/4)) > 0));
+
+ % Std differential (used for freq offset est and BPSK sync) and psuedo coherent detection -----------------------
+
+ [rx_bits_unused sync ferr pd] = qpsk_to_bits(prev_rx_symbols, rx_symbols, 'dqpsk');
+ [rx_bits sync_unused ferr_unused pd] = qpsk_to_bits(prev_rx_symbols_ph, rx_symbols_ph, 'dqpsk');
+
+ %----------------------------------------------------------------------
foff -= 0.5*ferr;
prev_rx_symbols = rx_symbols;
+ prev_rx_symbols_ph = rx_symbols_ph;
sync_log = [sync_log sync];
% freq est state machine
[sig_est noise_est] = snr_update(sig_est, noise_est, pd);
snr_log = [snr_log calc_snr(sig_est, noise_est)];
- % coherent PSK demod --------------------------------------
-
- % update memory of the last Nph rx symbols
-
- rx_symbols_mem(:,1:Npom-1) = rx_symbols_mem(:,2:Npom);
- rx_symbols_mem(:,Npom) = rx_symbols .* exp(j*pi/4);
-
- % estimate and correct phase of each received symbol in 16 symbol frame
-
- for c=1:Nc
- for s=Nph2+1:Nph2+Nstates % symbols in middle of matrix that we correct
- st = s - Nph2; en = s + Nph2;
- corr = rx_symbols_mem(c,st:en) * tx_symbols_mem(c,st:en)';
- rx_symbols_mem_ph(c,s) = rx_symbols_mem(c,s) .* exp(-j*angle(corr));
- end
- end
-
- % count bit errors in centre of window
-
- psk_bit_errors = 0;
- psk_bits = Nc*Nstates;
- st = Nph2+1; en = Nph2+Nstates;
- for c=1:Nc
- psk_bit_errors += sum(real(rx_symbols_mem_ph(c,st:en)) .* real(tx_symbols_mem(c,st:en)) < 0);
- psk_bit_errors += sum(imag(rx_symbols_mem_ph(c,st:en)) .* imag(tx_symbols_mem(c,st:en)) < 0);
- end
- printf("psk_bit_errors: %d corr: %f\n", psk_bit_errors, angle(corr));
- ber = psk_bit_errors/psk_bits;
- if ber < 0.1
- psk_total_bit_errors += psk_bit_errors;
- psk_total_bits += psk_bits;
- % phase correction only valid when we are in sync
- rx_symbols_ph_log = [rx_symbols_ph_log rx_symbols_mem_ph(:,st:en)];
- end
-
- % ---------------------------------------------------------
-
% count bit errors if we find a test frame
[test_frame_sync bit_errors] = put_test_bits(test_bits, rx_bits);
- if test_frame_sync == 1
+ if (test_frame_sync == 1) && (f > 15)
total_bit_errors = total_bit_errors + bit_errors;
total_bits = total_bits + Ntest_bits;
bit_errors_log = [bit_errors_log bit_errors];
printf(" bits......: %d\n", total_bits);
printf(" errors....: %d\n", total_bit_errors);
printf(" BER.......: %1.4f\n", ber);
-printf("PSK\n");
-printf(" bits......: %d\n", psk_total_bits);
-printf(" errors....: %d\n", psk_total_bit_errors);
-printf(" BER.......: %1.4f\n", psk_total_bit_errors/psk_total_bits);
% ---------------------------------------------------------------------
% Plots
figure(3)
clf;
-subplot(211)
-plot(real(tx_fdm_log));
-title('FDM Tx Signal');
-subplot(212)
-Nfft=Fs;
-S=fft(rx_fdm_log,Nfft);
-SdB=20*log10(abs(S));
-plot(SdB(1:Fs/4))
-title('FDM Rx Spectrum');
-
-figure(4)
-clf;
subplot(311)
stem(sync_log)
axis([0 frames 0 1.5]);
axis([0 frames 0 1.5]);
title('Test Frame Sync')
-figure(5)
-clf;
-subplot(211)
-stem(real(rx_symbols_mem(2,:)))
-subplot(212)
-stem(imag(rx_symbols_mem(2,:)))
-
-figure(6)
-clf;
-subplot(211)
-stem(real(tx_symbols_mem(2,:)))
-subplot(212)
-stem(imag(tx_symbols_mem(2,:)))
-
-figure(7)
+figure(4)
clf;
[n m] = size(rx_symbols_ph_log);
plot(real(rx_symbols_ph_log(1:Nc+1,15:m)),imag(rx_symbols_ph_log(1:Nc+1,15:m)),'+')
%plot(real(rx_symbols_ph_log(2,15:m)),imag(rx_symbols_ph_log(2,15:m)),'+')
axis([-3 3 -3 3]);
title('Scatter Diagram - after phase correction');
+
+figure(5)
+clf;
+subplot(211)
+plot(rx_phase_offsets_log(1,:))
+subplot(212)
+plot(phase_amb_log(1,:))
+title('Rx Phase Offset Est')
projects / freetel-svn-tracking.git / commitdiff