I'm trying to reimplement FDMDV in GNU Radio, largely for educational purposes. It consists of 15 PSK subcarriers, 75 Hz apart, 50 symbols/sec each. The center one is BPSK and alternates phase each symbol, and the rest are QPSK and carry the data.
I'm somewhat stuck on how to perform timing recovery if this signal. Since each subcarrier has the same clock, it seems this should be worth something. The existing implementation seems to work by (full implementation at the end):
- sum the envelopes of each subcarrier
- this sum has a frequency component at the symbol rate, so determine the phase of that
- drive this phase to a fixed value through a control loop
How might I implement this in GNU Radio? Is this something that can be accomplished with the standard blocks or will I need to implement a new block of my own? Or is there another method I could use to perform the same task?
function [rx_symbols rx_timing_M env fdmdv] = rx_est_timing(fdmdv, rx_filt, nin)
samples_per_symbol = fdmdv.M;
Nt = fdmdv.Nt;
num_carriers = fdmdv.Nc;
rx_filter_mem_timing = fdmdv.rx_filter_mem_timing;
P = fdmdv.P;
Nfilter = fdmdv.Nfilter;
Nfiltertiming = fdmdv.Nfiltertiming;
% nin adjust
% --------------------------------
% 120 -1 (one less rate P sample)
% 160 0 (nominal)
% 200 1 (one more rate P sample)
adjust = P - nin*P/samples_per_symbol;
% update buffer of Nt rate P filtered symbols
rx_filter_mem_timing(:,1:(Nt-1)*P+adjust) = rx_filter_mem_timing(:,P+1-adjust:Nt*P);
rx_filter_mem_timing(:,(Nt-1)*P+1+adjust:Nt*P) = rx_filt(:,:);
% sum envelopes of all carriers
env = sum(abs(rx_filter_mem_timing(:,:))); % use all num_carriers+1 carriers for timing
%env = abs(rx_filter_mem_timing(num_carriers+1,:)); % just use BPSK pilot
[n m] = size(env);
% The envelope has a frequency component at the symbol rate. The
% phase of this frequency component indicates the timing. So work out
% single DFT at frequency 2*pi/P
x = env * exp(-j*2*pi*(0:m-1)/P)';
norm_rx_timing = angle(x)/(2*pi);
rx_timing = norm_rx_timing*P + P/4;
if (rx_timing > P)
rx_timing -= P;
end
if (rx_timing < -P)
rx_timing += P;
end
% rx_filter_mem_timing contains Nt*P samples (Nt symbols at rate P),
% where Nt is odd. Lets use linear interpolation to resample in the
% centre of the timing estimation window
rx_timing += floor(Nt/2)*P;
low_sample = floor(rx_timing);
fract = rx_timing - low_sample;
high_sample = ceil(rx_timing);
%printf("rx_timing: %f low_sample: %f high_sample: %f fract: %f\n", rx_timing, low_sample, high_sample, fract);
rx_symbols = rx_filter_mem_timing(:,low_sample)*(1-fract) + rx_filter_mem_timing(:,high_sample)*fract;
% rx_symbols = rx_filter_mem_timing(:,high_sample+1);
rx_timing_M = norm_rx_timing*samples_per_symbol;
fdmdv.rx_filter_mem_timing = rx_filter_mem_timing;
endfunction