How can I perform timing recovery from the envelope of PSK signals?

Question

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):

1. sum the envelopes of each subcarrier
2. this sum has a frequency component at the symbol rate, so determine the phase of that
3. 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