2
$\begingroup$

I am currently working on a FMCW radar using a SDR I have done the simulations in python and am moving towards simulations on GNU radio. In python it was easy to generate a chirp, just use the scipy.signal.chirp function and give the necessary parameters as input. But in GNU radio while creating my own block I was told that I cannot use the above function. So while trying to find the source code I came to know that scipy.signal.chirp was actually written in c and python bindings were used to run the function in python, So i tried to write the code myself but I am not able to get the correct output. Below is the link for the code I have written. This code is supposed to be run on PLUTO SDR.

import numpy as np
import matplotlib.pyplot as plt
from math import pi

bandwidth = 160e6
#samp_freq = 352e6
chirp_duration = 0.000000002  #ns0.000000002
target_range = 100
rcs = 0.0377
f0 = 76.92e9
f1 = 77.08e9
c = 3e8

slope = bandwidth/chirp_duration
time_delay = 2*target_range/c
#no_samp = samp_freq * chirp_duration
no_samp = 16384
change_in_frequency = slope * time_delay
t1 = chirp_duration

t_axis = np.linspace(0, chirp_duration, int(no_samp))
##change would be done in instantaneous freq addition with change in 
frequency.
##t_axis will be added with time_delay then plottin will be done

def chirp(t, f0, t1, f1, method='linear', phi=0, vertex_zero=True):
    # _chirp_phase implementation
    t = np.asarray(t)
    f0 = float(f0)
    t1 = float(t1)
    f1 = float(f1)
    if method in ['linear', 'lin', 'li']:
        beta = (f1 - f0) / t1
        phase = 2 * np.pi * (f0 * t + 0.5 * beta * t * t)
    else:
       raise ValueError("Invalid method specified.")

# Convert phi to radians
phi *= np.pi / 180

# Generate the chirp signal
signal = np.cos(phase + phi)

return signal

chirp_signal = chirp(t = t_axis , f0 = f0, t1 = chirp_duration, f1 = f1, 
method='linear')

plt.figure(figsize=(12, 3))
plt.xlabel("Time")
plt.ylabel("Amplitude")
plt.title('chirp')
plt.plot(t_axis, chirp_signal)

This is an example of what I should get in the output

This is the output I am getting with my code

$\endgroup$
5
  • $\begingroup$ Can you share your thoughts on why do you think your results are wrong? What did you expect to see and how that differs from what you actually got? $\endgroup$
    – ZaellixA
    Jul 21, 2023 at 9:31
  • $\begingroup$ @ZaellixA Instead of getting a chirp which starts at a lower frequency and ends at a lower frequency, I am getting a signal with constant frequency $\endgroup$
    – Raj Patil
    Jul 21, 2023 at 16:23
  • $\begingroup$ Are you sure? I mean, is your time axis adequate to spot the frequency difference by eye? $\endgroup$
    – ZaellixA
    Jul 21, 2023 at 17:07
  • $\begingroup$ @ZaellixA I even tried plotting a lesser number of samples to check the plot that also showed a constant freq $\endgroup$
    – Raj Patil
    Jul 21, 2023 at 18:08
  • $\begingroup$ Post the applicable c code. And just the applicable c code. $\endgroup$ Jul 21, 2023 at 18:39

1 Answer 1

2
$\begingroup$

Your code works correctly. I have validated it against MATLAB's implementation and the results are the same down to numerical errors. You can see the results in the following figure.

Chirps and their difference

The only difference here is some parameter values. These are $f_{0} = 10^{9} ~ Hz$, $f_{1} = 9 \cdot 10^{9} ~ Hz$ and the total duration $t_{dur} = 5 \cdot 10^{-9} ~ s$ from $2 \cdot 10^{-9} ~ s$ which was in your test (I used this to make the result more obvious).

As already mentioned in a comment, the only problem was your values which, for the set frequency difference did not result in a discernible difference in the time domain representation by eye.

Other than that, the algorithm I used was a copy-pasted version of yours in MATLAB (I had to get rid of the numpy stuff and convert the code for element-wise multiplication but this is more technical than algorithmic).


Update

After some more testing, I found out that the errors accumulate with time in, what seems to be (haven't verified that mathematically so take it with a grain of salt) an unbounded way. This is most probably attributed to the additions you have in your function for the instantaneous phase calculations but I haven't verified that (either).

So, unless you don't intend to use long sweeps or you don't care that much about "small" errors (the order of magnitude $\approx 7.5 \cdot 10^{-6}$ for a $1 ~ s$ sweep) you could very well keep your function, otherwise I suggest you use one from a dedicated library (optimised to avoid such issues).


Code update

The OP has asked for the code changes so below is the code I have used to create the plots of the figure above.

bandwidth = 160e6;
chirp_duration = 0.000000005;  % Old value was 0.000000002
target_range = 100;
rcs = 0.0377;
f0 = 1e9; % Old value was 76.92e9
f1 = 9e9; % Old value was 77.08e9
c = 3e8;

slope = bandwidth/chirp_duration;
time_delay = 2*target_range/c;
no_samp = 16384;
change_in_frequency = slope * time_delay;
t1 = chirp_duration;

t_axis = linspace(0, chirp_duration, no_samp);

chirp_signal = customChirp(t_axis, f0, chirp_duration, f1, 0); % Removed not-needed parameters
matlab_signal = chirp(t_axis, f0, chirp_duration, f1);

% Plot
figure()
subplot(2, 1, 1)
plot(t_axis, chirp_signal, '-', 'LineWidth', 3); hold on;
plot(t_axis, matlab_signal, '-.', 'LineWidth', 3); hold off; grid minor;
ax = gca; ax.FontSize = 16;
xlabel("Time (s)", 'FontSize', 24, 'Interpreter', 'latex')
ylabel("Amplitude (a.u.)", 'FontSize', 24, 'Interpreter', 'latex')
title("Linear chirp signals (time domain)", 'FontSize', 32, 'Interpreter', 'latex')
legend({"Custom", "MATLAB "}, 'FontSize', 16, 'Interpreter', 'latex')

subplot(2, 1, 2)
plot(t_axis, matlab_signal - chirp_signal, 'LineWidth', 3)
ax = gca; ax.FontSize = 16;
xlabel("Time (s)", 'FontSize', 24, 'Interpreter', 'latex')
ylabel("Error", 'FontSize', 24, 'Interpreter', 'latex')
title("Difference between custom and MATLAB implementation", 'FontSize', 32, 'Interpreter', 'latex')


%% Here goes the custom chirp calculation function
function signal = customChirp(t, f0, t1, f1, phi)
    beta = (f1 - f0) / t1;
    phase = 2 * pi * (f0 .* t + 0.5 .* beta .* t .* t);

    phi = phi * pi / 180; % Or use rad2deg()

    signal = cos(phase + phi);
end

Please note that this is MATLAB code. The OP's solution was coded in Python and below are the changes in OP's code to match the code I have used to generate the plots. I haven't tested the code but since I do believe that the numpy functions used do pretty much the same as their MATLAB counterparts I believe you should get the same results (different plots though, unless you add some titles, labels, and MATLAB's implementation on top of yours).

import numpy as np
import matplotlib.pyplot as plt
from math import pi

bandwidth = 160e6
#samp_freq = 352e6
chirp_duration = 0.000000005 # THIS WAS 0.000000002
target_range = 100
rcs = 0.0377
f0 = 1e9 # THIS WAS 76.92e9
f1 = 9e9 # THIS WAS 77.08e9
c = 3e8

slope = bandwidth/chirp_duration
time_delay = 2*target_range/c
#no_samp = samp_freq * chirp_duration
no_samp = 16384
change_in_frequency = slope * time_delay
t1 = chirp_duration

t_axis = np.linspace(0, chirp_duration, int(no_samp))
##change would be done in instantaneous freq addition with change in 
frequency.
##t_axis will be added with time_delay then plottin will be done

def chirp(t, f0, t1, f1, method='linear', phi=0, vertex_zero=True):
    # _chirp_phase implementation
    t = np.asarray(t)
    f0 = float(f0)
    t1 = float(t1)
    f1 = float(f1)
    if method in ['linear', 'lin', 'li']:
        beta = (f1 - f0) / t1
        phase = 2 * np.pi * (f0 * t + 0.5 * beta * t * t)
    else:
       raise ValueError("Invalid method specified.")

    # Convert phi to radians
    phi *= np.pi / 180

    # Generate the chirp signal
    signal = np.cos(phase + phi)

    return signal

chirp_signal = chirp(t = t_axis , f0 = f0, t1 = chirp_duration, f1 = f1, 
method='linear')

plt.figure(figsize=(12, 3))
plt.xlabel("Time")
plt.ylabel("Amplitude")
plt.title('chirp')
plt.plot(t_axis, chirp_signal)
$\endgroup$
8
  • $\begingroup$ I am sorry I could not understand the changes I am supposed to do can you please share .py file or something. $\endgroup$
    – Raj Patil
    Jul 22, 2023 at 8:53
  • 1
    $\begingroup$ @RajPatil here you go. I included the Python code I copied and pasted from your question above with the different numbers. But, as I already said, all you had to do was to change the parameters of your problem (see the answer again, I have even labelled the parameters with the same names you use in your code) and not your algorithm. Have another look and let me know what is confusing you (if anything) $\endgroup$
    – ZaellixA
    Jul 22, 2023 at 9:35
  • $\begingroup$ Thankyou for the code now I understand that the time axis in not adequate to view the change in frequency. But I was using the specific frequency of 77GHz as I was going to actually implement this FMCW radar. So I wanted to ask you that for 77ghz frequency should I increase the chirp duration or what parameter should I change? $\endgroup$
    – Raj Patil
    Jul 22, 2023 at 17:09
  • $\begingroup$ I am not sure I understand your question. 77 GHz is just one frequency, it's not a range of frequencies. Furthermore, the problem was not the duration of the signal but the frequencies. They lie very close to each other to be discernible in your plots. You can increase the duration as much as you want but with the specified parameters your plot will always look like having constant frequency since the change in frequency (for the specified parameters) is not easy to see. But, as it is showcased in the answer, your code is working properly, so you'll get the results you want. $\endgroup$
    – ZaellixA
    Jul 23, 2023 at 2:04
  • $\begingroup$ What is also quite confusing in your code is that you define a whole bunch of variables that do absolutely nothing. They may be useful in other places in your codebase (which we haven't really seen) but they do nothing here. The only variables that actually affect the results are the chirp duration and the two frequencies (initial and final). $\endgroup$
    – ZaellixA
    Jul 23, 2023 at 2:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.