# Can I search Doppler shift in freqency domain

I would like to know if it is possible to search Doppler shift of a carrier with sufficient oversampling while performing code search in the freqency domain. For instance can I selectively drop every X time domain sample bin from the fft processing? Will that not effectively change the input frequency being compared to the locally generated code? I would like to perform the Doppler shift search of digital satellite broadcasts without having to change the lo supplying the IF in hardware. Can this be done this way?

• A tried and true brute force method is to implement multiple parallel code search branches spaced in frequency by the baud rate. You need enough branches to cover your frequency uncertainty.
– John
Jan 6, 2014 at 18:48
• OK so if I generate several code arrays each with varied timedelays between chips. The point of greatest correlation when i multiply my generated code ffts each to the received fft should represent the actual frequency that the received carrier is Doppler shifted to? Jan 8, 2014 at 13:44
• What I was referring to is a set of PN codes separated in frequency by mixing with a complex exponential. Correlate with each of them and the one that peaks the most gives an indication of frequency. You can do this with one forward FFT and multiple inverse FFTs.
– John
Jan 8, 2014 at 14:59
• Can I runs this by you then? If I understand correctly I would take my FFT generated freq domain bins of local code and multiply by e(-1*sqrt(-1)*2*pi* shift_factor / N *k) to get a frequency shifted code? Do I need to limit the value of k to just a few values? I have read some posts where k needed to be symetrical values about 0, for instance a set of [0, 1, 2, -2, -1]? What effect does increasing the number of k values have? Jan 9, 2014 at 12:34
• This is the way I do it. Start with the basic PN code with one sample per chip. Time reverse it, and oversample it to 4 samples per chip. Define a set of search frequencies, say f=[-456,-123,0,123,456] Hz. Now, multiply the time-reversed, oversampled PN code by exp(j*2*pi*f*[0:N-1]/fs) for each of the f values, where fs is 4x the chip rate. Take the FFT of each freq-shifted PN codes and save the FFTs off in a table somewhere. Use the FFT table to multiply with the forward FFT of your incoming data. Notice that the specific search frequencies are arbitrary.
– John
Jan 9, 2014 at 13:44

## 3 Answers

Note my function below showing how I simultaneously searched for code and carrier (Doppler) offset for spread spectrum signals such as GPS in one 2D IFFT computation*. The 2D IFFT is processing intensive, but this approach has a lot of merit if you need very fast ("one-shot") acquisition. (*Note the total processing involves two 1D FFT's, an element by element matrix multiplication for all Doppler and Code bins, and a 2D IFFT).

The basic premise is using fft's to cross-correlate two sequences (your waveform that has code and carrier offsets, and your reference code):

$\textrm{XCORR} = \textrm{ifft}(\textrm{fft}(a)\textrm{fft}(b)^*)$

Specifically the circular cross correlation between two sequences is equal to the inverse fft of the complex conjugate multiplication of their respective fft's.

If you circularly shift one of the fft's before multiplying, you calculate the result for a different Doppler bin! It is that easy. For actual implementation detail I included my Matlab function below. Here is an example plot of the result:

One FFT is computed for the received sequence, and one FFT is computed for the reference code. A matrix of delayed FFT's is created for one of the two that only involves replicating columns with a circular shift of the FFT already computed (which is a Hankel matrix) , and then a matrix of the other is created by simply replicating the FFT computed in each column, matching the number of Doppler columns created. The two matrices are complex conjugate multiplied element by element ($M_{ij}= A_{ij}B_{ij}^*$) Then (this is where all the processing is needed), a 2D inverse FFT is computed of the product matrix resulting in the solution depicted in the surface plot above.

On one axis for the figure above is the normalized frequency offset (Doppler) and the other axis is the relative code delay, and the vertical axis is the correlation where the delay and Doppler were simultaneously determined. So with one capture (whose length depends only on post correlation SNR required of the received sequence) you can complete acquisition processing with the use of two ffts, an element by element matrix multiplication and a 2D IFFT.

function [fdout,freq_axis]= fdcorr(a,b,frange)
%FDCORR Joint frequency and delay correlation
%   Performs an efficient circular cross correlation between
%   two vectors over all possible frequency offsets (Dopplers)
%   between -pi and pi where 2pi is the sampling rate.
%
%[FDOUT, FREQ_AXIS] = FDCORR(A,B)
%
%   A:   vector representing received signal sampled at FS = 2pi
%   B:   vector representing correlation code sampled at FS = 2pi
%
%   FDOUT: 2D Correlation surface vs delay and frequency offset
%   FREQ_AXIS: Frequency values aligned to frequency axis (-pi to pi)
%
%[FDOUT, FREQ_AXIS] = FDCORR(A,B,FRANGE)
%
%   Limits frequency range as defined in FRANGE.    Typically
%   frequency offsets are much smaller than the Nyquist bandwidth, so
%   using a reduced range for FRANGE from -pi to pi will be more
%   efficient.
%
%  FRANGE: Two element vector from lowest to highest frequency to
%          search over within a frequency range from -pi to pi.
%          Example [-pi/10 pi/10]
%
%FDCORR(A,B,<FRANGE>) with no output arguments plots the surface
%       correlation magnitude vs frequency and delay offset in the
%       current figure window.  FRANGE optional.
%
%
% Dan Boschen 11/22/2009

%==========================================================
% VERSION HISTORY
%
% 11/22/2009     Inital Release
%
%
%==========================================================

%==========================================================
% APPLICATION NOTES:
%
% Vector "A" represents a received signal that has been encoded with
% a correlation code (such as GPS and CDMA applications) or a training
% sequence.  It has an unknown frequency offset due to Doppler offsets
% and frequency offsets in the receiver, and an unkown delay shift
% relative to the start of the code sequence represented by vector "B".
% FDCORR will efficiently find the correlation between the received
% sequence and reference code over frequency and delay offsets.
% Vectors A and B are interchangable.
%
% To detect the maximum correlation value and delay and freq location use:
% [temp1, temp2]=max(fdout);
% [maxcorr,freq_index]=max(temp1);
% delay =temp2(freq_index)-1;
% freq= freq_axis(freq_index);
%
% The basic operation uses the Cross-Correlation theorem:
% xcorr= ifft(fft(a).*conj(fft(b)))
%
% Note the detected frequency offset is course and is the
% frequency of the closest Doppler bin where each Doppler bin
% is seperated by 2pi/N
%
%==========================================================

%
%error handling
if nargin==3
if length(frange)~= 2
error('FRANGE must be a two element vector')
elseif abs(frange)> pi
error('FRANGE values must be between -pi and pi')
end
end

%force column vectors:
a=a(:);
b=b(:);

%============================================================
%create all Doppler versions of a by shifting the fft:
%============================================================

N=max(length(a),length(b));
fft_a= fft(a,N);

%frequency axis for all doppler bins -pi to pi

%index is the "fft shifted" freq axis. Method below faster than
%using fftshift(0:N-1), but same result.
freq_index=[-floor((N-1)/2):floor(N/2)];
freq_axis= freq_index*2*pi/N;   %convert to -pi to pi range

if nargin==3     %reduced doppler range
freq_include=find(freq_axis>=frange(1) & freq_axis<=frange(2));
freq_axis=freq_axis(freq_include);
freq_index=freq_index(freq_include);
num_freqs=length(freq_index);
else
num_freqs=N;
end

%create Hankel matrix for all Dopplers (each column is the fft
%of vector a at a different Doppler offset)
fft_idx=(1:N)';
fft_dopplers = fft_idx(:,ones(num_freqs,1))+freq_index(ones(N,1),:);
fft_dopplers= mod(fft_dopplers-1,N)+1;    %Hankel subscripts
fft_dopplers(:)=  fft_a(fft_dopplers);      %fill data

%alternate method if solving for all Doppler bins (processing intensive):
%fft_all_dopplers= hankel(fft_a,[fft_a(end); fft_a(2:end)]);

%==================================================================
%Circular Correlation using circ_corr= ifft(fft(a).*conj(fft(b))
%==================================================================
fft_code_mtx= fft(b,N)*ones(1,num_freqs);
fdout=(ifft(fft_code_mtx.*conj(fft_dopplers)));

%===============================================================
%plot result
%===============================================================
if nargout ==0
surf(freq_axis,1:N,(abs(fdout)));
shading interp;
xlabel('Normalized Frequency (xpi rad/sample)');
ylabel('Delay Index');
end


The point of doppler shift would be a change in frequency. So if you compare to a fixed frequency, that's not going to work.

I assume you try to ,,correct'' the doppler shift so that you can compare with the original frequency as well. That depends on the circumstances. If you just strike samples, you will introduce aliases. These might interfere with what you try to measure. In the worst case, they could introduce ,,beating'' and ruin your signal.

I would rather try to measure different bins, or maybe use a window that is wider than the range of the doppler shift. It depends on the parameters if latter would be feasible. Given that satellites move at a tiny fraction of the speed of light, and that most radio analysis mixes down to MHz range, your doppler shift might be very small (measured in bins) actually.

• I want to remove the PN spreading code from the data. The normal range would be +/-5Khz from center. I need to keep track of the adjustment wherever its made so I can log the number of cycles of shift of the carrier. Jan 7, 2014 at 1:03
• The doppler shift is not reduced by mixing the RF down to IF, if that is what is being suggested above.
– John
Jan 7, 2014 at 12:29

Check this article:
http://asp.eurasipjournals.com/content/2014/1/143

We used energy detector for fft-based acquisition in GNSS systems but you can use it in similar situation.