1
$\begingroup$

I'm trying to compute the Fourier Transform of a lightning impulse waveform. The python code I'm using is below. The problem that I am running into is that the magnitude of the DFT and the magnitude of the analytically derived Fourier Transform are off by a factor equal to 1 / time period, where the time period is the length of time that I am using to model the signal for the DFT.

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


def Impulse(t, V_0=1.04, alpha=0.0146e6, beta=2.467e6, t_chop=float("inf")):
    if t<=t_chop:
        return V_0 * (exp(-alpha*t)-exp(-beta*t))
    else:
        return 0


npImpulse = np.vectorize(Impulse, excluded=['V_0', 'alpha', 'beta', 't_chop'])

# Sampling rate and time vector
total_time = 500e-6 # seconds
n_samples = 2**24
sampling_rate = n_samples/total_time # Hz (samples/s)

t_values = np.linspace(0, total_time, n_samples)

# Time constants
tau1 = 68.2e-6
tau2 = 0.405e-6

y_values = npImpulse(t_values, V_0=1.0, alpha=1/tau1, beta=1/tau2)

# plot with various axes scales
plt.figure(1)
plt.plot(t_values, y_values)
plt.grid(True)
plt.show()

# FFT calculation
sp = np.fft.rfft(y_values) # real FFT, + freq only
freq = np.fft.rfftfreq(y_values.size, t_values[1]-t_values[0])

# Normalization
# This is where things are not right.  Should be dividing by n_samples, not samples/sec
sp = sp / n_samples

# Calculation of Fourier Transform analytically
sp_analytic_real = tau1/(1+(2*pi*freq*tau1)**2)-tau2/(1+(2*pi*freq*tau2)**2)
sp_analytic_imag = -(2*pi*freq*(tau1**2)/(1+(2*pi*freq*tau1)**2)-2*pi*freq*(tau2**2)/(1+(2*pi*freq*tau2)**2))
sp_analytic_mag = np.sqrt(sp_analytic_real**2 + sp_analytic_imag**2)

plt.figure(2)
plt.plot(freq, np.abs(sp), label='FFT')
plt.plot(freq, sp_analytic_mag, label='Analytic')
plt.xscale("log")
plt.yscale("log")
plt.legend()
plt.grid(True)
plt.show()

If I set my calculation time period to 1s, the two curves match. If I divide the magnitude of the DFT return by rfft by the sampling rate (or alternately multiply by the 1 / time period), they match. I can't for the life of me figure out why that should be. Can someone clue me in to the obvious that I am missing?

$\endgroup$

1 Answer 1

1
$\begingroup$

That's the case. Assuming you have a band-limited continuous-time signal $x_c(t)$, then the relation between its (analytically derived) CTFT, $X_c(j\Omega)$, and the (numerically computed) DTFT, $X_d(e^{j\omega})$, (or actually its DFT $X_d[k]$ ) of the samples, $x_d[n] = x_c(n T_s)$, of it taken at a period $T_s$ will be :

$$ X_d(e^{j\omega}) = \frac{1}{T_s} X_c( \frac{\omega}{T_s}) ~~~, |\omega| < \pi $$

In practice you are computing the N-point DFT $X[k]$ of $x_d[n]$ which is the uniform samples of $X_d(e^{j \omega})|_{\omega = \frac{2\pi}{N} k}$, for $k=0,1,...,N-1$. And they have the same magnitude relation too.

$$ X[k] =X_d(e^{j\frac{2\pi}{N} k}) = \frac{1}{T_s} X_c( \frac{\omega}{T_s}) ~~~, |\omega| < \pi $$

In other words, you will have a weight of $$ K = \frac{1}{T_s}$$

between the computed DFT / DTFT and the theoretical CTFT of the signal.

$\endgroup$

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.