How to correctly scale the FFT python function in order to validate it with the rect() and sinc() function pair?

Question

I am currently working on an python implementation that is using the FFT to convert signals in the time domain to signals in the frequency domain and the other way around. In order to validate my FFT function, I’ve tried to use the rectangular and sinc function pair given by

def rect(t, T): #rectangular function
    ans = np.zeros(len(t))
    for i in range(len(t)):
        if abs(t[i]) < T/2:
            ans[i] = 1
    return ans

and

def Tsinc_hz(f, T): # sinc function
    return np.sin(np.pi * f * T) / (np.pi * f)

Unfortunately, the results are unsatisfying and I don’t quit get why. I hope someone is able to help me.

Here is my FFT code:

def to_frequency_domain(t, ys):
    last = t[-1]  # length of time sample, so is eual to limit
    N = len(t)  # number of samples in the time domain
    T = last/N  # 1/Sample rate
    xf = fftfreq(N, T) 
    xf = fftshift(xf)
    yf = fft(ys, N) 
    yf = fftshift(yf) 
    return xf, yf
 
def to_time_domain(x, y):
    N = len(x)  # Number of samples in the time domain
    last = (len(x)-1) / (2*x[-1]) # length in the time domain 
    T = last/(N-1) # 1/sample rate
    ys =ifft(ifftshift(y), N)   
    return ys

The following code is supposed to validate the FFT implementation by comparing the analytical functions given above with the FFT output in time and frequency domain.

t = np.linspace(0, 10, 250*10) 
ans = rect(t, 1)

x_val,y_val = to_frequency_domain(t, ans)

f = np.linspace(-125, 125, 2500)
proof = Tsinc_hz(f, 1)

# sinc function comparison plot
plt.plot(x_val, y_val.real)
plt.plot(f,proof.real, c='r')
plt.show()

proof_time = to_time_domain(f, proof)

# rectangular function comparison plot
plt.plot(t, proof_time, linewidth = 3, c = 'r')
plt.plot(t, ans)
plt.show()

Running the code gives the following plot: sinc function comparison plot and rectangular function comparison plot

It is obvious that the there is a scaling problem. I know that the frequency domain is dependent on the sample rate of my time domain. In this case I have used a sample rate of 250 and 2500 data points, meaning that I have 0,1 hz per frequency bin and thus my x-axis in the frequency domain reach form - 125 to 125. I was wondering if I can also formulate a relation between the frequency power and the sample rate. I was thinking that if I keep my data points constant and reduce my sampling rate the function is obviously jammed along the x-axis. Can this somehow result in stretching along the y-axis? Furthermore, the transfer of the sinc-function in the second plot (proof_time) is mirrored.

Following, I have tried to fix my first problem by dividing the FFT output by the sampling rate

def to_frequency_domain(t, ys):
    last = t[-1]  # length of time sample, so is eual to limit
    N = len(t)  # number of samples in the time domain
    T = last/N  # 1/Sample rate
    xf = fftfreq(N, T) 
    xf = fftshift(xf)
    yf = fft(ys, N)*T
    yf = fftshift(yf) 
    return xf, yf
 
def to_time_domain(x, y):
    N = len(x)  # Number of samples in the time domain
    last = (len(x)-1) / (2*x[-1]) # length in the time domain 
    T = last/(N-1) # 1/sample rate
    ys =ifft(ifftshift(y)/T, N)  
    return ys

which gives: sinc function comparison plot_scaled, rectangular function comparison plot_scaled

This time both functions are obviously closer together but the result is still not perfect.

Additional, I am a little confused by the output of the FFT when I define my time domain to be symmetrical about zero t =np.linespace(-10,10, 250*20))

sinc function comparison plot

Why is the blue curve doubled like that? Supposedly because we have a negative and positive frequency component for every value in time, right? But how do I fix that?

I have tried figuring it out for a while now but just can’t seem to solve the problem, so I am very grateful for every tip! Thanks in advance!

Answer 1

My guess is that you have to scale by d_omega= d_f/(2*pi) as the fft assumes a time step of 1 by default while you are using a different time step. If i add

df = xf[1] - xf[0]
domega = df / 2 / np.pi
yf = fft(ys, N) * domega

it scales.