# GPS C/A Code Time Delay Using Ultrasound

I am undertaking a project which requires to measure a distance from source to receiver. I am doing this using two ultrasonic transducers operating at 40kHz and using FFT search to determine the time delay and from that the distance. These are the following steps I'm taking to achieve this. All the processing is done in MATLAB.

1. Transmitting a C/A code at 1kHz using a carrier frequency of 40kHz
2. At the receiving end I sample the incoming signal at 160kHz
3. Generate a C/A code locally with a sampling frequency (fs) of 160kHz
4. Multiply the received signal with the one generated in step (3)
5. Multiply the signal from step (4) with sin(2*pifct) and cos(2*pifct), where fc=40kHz and t=0:1/fs:1.023-(1/fs) to obtain the I and Q components down-mixed to 40khz.
6. Then, compute the fft down-mixed signal, using fft(I+iQ), where i=imaginary
7. Compute the fft of the localy generated C/A code and then take the conjugate of that: conj(fft(C/A))
8. Compute the ifft,square and take the absolute magnitude of the result, as so: abs(ifft(step7*step6)).ˆ2 to get the time delat

Are the above steps correct? Is there anything else I have to do or am I missing something? It works fine in matlab but when I use 'real' data the correlation is much weaker. I haven't actually used the ultrasonic yet I'm simply using one microcontroller to generate and send signal and another to sample it (connected with a wire). When I do the processing in matlab the results aren't as good as the pure simulation. Note: DSP is very new to me, I'm just stuck with this project and learning as I go along. C/A Code generator [taken from][1]

tap=[2 6; 3 7; 4 8; 5 9; 1 9; 2 10; 1 8; 2 9; 3 10; 2 3; 3 4;
5 6; 6 7; 7 8; 8 9; 9 10; 1 4; 2 5; 3 6; 4 7; 5 8; 6 9; 1 3;
4 6; 5 7; 6 8; 7 9; 8 10; 1 6; 2 7; 3 8; 4 9; 5 10;4 10;
1 7; 2 8; 4 10];

G1 = -1*ones(1,10);
G2 = G1;

s1 = tap(n,1);
s2 = tap(n,2);

tmp = 0;

for i = 1:1023
G(i) = G2(s1)*G2(s2)*G1(10);
tmp = G1(1);
G1(1) = G1(3)*G1(10);
G1(2:10) = [tmp G1(2:9)];

tmp = G2(1);
G2(1) = G2(2)*G2(3)*G2(6)*G2(8)*G2(9)*G2(10);
G2(2:10) = [tmp G2(2:9)];
end

%resample
x=1:n_samples
CAa(x) = G(ceil(x*chiprate/fs));
end


  function send = s()
clear all;
%        CAa = ca(n,chiprate,fs,n_samples)
CAtest = CAa(1,1e3,160e3,(160e3*1.023)); %generate the outgoing signal
Fs=160e3;
n_data = 160e3*1.023;
t =(0:(n_data-1))/Fs
carrier = cos(2*pi*40e3*t); %carrier @40khz.
send = carrier.*CAtest %modulate

end


Code for xcorrr

close all;
modulated =  send();                          %generate incoming signal
%modulated = circshift(modulated,80000,2);    %simulate a time delay
fc = 40e3;                                    %carrier frequency
Fs = 160e3;
Ts = 1/(Fs);
n_samples = Fs*1.023;                         %sampling period
t = (0:(n_samples-1))/Fs;                     %generate time vector

b = CAa(1,1e3,160e3,(160e3*1.023));           %generate a C/A code locally

rx = modulated .* exp(-j*2*pi*fc.*t);            %downconvert to baseband

corrOut = abs(ifft(fft(rx).*conj(fft(b)))).^2;
figure
plot(corrOut);

• Yes your process is correct. If you would like to see an example of this, I have code here you could follow to compare to yours: mathworks.com/matlabcentral/fileexchange/…. Do you have a carrier offset to deal with due to clock offsets between your transmitter and receiver? Your correlation vs frequency offset would be a sinc function with the first null at 1/T where T is the length of the data you correlate over. – Dan Boschen Apr 14 '17 at 21:57
• Can you post a plot of your actual results for your circular correlation, as well as the standard deviation of your received signal? – Dan Boschen Apr 14 '17 at 22:13
• @DanBoschen Thank you very much for your help. Unfortunately I don't have access to matlab right now. I'll try to get the plots during the weekend if not I'll upload them on Monday. As for the carrier offset I have nothing to deal with it...How should I go about correcting it? I'll also post some raw data from my adc. – smokingRooster Apr 14 '17 at 22:52
• Consider downloading and using Octave as it is free and 99% MATLAB compatible. Look at my answers to other posts as a few deal with carrier recovery – Dan Boschen Apr 14 '17 at 23:02
• I don't understand your use of "square" instead of sine and cosine. -1 square is not 90° shifted so I don't think it is doing what you think. I simply multiply by $exp(-j 2 \pi f_c t)$ to do the quadrature downconversion on one line. Perhaps this is where your issue is. Review your spectrums using FFT's to see if the signal is where you think it is. Also why do you use "Square" for your carrier in your generation instead of cosine? – Dan Boschen Apr 16 '17 at 14:19

Also the -1 square in your receiver will not generate a quadrature signal from + square as you have used, copied below:

This does not produce a sine and cosine:

x2 = square(2*pi*fc*t); %sin
x3 = -1*square(2*pi*fc*t);%cos


Here as well sin and cos directly would be more appropriate, but the following is how I would approach your receiver, using exp instead of sines and cosines since $e^{-j\omega t} = cos(\omega t)-jsin(\omega t)$:

modulated =  send();  %generate incoming signal at carrier fc
rx = modulated * exp(-j*2*pi*fc*t)   %downconvert to baseband
%the above will have a high frequency image at -2fc that must be filtered
coeff= firls(31, [0 .2 .5 1], [1 1 0 0]);
rxfilt= filter(coeff,1,rx)
corrOut = ifft(fft(rxfilt).*conj(fft(CA_local));


The $e^{-j\omega t}$ will shift the positve half of your modulated signal to baseband (0) and the negative half to $-2f_c$, so the example filter shown is to filter out the image. (Note there was no consideration given for actual filter design, this was quickly "pulled out of the hat" just to demonstrate what will work).

• What you are telling me makes sense. I tried your code...but I'm still getting the wrong results. What's even more frustrating is the fact that xcorr(modulated,CA_Local) works perfectly... [link] (imgur.com/a/9LJEe) – smokingRooster Apr 16 '17 at 14:55
• If xcorr(modulated, CA_local) works it means that "modulated" is not at the carrier frequency like you think. You can also check then that corrOut=ifft(fft(modulated).*conj(fft(CA_local)) should also work. I will look at why modulated is not at the carrier.... Note I also updated some signs in what I previously wrote, but what I had before would only cause a spectral inversion which would not matter in this case (BPSK). Can you change your "square" in your modulated to be "cos"? – Dan Boschen Apr 16 '17 at 15:01
• I've made all the changes you suggested, square is now cos. I also double checked my carrier frequency with an fft and it seems correct (40khz) . I'll post the updated codes in the main question – smokingRooster Apr 16 '17 at 15:06
• I looked at the plot and it is not working perfectly in the sense that the signal has peak positve and negative correlations--- if both signals were both at baseband then you would see a positive correlation only. Also the peak correlation should be N times your signal level where N is the total number of samples. I am not sure what your second plot is but shouldn't both signals be the same length based on your code? And isn't your code going +1/-1? (If your code generator is returning 0 1 you should multiply it by 2 and subtract 1 to make it bipolar). If so then you should see 160K for max – Dan Boschen Apr 16 '17 at 15:08
• I suggest doing the following: Make your xcorr work without using a carrier on either side first. This should result in a peak correlation of 160,000 and only positive. Once that works, introduce the carrier on the transmit side and carrier wipe-off on the receiver side as you have done. Without seeing the details of your code generator the problem may be there, perhaps out of sync in what we think the actual sampling rate is--but the above process will clarify that. When you do the correlation, confirm that the corr drops linearly to near zero over the duration of 1 chip in your timebase. – Dan Boschen Apr 16 '17 at 15:09