# Why do problems arise in FFT for smaller value of df in Python?

import numpy as np
import matplotlib.pyplot as plt

Fs = 200                        # sampling rate
Ts = 1.0/Fs                      # sampling interval
t = np.arange(0, 1, Ts)            # time vector

f_f = np.arange(30, 80, 0.001)
a = np.arange(0, len(f_f), 1)
y2 = 0.0
for i in range(len(f_f)):
y2 = y2 +  np.sin(2 * np.pi * (f_f[i]) * t)

n = len(y2)
k = np.arange(n)
T = n/Fs
frq = k/T               # two sides frequency range
freq = frq[range(int(n/2))]     # one side frequency range

Y = np.fft.fft(y2)/n                # fft computing and normalization
#Y = np.fft.fft(y2)
Y = Y[range(int(n/2))]
plt.plot(freq, abs(Y), 'r-')    # Fourier data
plt.show()


Here, as you may have noticed I have chosen df = 0.001 and the corresponding DFT is shown in the attached pic. However, while I am considering df = 1 instead, DFT comes out as a sort of continuous square wave over the frequency domain of the signal (expected since the amplitude of all the waves is equal). Why is it then getting distorted for df = 0.1 or smaller? I mean why two peaks are appearing at f = 30 Hz & f = 80 Hz and not the entire band?

Even, the sampling rate has been so chosen so that Nyquist rate exceeds maximum frequency component of the pulse.

The discrete Fourier transform for a signal of period $$T$$ with $$N$$ samples reads in its inverse or reconstruction form as $$y(t)=\frac1{N}\sum_{k=-N/2}^{N/2}c_k e^{i2\pi k\frac{t}{T}}$$ with redundancy in $$c_N$$ and $$c_{-N}$$ if $$N$$ is even. Sampling this at points $$t_m=\frac{mT}{N}$$ gives a completely determined linear system whose solution is given by the forward DFT, which is the basis for the Nyquist sampling theorem.

The cycle frequency $$f$$ occurs in the exponential as $$e^{i2\pi ft}$$, so that $$f_k=\frac{k}{T}$$ and the frequency spacing is $$df=\frac1T$$.

The sampling theorem in the form used (more Nyquist than Shannon) only applies to periodic functions. As your sampling interval has length $$T=1$$, thus $$df=1$$, this only allows for frequencies that have $$1$$ as one of their periods, that is, for angular frequencies $$2\pi f$$ where the cyclic frequency $$f$$ is an integer.

If you use $$df=1/1000$$ as spacing between frequencies, then the period to be used has length $$T=1000$$ and correspondingly more samples have to be taken for the same sampling frequency.

In both cases you will get the box shape in the amplitude picture, taking a quadratically increasing time to construct the signal from about $$\sim f_s·T^2$$ sine evaluations.

programming details

• floating point arithmetic is unreliable, so in Ts = 1.0/Fs; t = np.arange(0, 1, Ts) it is uncertain what the number of points in t will be in the result. Change that to t = np.linspace(0,1,Fs+1)[:-1] for a deterministic outcome. Note that if $$x_k=kh$$ and $$x_N=T$$, then the fundamental period of $$T$$-periodic signal is sampled with $$x_0,x_1,...,x_{N-1}$$, at $$x_N,...$$ starts a new period that repeats $$x_0,...$$ It is easier to construct the larger array from $$0$$ to $$T$$ and reduce it than to directly construct the shorter array.
• one can shorten the construction of the signal to y2 = sum(np.sin(2 * np.pi * f * t) for f in f_f). This also solves the missing declaration of y2.

The signal construction for the frequency spacing of $$df=0.001$$ is prohibitively expensive, it requires the evaluation of the sine at $$10\,000\,000\,000$$ points, for $$50/df=50\,000$$ frequencies at $$f_s·T=200\,000$$ time points. Waiting 5 to 20 minutes (pure computation cost plus estimated python structures overhead) for a small experiment is too much. Reduced by a factor of 10, this becomes bearable, the signal constructed as

Fs = 200                  # sampling rate
Ts = 100                  # end of sampling interval
Ns = Ts*Fs
dF = 1/Ts
t = np.linspace(0, Ts, Ns+1)[:-1]            # time vector

f_f = np.concatenate([np.arange(1,10,dF),np.arange(30, 80, dF)])

y2 = sum(np.sin(2 * np.pi * f * t) for f in f_f)


gives the expected amplitude plot

One can also see the summation as kind-of Riemann sums and replace it by the exact integral values. As the result is a combination of functions of the type $$\frac{\sin x}{x}$$, one would have to center the time interval at $$t=0$$ so that the symmetry causes the signal to be continuous if cut to the time interval and periodically continued.

With all that, the signal construction code

Fs = 200                  # sampling rate
Ts = 1000                  # end of sampling interval
Ns = Ts*Fs
t = np.linspace(0, Ts, Ns+1)[:-1] + 1e-10 - Ts/2    # zero-centered time vector, avoid zero-division

y2 = (-np.cos(2*np.pi*80*t)+np.cos(2*np.pi*30*t)-np.cos(2*np.pi*10*t)+np.cos(2*np.pi*1*t))/t


is again fast and results in the FFT amplitudes

The Gibbs-oscillations, visible as overshoot at the corners, are typical, they result from the filled-in frequencies of the integration.

• in your last code where you have mentioned about Riemann sums, Why have you chosen Ts so large ? Since,here df = is not so small !! And in the output there shall be values at 4 different frequencies i.e. 80, 30, 10, 1 Hz, Then why is not there only 4 peaks ? Kindly explain this. Sep 26 at 9:52
• As $df$ no longer has any influence, $T$ could also be chosen smaller. Then the peaks at the corners will show more of the Runge phenomenon, the peaks will become wider. As to there being more than 4 frequencies, kindly observe the division by $t$ and compare with the (inverse) Fourier transform of the box function. Sep 26 at 9:58
• It is an artefact of partial sums of Fourier series. As the DFT approximates these somewhat, you will get it also in this more discrete situation. The usual way to avoid or reduce it is to use windowing with sufficiently smooth window functions. Sep 26 at 10:31
• Ok. I would like to request yout to suggest me some relevant illustrations or references on this as I would require this in my work and thus I am keener to learn these things in details and earnestly hope for your kind aid and coopertion. Thanks. And just to ascertain t = np.linspace(0, Ts, Ns+1)[:-1] + 1e-10 - Ts/2 shall be a better choice than t = np.linspace(0, Ts, Ns+1)[:-1] used earlier right ? Sep 26 at 10:41
• Sorry, used the wrong name, Runge is overshoot in polynomial interpolation, here it is en.wikipedia.org/wiki/Gibbs_phenomenon. Sep 26 at 10:56