Changing the amount of points changes the result of the fft

Question

I wrote a script for generating a function (sech(x)), generating the fft, and afterwards measuring the FWTH of the result. Now the strange thing is that when I change the amount of points I use for generating the function, the result of the fft and the FWTH-measurement changes also quite a lot. Is that a result of increased accuracy, or something else? The script I used is:

from __future__ import division
import numpy as np
from math import *
import numpy.fft as fft
from matplotlib import pyplot as plt


def formula(t, fwhm = 0):
    if fwhm == 0:
        t0 = 68e-15#15#68e-15#*1.76#25e-15#68e-15
    else:
        t0 = fwhm/(2*log(2+sqrt(3)))
    #print t0*84e-9
    E0 = 0.68e-9
    p = 1e10
    #if E0 != 0:
#        p = E0/t0
    print("p is: " + str(p))
    inifield = []
    inifield = np.sqrt(p) * 1/np.cosh(t/t0)
    return inifield

def find_nearest(array, values):
    values = np.atleast_1d(values)
    indices = np.abs(np.subtract.outer(array, values)).argmin(0)
    out = array[indices]
    return out if len(out) > 1 else out[0]

def get_FWHM(data, time, inverse):
    epsilon = 1e-12
    max_val = max(data)
    h_max = max_val/2
    idx_max = np.where(data == max_val)
    #idx_h = find_nearest(data, h_max)
    idx_h = []
    print "H_max/2 = " + str(h_max)
    print "Highest values is: " + str(find_nearest(data, h_max * 2))
    subdata = np.split(data, 2)
    for elem in subdata:
        idx_h.append(np.where(elem == find_nearest(elem, h_max)))
    for i, elem in enumerate(idx_h):
        idx_h[i] = int(elem[0])
    idx_h[1] += len(data)/2
    print idx_max
    print "idx_max = " + str(data[int(idx_max[0][0])])
    for elem in idx_h:
        print "idx_h = " + str(elem) + ' ' + str(data[elem])
    if inverse:
        print "Frequency is: " + str(1/time[idx_h[0]]) + " " + str(1/time[idx_h[1]]) + " " + str((1/time[idx_h[1]] - 1/time[idx_h[0]]))
    else:
        print "Time is: " + str(time[idx_h[0]]) + " " + str(time[idx_h[1]]) + " " + str((time[idx_h[1]] - time[idx_h[0]]))

Ntime = 2**17
Tmax = Ntime/2
dt = 2*Tmax/(Ntime-1)
c = 2.99792458e8
t = np.linspace(-1.633e-12, 1.633e-12, Ntime)
#for i, elem in enumerate(t):
#    t[i] *= dt
#for elem in t:
#    print(str(elem/T0))
#    y.append(1/cosh(elem/T0))
y = formula(t, 68e-15)
#get_FWHM(y, t, False)
h = fft.fftshift(fft.ifft(fft.fftshift(y)))
print "get_FWHM for h"
get_FWHM(h, t, True)
print "Target FWHM is: " + str(c/88e-9)
plt.plot(t, y, 'b:', t, h)
plt.show()

When running it with Ntime=2**13, the output is:

Frequency is: -1.51997624747e+14 1.4331204619e+14 2.95309670937e+14
Target FWHM is: 3.40673247727e+15

When using 2**17 points, the output is:

Frequency is: -2.4322403459e+15 2.29325518327e+15 4.72549552917e+15
Target FWHM is: 3.40673247727e+15

Answer 1

Your issue is in the inverse FFT of the shifted signal.

By applying ifft to your shifted signal, you are interpreting your signal as being in the frequency domain. By adding points, you are effectively increasing your Nyquist frequency. But you hold your shape the "same", smearing it across a much wider frequency range. Higher frequency content makes the ifft peak much narrower (narrow peaks in time/frequency domain means broad peaks in corresponding time/frequency domain).

Ultimately, this is more of a math question than a code question. fftshift/fft/ifft is behaving correctly. I'm not certain exactly what you're trying to accomplish, as you are using a frequency modification function (fftshift) on what appears to be a time-domain signal, then applying ifft to the result, and modifying the result as a frequency signal again.

Run the following code at the bottom of your script to see where the number-of-points dependence appears.

plt.figure()

for Ntime in 2**np.arange(14,18):

    # Ntime = 2**14
    Tmax = Ntime/2
    dt = 2*Tmax/(Ntime-1)
    c = 2.99792458e8
    t = np.linspace(-1.633e-12, 1.633e-12, Ntime)
    #for i, elem in enumerate(t):
    #    t[i] *= dt
    #for elem in t:
    #    print(str(elem/T0))
    #    y.append(1/cosh(elem/T0))
    y = formula(t, 68e-15)
    #get_FWHM(y, t, False)
    h1 = fft.fftshift(y)
    h2 = ifft(h1)
    h3 = fft.fftshift(h2)
    h = h3
    # h = fft.fftshift(fft.ifft(fft.fftshift(y)))
    print "get_FWHM for h"
    get_FWHM(h, t, True)
    print "Target FWHM is: " + str(c/88e-9)

    plt.subplot(4,1,1)
    plt.plot(t, y)
    plt.subplot(4,1,2)
    plt.plot(t, h1)
    plt.subplot(4,1,3)
    plt.plot(t, h2)
    plt.subplot(4,1,4)
    plt.plot(t, h3,label='%e'%Ntime)
    plt.show()

plt.legend()