Frequency (Doppler) shift problem

Started by June 20, 2018
Hi all,

I am processing satellite transmissions and assume no prior knowledge of the Doppler
shift. Using a brute force technique I have been able to extract coarse Doppler and
fit a polynomial. I wish to be able to use the polynomial to remove the Doppler
components from the signal but have come up against a problem.

Below is some (modified) code which illustrates the issue I'm having. 
I have created a sine wave and I'm trying to frequency shift it over a range of
frequencies which may not make much sense in this context but I need the process to
work for a time varying Dopplered signal.
In this example I have specified a frequency shift from 200Hz to 300Hz but the
result comes out as a chirp from 200Hz to 400Hz.
I can try a range of frequencies, the lower frequency is always right but the upper
frequency is  usually twice what I specify.

I know I'm not doing this in the most efficient way but it is my best effort to date
so I would love if someone could shed some light on where this method is failing
me.

Also the sine below is not complex as the method seems to completely fail if I make
it so.. My satellite transmissions have been complex sampled!

Thanks in advance
Enda



Fs = 4e6;
deltaT = 1 ./ Fs;
 
    sine1 = dsp.SineWave(1, 2e6);
    sine1.SamplesPerFrame = 4e6;
    sine1.SampleRate = Fs;
    sine1.Amplitude = 1;  
    sine1.ComplexOutput = false;     % Needs to be complex to supress negative
frequencies
    sine1.PhaseOffset = [0];
    y = sine1()';
 
SeqLen = length(y);
    
tme = SeqLen/Fs;        % Time span Seq covers
bit_tme = tme/SeqLen;
Seq_tme = (bit_tme : bit_tme: tme);
 
      doppler = linspace(200, 300, SeqLen);
 
    for i = 1:SeqLen
        phase(i) = 2*pi*doppler(i)*deltaT*(i-1);
        y_dopplered(i) = y(i)*(cos(phase(i))+1i*sin(phase(i))); 
    end  
    
 
NFFT = 2 .^ ceil(log2( length(y) ));
X1 = fftshift(fft(y,NFFT)); %compute FFT     
y_f = (Fs*(-NFFT/2:NFFT/2-1))/NFFT;   %FFT Sample points    
figure; plot(y_f, 20*log10(abs(X1)/max(abs(X1)))); grid on; hold on;
 
NFFT = 2 .^ ceil(log2( length(y_dopplered) ));
X2 = fftshift(fft(y_dopplered,NFFT)); %compute FFT   
f_Filtered_Y3 = (Fs*(-NFFT/2:NFFT/2-1))/NFFT;   %FFT Sample points  
plot(f_Filtered_Y3, 20*log10(abs(X2)/max(abs(X2)))); grid on;    xlim([-100 500]);
xlabel('Frequency offset (Hz)');
 

On 2018-06-20, endakenna@gmail.com <endakenna@gmail.com> wrote:
> Hi all, > > I am processing satellite transmissions and assume no prior knowledge of the
Doppler shift. Using a brute force technique I have been able to extract coarse Doppler and fit a polynomial. I wish to be able to use the polynomial to remove the Doppler components from the signal but have come up against a problem.
> > Below is some (modified) code which illustrates the issue I'm having. > I have created a sine wave and I'm trying to frequency shift it over a range of
frequencies which may not make much sense in this context but I need the process to work for a time varying Dopplered signal.
> In this example I have specified a frequency shift from 200Hz to 300Hz but the
result comes out as a chirp from 200Hz to 400Hz.
> I can try a range of frequencies, the lower frequency is always right but the
upper frequency is usually twice what I specify.
> > I know I'm not doing this in the most efficient way but it is my best effort to
date so I would love if someone could shed some light on where this method is failing me.
> > Also the sine below is not complex as the method seems to completely fail if I
make it so.. My satellite transmissions have been complex sampled!
> > Thanks in advance > Enda > > > > Fs = 4e6; > deltaT = 1 ./ Fs; > > sine1 = dsp.SineWave(1, 2e6); > sine1.SamplesPerFrame = 4e6; > sine1.SampleRate = Fs; > sine1.Amplitude = 1; > sine1.ComplexOutput = false; % Needs to be complex to supress negative
frequencies
> sine1.PhaseOffset = [0]; > y = sine1()'; > > SeqLen = length(y); > > tme = SeqLen/Fs; % Time span Seq covers > bit_tme = tme/SeqLen; > Seq_tme = (bit_tme : bit_tme: tme); > > doppler = linspace(200, 300, SeqLen); > > for i = 1:SeqLen > phase(i) = 2*pi*doppler(i)*deltaT*(i-1); > y_dopplered(i) = y(i)*(cos(phase(i))+1i*sin(phase(i))); > end > > > NFFT = 2 .^ ceil(log2( length(y) )); > X1 = fftshift(fft(y,NFFT)); %compute FFT > y_f = (Fs*(-NFFT/2:NFFT/2-1))/NFFT; %FFT Sample points > figure; plot(y_f, 20*log10(abs(X1)/max(abs(X1)))); grid on; hold on; > > NFFT = 2 .^ ceil(log2( length(y_dopplered) )); > X2 = fftshift(fft(y_dopplered,NFFT)); %compute FFT > f_Filtered_Y3 = (Fs*(-NFFT/2:NFFT/2-1))/NFFT; %FFT Sample points > plot(f_Filtered_Y3, 20*log10(abs(X2)/max(abs(X2)))); grid on; xlim([-100
500]);
> xlabel('Frequency offset (Hz)'); > >
I have not read your code in detail, but the factor of 2 is likely due to: Vout=sin((w0+kt)*t) Instantaneous frequency = d(phase)/d(time) = d/dt (w0+kt)*t = w0 + 2*k*t Note the factor of two from the derivative.
On Thursday, 21 June 2018 01:38:00 UTC+1, Wendell  wrote:
> On 2018-06-20, endakenna@gmail.com <endakenna@gmail.com> wrote: > > Hi all, > > > > I am processing satellite transmissions and assume no prior knowledge of the
Doppler shift. Using a brute force technique I have been able to extract coarse Doppler and fit a polynomial. I wish to be able to use the polynomial to remove the Doppler components from the signal but have come up against a problem.
> > > > Below is some (modified) code which illustrates the issue I'm having. > > I have created a sine wave and I'm trying to frequency shift it over a range of
frequencies which may not make much sense in this context but I need the process to work for a time varying Dopplered signal.
> > In this example I have specified a frequency shift from 200Hz to 300Hz but the
result comes out as a chirp from 200Hz to 400Hz.
> > I can try a range of frequencies, the lower frequency is always right but the
upper frequency is usually twice what I specify.
> > > > I know I'm not doing this in the most efficient way but it is my best effort to
date so I would love if someone could shed some light on where this method is failing me.
> > > > Also the sine below is not complex as the method seems to completely fail if I
make it so.. My satellite transmissions have been complex sampled!
> > > > Thanks in advance > > Enda > > > > > > > > Fs = 4e6; > > deltaT = 1 ./ Fs; > > > > sine1 = dsp.SineWave(1, 2e6); > > sine1.SamplesPerFrame = 4e6; > > sine1.SampleRate = Fs; > > sine1.Amplitude = 1; > > sine1.ComplexOutput = false; % Needs to be complex to supress negative
frequencies
> > sine1.PhaseOffset = [0]; > > y = sine1()'; > > > > SeqLen = length(y); > > > > tme = SeqLen/Fs; % Time span Seq covers > > bit_tme = tme/SeqLen; > > Seq_tme = (bit_tme : bit_tme: tme); > > > > doppler = linspace(200, 300, SeqLen); > > > > for i = 1:SeqLen > > phase(i) = 2*pi*doppler(i)*deltaT*(i-1); > > y_dopplered(i) = y(i)*(cos(phase(i))+1i*sin(phase(i))); > > end > > > > > > NFFT = 2 .^ ceil(log2( length(y) )); > > X1 = fftshift(fft(y,NFFT)); %compute FFT > > y_f = (Fs*(-NFFT/2:NFFT/2-1))/NFFT; %FFT Sample points > > figure; plot(y_f, 20*log10(abs(X1)/max(abs(X1)))); grid on; hold on; > > > > NFFT = 2 .^ ceil(log2( length(y_dopplered) )); > > X2 = fftshift(fft(y_dopplered,NFFT)); %compute FFT > > f_Filtered_Y3 = (Fs*(-NFFT/2:NFFT/2-1))/NFFT; %FFT Sample points > > plot(f_Filtered_Y3, 20*log10(abs(X2)/max(abs(X2)))); grid on; xlim([-100
500]);
> > xlabel('Frequency offset (Hz)'); > > > > > I have not read your code in detail, but the factor of 2 is > likely due to: > Vout=sin((w0+kt)*t) > Instantaneous frequency = d(phase)/d(time) > = d/dt (w0+kt)*t = w0 + 2*k*t > > Note the factor of two from the derivative.
On Thursday, 21 June 2018 01:38:00 UTC+1, Wendell wrote:
> On 2018-06-20, endakenna@gmail.com <endakenna@gmail.com> wrote: > > Hi all, > > > > I am processing satellite transmissions and assume no prior knowledge of the
Doppler shift. Using a brute force technique I have been able to extract coarse Doppler and fit a polynomial. I wish to be able to use the polynomial to remove the Doppler components from the signal but have come up against a problem.
> > > > Below is some (modified) code which illustrates the issue I'm having. > > I have created a sine wave and I'm trying to frequency shift it over a range of
frequencies which may not make much sense in this context but I need the process to work for a time varying Dopplered signal.
> > In this example I have specified a frequency shift from 200Hz to 300Hz but the
result comes out as a chirp from 200Hz to 400Hz.
> > I can try a range of frequencies, the lower frequency is always right but the
upper frequency is usually twice what I specify.
> > > > I know I'm not doing this in the most efficient way but it is my best effort to
date so I would love if someone could shed some light on where this method is failing me.
> > > > Also the sine below is not complex as the method seems to completely fail if I
make it so.. My satellite transmissions have been complex sampled!
> > > > Thanks in advance > > Enda > > > > > > > > Fs = 4e6; > > deltaT = 1 ./ Fs; > > > > sine1 = dsp.SineWave(1, 2e6); > > sine1.SamplesPerFrame = 4e6; > > sine1.SampleRate = Fs; > > sine1.Amplitude = 1; > > sine1.ComplexOutput = false; % Needs to be complex to supress negative
frequencies
> > sine1.PhaseOffset = [0]; > > y = sine1()'; > > > > SeqLen = length(y); > > > > tme = SeqLen/Fs; % Time span Seq covers > > bit_tme = tme/SeqLen; > > Seq_tme = (bit_tme : bit_tme: tme); > > > > doppler = linspace(200, 300, SeqLen); > > > > for i = 1:SeqLen > > phase(i) = 2*pi*doppler(i)*deltaT*(i-1); > > y_dopplered(i) = y(i)*(cos(phase(i))+1i*sin(phase(i))); > > end > > > > > > NFFT = 2 .^ ceil(log2( length(y) )); > > X1 = fftshift(fft(y,NFFT)); %compute FFT > > y_f = (Fs*(-NFFT/2:NFFT/2-1))/NFFT; %FFT Sample points > > figure; plot(y_f, 20*log10(abs(X1)/max(abs(X1)))); grid on; hold on; > > > > NFFT = 2 .^ ceil(log2( length(y_dopplered) )); > > X2 = fftshift(fft(y_dopplered,NFFT)); %compute FFT > > f_Filtered_Y3 = (Fs*(-NFFT/2:NFFT/2-1))/NFFT; %FFT Sample points > > plot(f_Filtered_Y3, 20*log10(abs(X2)/max(abs(X2)))); grid on; xlim([-100
500]);
> > xlabel('Frequency offset (Hz)'); > > > > > I have not read your code in detail, but the factor of 2 is > likely due to: > Vout=sin((w0+kt)*t) > Instantaneous frequency = d(phase)/d(time) > = d/dt (w0+kt)*t = w0 + 2*k*t > > Note the factor of two from the derivative.
I am seeing something similar when i try this with a sinusoid. Below is my code where I am wanting to sweep a 400MHz tome between -50Mhz and 50Mhz. Meaning the sweep should span between 350MHz and 450MHz. But that is no the case. Don&rsquo;t understand what the reason for this is. import numpy as np import matplotlib.pyplot as plt def gen_sine(freq=380e6, fs=6080e6, n_samp=6400): ts = 1/fs t_arr = np.arange(0, n_samp)*ts #sine_sig = np.exp(2 * 1j * np.pi * freq * t_arr) sine_sig = np.sin(2 * np.pi * freq * t_arr) return sine_sig,ts def freq_shift_sine(sine, ts, shift_freq = 50e6): tx_sig_t = np.linspace(0, ts*sine.shape[-1], num=sine.shape[-1]) #tx_sig_sqrd = np.square(tx_sig_t) #hift the sine freq_shftd_sig = sine * np.exp(1.0j * 2 * np.pi * (shift_freq * tx_sig_t)) #freq_shftd_sig = sine * np.exp(1.0j * np.pi * (shift_freq * tx_sig_sqrd)) return freq_shftd_sig def sweep_sine(sine, ts, up_lim = 50e6, low_lim = -50e6): tx_sig_t = np.arange(0, sine.shape[-1])*ts #tx_sig_sqrd = np.square(tx_sig_t) freq_step_arr = np.linspace(low_lim, up_lim, sine.shape[-1]) dopp_shftd_sig = sine * np.exp(1.0j * 2 * np.pi * (freq_step_arr * tx_sig_t)) #dopp_shftd_sig = sine * np.exp(1.0j * np.pi * (freq_step_arr * tx_sig_sqrd)) return dopp_shftd_sig if __name__=='__main__': #generate a sine wave 16 times over sampled tx_sig, t_samp = gen_sine(freq=400e6, fs=6400e6, n_samp=6400) #do an fft tx_sig_fft = np.fft.fft(tx_sig) #generate freqency axis for fft freq_arr = np.fft.fftfreq(tx_sig.shape[-1], t_samp) #shift sine wave tx_sig_shifted = freq_shift_sine(tx_sig, t_samp, shift_freq = -100e6) #fft the shifted sine tx_sig_shftd_fft = np.fft.fft(tx_sig_shifted) #sweep sine wave by up_lim+low_lim Hz tx_sig_swept = sweep_sine(tx_sig, t_samp, up_lim = 50e6, low_lim = -50e6) #fft the swept sine tx_sig_swept_fft = np.fft.fft(tx_sig_swept) #plt.figure() #plt.plot(freq_arr, abs(tx_sig_swept_fft)) #plot sine wave fft #plt.figure() plt.figure() plt.plot(freq_arr, abs(tx_sig_fft)) plt.plot(freq_arr, abs(tx_sig_shftd_fft)) plt.plot(freq_arr, abs(tx_sig_swept_fft)) plt.axis([0,1e9, 0, 2e3]) plt.figure() plt.specgram(tx_sig_swept, NFFT=80, Fs=6400e6, noverlap=16) #plt.axis([0,0.000001, 0, 5e6]) plt.show()
be aware Doppler "SHIFT" is a misnomer

the effect does NOT create an equal shift to all frequencies.

A more correct term would be Doppler stretch in that higher frequencies are shifted
more compared to lower ones.  

This sometimes doesn't matter for narrow band signals but it can matter for wideband
signals.

m


On Thursday, June 21, 2018 at 7:14:06 PM UTC+1, mako...@yahoo.com wrote:
> be aware Doppler "SHIFT" is a misnomer > > the effect does NOT create an equal shift to all frequencies. > > A more correct term would be Doppler stretch in that higher frequencies are
shifted more compared to lower ones.
> > This sometimes doesn't matter for narrow band signals but it can matter for
wideband signals.
> > m
Thanks but that is actually the point of this function; to 'Doppler stretch' across a range of frequencies!
thanks for your input, it took me a while to understand but I think I have it
working now.
On 21.06.2018 22:43, endakenna@gmail.com wrote:
> On Thursday, June 21, 2018 at 7:14:06 PM UTC+1, mako...@yahoo.com wrote: >> be aware Doppler "SHIFT" is a misnomer >> the effect does NOT create an equal shift to all frequencies. >> >> A more correct term would be Doppler stretch in that higher frequencies are
shifted more compared to lower ones.
>> >> This sometimes doesn't matter for narrow band signals but it can matter for
wideband signals.
>> >> m > > Thanks but that is actually the point of this function; to 'Doppler stretch'
across a range of frequencies!
>
Besides, do you have an idea of how does the Doppler effect work on radio signals? Suppose you've started with a complex-valued baseband signal s(t). And you've up-converted it to the carrier frequency f. After that you've ended up with the signal s^(t) which equals s^(t) = Real(s(t)) * cos(2*pi*f*t) - Imag(s(t)) * sin(2*pi*f*t) Now, the Doppler effect (the part of it that matters) is nothing else but a linear stretch of the time which arises due to Lorentz transformations, when you move to an inertial reference frame that moves with a constant velocity. So, to account for the Doppler effect, we need to stretch the time. Say, we get tn = t*alpha (where alpha is close to unity). So, the received signal is s_r(tn) = Real(s(t*alpha)) * cos(2*pi*f*alpha*t) - Imag(s(t*alpha)) * sin(2*pi*f*alpha*t) The part of the "Doppler stretch" which affects the carrier frequency is usually called the Doppler shift. Just keep in mind that your baseband signal is also affected by the Doppler effect -- which might matter, or not. ^_^ Gene
On 06.07.2018 11:12, Gene Filatov wrote:
> On 21.06.2018 22:43, endakenna@gmail.com wrote: >> On Thursday, June 21, 2018 at 7:14:06 PM UTC+1, mako...@yahoo.com wrote: >>> be aware Doppler "SHIFT" is a misnomer >>> the effect does NOT create an equal shift to all frequencies. >>> >>> A more correct term would be Doppler stretch in that higher >>> frequencies are shifted more compared to lower ones. >>> >>> This sometimes doesn't matter for narrow band signals but it can >>> matter for wideband signals. >>> >>> m >> >> Thanks but that is actually the point of this function; to 'Doppler >> stretch' across a range of frequencies! >> > > Besides, do you have an idea of how does the Doppler effect work on > radio signals? > > Suppose you've started with a complex-valued baseband signal s(t). And > you've up-converted it to the carrier frequency f. After that you've > ended up with the signal s^(t) which equals > > s^(t) = Real(s(t)) * cos(2*pi*f*t) - Imag(s(t)) * sin(2*pi*f*t) > > Now, the Doppler effect (the part of it that matters) is nothing else > but a linear stretch of the time which arises due to Lorentz > transformations, when you move to an inertial reference frame that moves > with a constant velocity. > > So, to account for the Doppler effect, we need to stretch the time. Say, > we get tn = t*alpha (where alpha is close to unity). So, the received > signal is > > s_r(tn) = Real(s(t*alpha)) * cos(2*pi*f*alpha*t) - Imag(s(t*alpha)) * > sin(2*pi*f*alpha*t) > > The part of the "Doppler stretch" which affects the carrier frequency is > usually called the Doppler shift. Just keep in mind that your baseband > signal is also affected by the Doppler effect -- which might matter, or > not. > > ^_^ > > Gene >
Also -- just in case you are building a GPS/GNSS receiver (like, who cares? just trying to be helpful), I strongly recommend Kaplan and Hegarty's "Understanding GPS/GNSS: Principles and Applications". Didn't have the opportunity to read the third edition yet -- but the second was quite good! Gene.