0
$\begingroup$

I have measurements of a transfer function $H(f)$ in the frequency domain between 1-6 GHz. I also have an input signal, which is pulsed chirp $x(t)$. The goal is to build the time-domain response of the system. To that end, I take a Fourier transform of the initial signal which gives me $X(f)$. Then I build the $H(f)$ by adding it to its conjugate mirror vector (in Python). Then I multiply it to $X(f)$ and then I take the Inverse Fourier Transform to find the time-domain response $y(t)$. First, I'd like to know if this process is correct. We had a discussion that we might just need to consider the positive frequencies by zeroing out the negative ones right before the IFFT. Which one of these two are technically correct?

Secondly, when I plot the real and imaginary parts of $y(t)$, the shapes make sense but there is an imaginary part about one order of of magnitude smaller than the real part. Shouldn't the imaginary part be much smaller than that?

Third, when I take the Fourier transform of $H(f)$, I find that there is an imaginary part. Is that Ok?

Finally, what is the significance of the small rise at the very end?

enter image description here

enter image description here

enter image description here

Here $g_21(t)$ is indeed $h(t)$ and $G_21(f)$ is $H(f)$ in the above description. $y_2(t)$ is the time domain response of my system. I changed the dataset because the imaginary part is more pronounced here.

Here's my entire code.

from scipy.signal import chirp import numpy as np import matplotlib.pyplot as plt N = 16000; sig_duration = 320e-9; dt = sig_duration/N; fs = 1/dt; density = 320; t = np.linspace(0, sig_duration, N); x = chirp(t, 1.5e9, sig_duration, 5.5e9, method='linear', phi=-90) x[:int(N/4)] = 0; x[int(3*N/4):] = 0;

import pandas as pd
import numpy as np

data = "/content/Meas S21 Complex_DARPA Cyl_EXT Excitn Apert 1_0 Absorbers.xlsx"
real_impulse_response = pd.read_excel(data, sheet_name="S21 Real", usecols="B", skiprows=5, nrows=1601, header=None)
imag_impulse_response = pd.read_excel(data, sheet_name="S21 Imag", usecols="B", skiprows=5, nrows=1601, header=None)
impulse_response = np.vectorize(complex)(real_impulse_response.values, imag_impulse_response.values)
experimental_frequencies = pd.read_excel(data, sheet_name="S21 Real", usecols="A", skiprows=5, nrows=1601, header=None)

plt.plot(experimental_frequencies, abs(impulse_response))
plt.yscale("log") 
plt.xlabel("Freq in Hz")
plt.ylabel("$G_{21}(f)$")

from numpy.fft.helper import fftshift
f = np.fft.fftfreq(N, dt)
G21 = np.zeros(len(f), dtype=np.complex_)
impulse_response = impulse_response.flatten()
G21[int(len(f)/2)+density: int(len(f)/2) +density+len(impulse_response)] = impulse_response
#G21_ds = G21 + np.conj(G21[::-1])
#G21_ds = np.fft.fftshift(G21_ds)
G21_ds = G21 + np.conj(G21[::-1])
G21_ds = np.fft.fftshift(G21_ds)
G21_ds[int(N/2):] = 0
G21_sym = (G21_ds + np.conj(G21_ds[::-1]))/2
g = np.fft.ifft(G21_sym)*N

plt.plot(f/1e9, abs(np.fft.fft(x, N)/N))
plt.yscale("log")
plt.xlim(0, 6)
plt.ylim(bottom=1e-4)
plt.xlabel("Freq(GHz)")
plt.ylabel("$X_1(f)$")

#G21_ds[int(N/2):] = 0
#g =np.fft.ifft(G21_ds)*N
plt.subplot(211)
plt.plot(1e9*t, g.real)
plt.xlabel("time(ns)")
plt.ylabel("$Re [g_{21}(t)]$")
plt.subplot(212)
plt.plot(1e9*t, g.imag)
plt.ylabel("$Im[g_{21}(t)]$")
plt.xlabel("time(ns)")
#plt.legend(bbox_to_anchor=(1.0, 1))
plt.show()

y2 = np.fft.ifft(G21_sym * np.fft.fft(x, N))
#plt.axvline(x = 0.8e2, color = 'r', label = 'Chirp Starts')
#plt.axvline(x = 2.4e2, color = 'g', label = 'Chirp Ends')
plt.subplot(311)
plt.plot(1e9*t, x)
plt.xlim([0,320])
plt.ylabel("$x(t)$")
plt.xlabel("time(ns)")

plt.subplot(312)
plt.plot(1e9*t, y2.real, label = "Re(y_2(t))")
plt.ylabel("$REAL~ y_2(t)$")
plt.xlabel("time(ns)")

plt.xlim([0,320])

plt.subplot(313)
plt.plot(1e9*t, y2.imag)
plt.ylabel("$IMAG ~ y_2(t)$")
plt.xlabel("time(ns)")

plt.xlim([0,320])
plt.legend(bbox_to_anchor=(1, 1))
plt.show()
$\endgroup$
4
  • $\begingroup$ Reza have you seen this post which also dealt with eliminating complex values by making H(f) properly from VNA data? dsp.stackexchange.com/questions/86900/… $\endgroup$ Commented Apr 10, 2023 at 11:50
  • $\begingroup$ And this may be of interest: optimized windowing for FFT of chirp for channel estimation: dsp.stackexchange.com/questions/66541/… $\endgroup$ Commented Apr 10, 2023 at 11:52
  • $\begingroup$ Will take a look Dan Boschen. $\endgroup$
    – Reza Afra
    Commented Apr 11, 2023 at 3:32
  • $\begingroup$ Dan! your first post helped yo get rid of the imaginary part. BIG THANKS 😊 $\endgroup$
    – Reza Afra
    Commented Apr 14, 2023 at 2:48

1 Answer 1

2
$\begingroup$

I have measurements of a transfer function $H(f)$ in the frequency domain between 1-6 GHz.

To do an inverse Fourier Transform you will need the entire range from 0Hz to the Nyquist Frequency. One option is to zero out the regions where you have no data but that can generate a lot of time domain ringing and potentially time domain aliasing.

Then I build the $H(f)$ by adding it to its conjugate mirror vector (in Python).

It would be helpful to post your code. This is not a trivial operation.

First, I'd like to know if this process is correct.

Depends on the details. Multiplication in frequency is equivalent to circular (not linear) convolution in time. All signals need to be properly zero padded to avoid time domain aliasing.

We had a discussion that we might just need to consider the positive frequencies by zeroing out the negative ones right before the IFFT. Which one of these two are technically correct?

You need both negative and positive frequency. If you zero out the negative frequency your output signal will be complex (and analytic).

Secondly, when I plot the real and imaginary parts of $y(t)$, the shapes make sense but there is an imaginary part about one order of of magnitude smaller than the real part. Shouldn't the imaginary part be much smaller than that?

If your input signal is real and $H(f)$ is supposed to represent a real world system, then your output signal should be real too. The imaginary part should be 0 or just some numerical noise.

Third, when I take the Fourier transform of $H(f)$, I find that there is an imaginary part. Is that Ok?

I assume you mean "inverse" Fourier transform of $H(f)$? No. That's not okay. Chances are you are doing something wrong here.

Finally, what is the significance of the small rise at the very end?

Most likely that's time domain aliasing or a noise problem.

$\endgroup$
4
  • $\begingroup$ "The imaginary part should be 0" -- well, it should. But typically when you do a DFT numerically, and especially an FFT, the resulting imaginary part will be non-zero because of numerical rounding. You need to purposely ignore (usually after checking that it's small enough to really be just rounding error) the imaginary part. Note that for the purposes of this discussion, an imaginary part 1/10 the amplitude of the real part definitely isn't just rounding error. $\endgroup$
    – TimWescott
    Commented Apr 9, 2023 at 20:05
  • 1
    $\begingroup$ Understood but in this case the imaginary part seems to be roughly in the same order of magnitude than the real part so something is definitely wrong here. $\endgroup$
    – Hilmar
    Commented Apr 10, 2023 at 12:17
  • $\begingroup$ TimWescott I understand that we will have some residual imaginary values but I have posted a new dataset and as you can the complex part is not residual here. In matlab that residual part can be zeroed by using the "symmetric" flag. $\endgroup$
    – Reza Afra
    Commented Apr 11, 2023 at 3:31
  • $\begingroup$ @Hilmar: I was nit picking in anticipation of some DSP newbie reading that out of context. $\endgroup$
    – TimWescott
    Commented Apr 11, 2023 at 18:35

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.