How to estimate the Doppler shift?

How to estimate the Doppler shift?

How to calculate the radial velocity for "jammer"?

clear variables
close all
c = physconst('LightSpeed');
f1 = 30e6;
f2 = 40e6;
fs = max(f1,f2)*2; 

channel = phased.WidebandFreeSpace(...
    'PropagationSpeed',c,...
    'OperatingFrequency',mean(f1,f2),...
    'SampleRate',f1,'NumSubbands',512);

T = 1/fs*5000;
t = 0:1/fs:T-1/fs;

xd = chirp(t,f1,t(end),f2,'complex');


rpos = [0;0;0];
rvel = [0;0;0];

tpos = [1e2;0;0];
tvelJ = [-3e6;0;0];
tvelS = [0;0;0];


add0 = 4000;
x0 = [xd zeros([1 add0])];
sigpower = pow2db(mean(abs(x0).^2));
SNR = 10000;
x0awgn = awgn(x0,SNR,sigpower);
tadd = [t t(end)+t(2:(add0+1))];

channel.reset
jammer = channel(x0awgn.',rpos,tpos,rvel,tvelJ);
channel.reset
desire = channel(x0awgn.',rpos,tpos,rvel,tvelS);

figure
hold on
n = 2^nextpow2(size(x0.',1));
x0_fft = fft(x0,n)/n;
jammer_fft = fft(jammer,n)/n;
desire_fft = fft(desire,n)/n;
f_ax = fs*(0:(n-1))/n;
x0_fftn = x0_fft/((max(abs(x0_fft))/max(abs(jammer_fft))));
plot(f_ax,abs(x0_fftn))
plot(f_ax,abs(jammer_fft))
plot(f_ax,abs(desire_fft))

hold off
legend ('x0norm','jammer','desire')
xlim([2e7 5e7])

Answer 1

[SOUMNATH PAUL]

Hi,

To my understanding you are trying to estimate the doppler shift and calculate the radial velocity for jammer.

Your script doesn't explicitly calculate the Doppler shift or the radial velocity; it simulates the effect of the jammer's motion on the received signal. To calculate the actual Doppler shift or radial velocity from the simulated data, you would need to compare the peak frequency of the jammer's FFT with the original signal's frequency and then apply the Doppler shift formula.

In the below code I have considered the jammer is moving directly towards or away from the receiver, which allows for a simplified calculation of the radial velocity. The Doppler shift is calculated by finding the peak frequency of the jammer's signal and comparing it with the original signal's frequency.

clear variables
close all

% Constants
c = physconst('LightSpeed'); % Speed of light in m/s
f1 = 30e6; % Start frequency of the chirp signal in Hz
f2 = 40e6; % End frequency of the chirp signal in Hz
operatingFrequency = mean([f1, f2]); % Operating frequency in Hz

% The SampleRate should be less than twice the OperatingFrequency
fs = 2 * operatingFrequency * 0.99; % Set SampleRate to slightly less than twice the OperatingFrequency

% Create a wideband free space channel
channel = phased.WidebandFreeSpace(...
    'PropagationSpeed',c,...
    'OperatingFrequency',operatingFrequency,...
    'SampleRate',fs,'NumSubbands',512);

% Time vector for the chirp signal
T = 1/fs*5000;
t = 0:1/fs:T-1/fs;

% Generate a complex chirp signal
xd = chirp(t,f1,t(end),f2,'complex');

% Positions and velocities
rpos = [0;0;0]; % Receiver position
rvel = [0;0;0]; % Receiver velocity
tpos = [1e2;0;0]; % Transmitter position
tvelJ = [-3e6;0;0]; % Jammer velocity
tvelS = [0;0;0]; % Desired signal velocity

% Add zeros to the signal to analyze it after the chirp ends
add0 = 4000;
x0 = [xd zeros([1 add0])];

% Add noise to the signal
sigpower = pow2db(mean(abs(x0).^2));
SNR = 10000;
x0awgn = awgn(x0,SNR,sigpower);

% Time vector including the added zeros
tadd = [t t(end)+t(2:(add0+1))];

% Reset the channel and simulate the signal from the jammer and the desired target
channel.reset
jammer = channel(x0awgn.',rpos,tpos,rvel,tvelJ);
channel.reset
desire = channel(x0awgn.',rpos,tpos,rvel,tvelS);

% Perform FFT on the signals
n = 2^nextpow2(size(x0.',1));
x0_fft = fft(x0,n)/n;
jammer_fft = fft(jammer,n)/n;
desire_fft = fft(desire,n)/n;

% Frequency axis for plotting
f_ax = fs*(0:(n-1))/n;

% Normalize the FFT of the original signal for comparison
x0_fftn = x0_fft/((max(abs(x0_fft))/max(abs(jammer_fft))));

% Plot the FFTs of the signals
figure
hold on
plot(f_ax,abs(x0_fftn))
plot(f_ax,abs(jammer_fft))
plot(f_ax,abs(desire_fft))
hold off
legend('x0norm','jammer','desire')
xlim([2e7 5e7])

% Calculate the Doppler shift for the jammer
[~, idx_jammer] = max(abs(jammer_fft)); % Find index of peak frequency of jammer's signal
f_jammer_peak = f_ax(idx_jammer); % Peak frequency of jammer's signal
f_original_peak = mean([f1, f2]); % Original frequency of the chirp signal

% Assuming the source is stationary, calculate the radial velocity
v_r = c * (f_jammer_peak - f_original_peak) / f_original_peak;

% Display the radial velocity
disp(['Radial velocity of the jammer: ' num2str(v_r) ' m/s'])

Radial velocity of the jammer: -2309558.3526 m/s

Hope it helps!

Regards,

Soumnath

[reincornator]

This is not a completely correct decision. I am using a broadband signal and I want to use all the energy of the signal, not the value of a single frequency.

[SOUMNATH PAUL]

Hi,

Here I have estimated the Doppler shift using the entire energy of the broadband signal.

clear variables
close all

% Constants
c = physconst('LightSpeed'); % Speed of light in m/s
f1 = 30e6; % Start frequency of the chirp signal in Hz
f2 = 40e6; % End frequency of the chirp signal in Hz
operatingFrequency = mean([f1, f2]); % Operating frequency in Hz

% The SampleRate should be less than twice the OperatingFrequency
fs = 2 * operatingFrequency * 0.99; % Set SampleRate to slightly less than twice the OperatingFrequency

% Create a wideband free space channel
channel = phased.WidebandFreeSpace(...
    'PropagationSpeed',c,...
    'OperatingFrequency',operatingFrequency,...
    'SampleRate',fs,'NumSubbands',512);

% Time vector for the chirp signal
T = 1/fs*5000;
t = 0:1/fs:T-1/fs;

% Generate a complex chirp signal
xd = chirp(t,f1,t(end),f2,'complex');

% Positions and velocities
rpos = [0;0;0]; % Receiver position
rvel = [0;0;0]; % Receiver velocity
tpos = [1e2;0;0]; % Transmitter position
tvelJ = [-3e6;0;0]; % Jammer velocity
tvelS = [0;0;0]; % Desired signal velocity

% Add zeros to the signal to analyze it after the chirp ends
add0 = 4000;
x0 = [xd zeros([1 add0])];

% Add noise to the signal
sigpower = pow2db(mean(abs(x0).^2));
SNR = 10000;
x0awgn = awgn(x0,SNR,sigpower);

% Time vector including the added zeros
tadd = [t t(end)+t(2:(add0+1))];

% Reset the channel and simulate the signal from the jammer and the desired target
channel.reset
jammer = channel(x0awgn.',rpos,tpos,rvel,tvelJ);
channel.reset
desire = channel(x0awgn.',rpos,tpos,rvel,tvelS);

% Calculate the cross-correlation between the transmitted and received signals
[correlation, lags] = xcorr(jammer, x0awgn);

% Find the lag at which the correlation is maximized
[~, max_idx] = max(abs(correlation));
time_delay = lags(max_idx) / fs;

% Estimate the Doppler shift using the time delay
doppler_shift = -time_delay * (f2 - f1) / T;

% Calculate the radial velocity of the jammer
radial_velocity = c * doppler_shift / operatingFrequency;

% Display the radial velocity
disp(['Radial velocity of the jammer: ' num2str(radial_velocity) ' m/s'])

Radial velocity of the jammer: 256964.964 m/s

% Plot the cross-correlation result
figure
plot(lags/fs, abs(correlation))
title('Cross-Correlation between Transmitted and Received Signals')
xlabel('Time delay (s)')
ylabel('Correlation magnitude')
grid on

% Indicate the peak correlation point on the plot
hold on
plot(time_delay, abs(correlation(max_idx)), 'ro', 'MarkerFaceColor', 'r')
legend('Cross-Correlation', 'Peak Correlation')
hold off