Am I supposed to normalize FFT in Python?

Question

I'm sorry if this question is posed often, I can't seem to understand the answers already out there. I'm working with several forms of the Fourier transform, including the FFT, PSD, and spectrograms. I'm not sure whether, upon computing the Fourier transform of my signal, I'm supposed to normalize the result by some factor. For context I don't have any interest in comparing signals, I just need my outputs to be mathematically accurate.

I've seen some answers that say the DFT has to be normalized by a 1/N factor - I'm interpreting this as the 1/N factor in the IFFT, is that correct? Or is there some other reason you would have to normalize by 1/N?

What's really confusing me though, is answers where I'm told to divide by the sampling frequency. This answer implies that it is equivalent to dividing by N- does this have something to do with the transform being computed over 1-second segments or something along those lines?

A tutorial I'm working from normalizes this way, and I'm not sure where it comes from. They construct a whitening procedure like so:

def whiten(strain, interp_psd, dt):     
    Nt = len(strain)
    freqs = np.fft.rfftfreq(Nt, dt)

    # whitening: transform to freq domain, divide by asd, then transform back, 
    # taking care to get normalization right.
    hf = np.fft.rfft(strain)
    norm = 1./np.sqrt(1./(dt*2))
    white_hf = hf / np.sqrt(interp_psd(freqs)) * norm
    white_ht = np.fft.irfft(white_hf, n=Nt)
    return white_ht

where interp_psd is a power spectral density interpolant of the noise, and is computed in the following way:

NFFT = 4*fs
Pxx_H1, freqs = mlab.psd(strain_H1, Fs = fs, NFFT = NFFT)
Pxx_L1, freqs = mlab.psd(strain_L1, Fs = fs, NFFT = NFFT)

# We will use interpolations of the ASDs computed above for whitening:
psd_H1 = interp1d(freqs, Pxx_H1)
psd_L1 = interp1d(freqs, Pxx_L1)

Why do they multiply the whitened strain in the frequency domain by dt*2? I understand that 2 could come from discarding negative frequencies, but wouldn't the RFFT take that into account?

The tutorial later on computes the spectrogram of the filtered data:

# Plot the H1 spectrogram:
plt.figure(figsize=(10,6))
spec_H1, freqs, bins, im = plt.specgram(strain_H1[indxt], NFFT=NFFT, Fs=fs, window=window, 
                                            noverlap=NOVL, cmap=spec_cmap, xextent=[-deltat,deltat])
plt.xlabel('time (s) since '+str(tevent))
plt.ylabel('Frequency (Hz)')
plt.colorbar()
plt.axis([-deltat, deltat, 0, 2000])
plt.title('aLIGO H1 strain data near '+eventname)
plt.savefig(eventname+'_H1_spectrogram.'+plottype)

Here, they don't normalize, but I recognize that they might not need to since they're only interested in the visual of the spectrogram.

Later they normalize by the sampling frequency when performing a matched-filter exercise, but then reverse it.

# Take the Fourier Transform (FFT) of the data and the template (with dwindow)
    data_fft = np.fft.fft(data*dwindow) / fs

    # -- Interpolate to get the PSD values at the needed frequencies
    power_vec = np.interp(np.abs(datafreq), freqs, data_psd)

    # -- Calculate the matched filter output in the time domain:
    # Multiply the Fourier Space template and data, and divide by the noise power in each frequency bin.
    # Taking the Inverse Fourier Transform (IFFT) of the filter output puts it back in the time domain,
    # so the result will be plotted as a function of time off-set between the template and the data:
    optimal = data_fft * template_fft.conjugate() / power_vec
    optimal_time = 2*np.fft.ifft(optimal)*fs

I apologize if this is too much information. All in all, my questions are these: do I have to normalize the output of a FFT in python (numpy, scipy, matplotlib) in order to be mathematically accurate, and by what factor? And does that normalization differ for transforms like the PSD or the spectrogram?

Answer 1

Your decision to normalize or not does not change the accuracy of your answer, as it is simply a scaling factor. If you use the common scaling of $1/N$, then the output for each DFT bin will represent the average of the portion of the input signal that is at the frequency defined by that bin, scaled to the same units as the input. So that is convenient and gives a certain meaning to what the magnitude of the output represents, but if you didn't divide by N it doesn't make the answer incorrect, as long as your approach is consistent.

For example, it is also common to use $1/\sqrt{N}$ as this would make the DFT and IDFT completely symmetric as the same scaling can be used for both.

To understand what is happening, consider the simplest case of the computation of the first bin of the DFT, which represents DC, and consider a DC signal that is 1 for all $N$ samples given. If this DC signal represents a magnitude in voltage, then this is samples of $f(t) = 1$ volt, given as $f[n]=1$, for all $N$ samples, $n = 0$ to $N-1$.

The first bin of the DFT without any normalization is simply the sum over all N samples:

$$F[k=0] = \sum_{n=0}^{N-1}v[n] = N$$

So if we were to scale this by $1/N$, then the result would represent the average for the frequency component of our signal within that bin in the same units as the input. In this case all values were at bin 0 "DC", so the average value was 1V, consistent with what we get when we scale by $1/N$. (It is clear from the formula above that when scaled by $1/N$ this is simply the average of $f[n]$).

Similarly all the other bins are the result of the average after multiplying by the time domain signal by $e^{-jk\omega_o n}$, which can be viewed as frequency translating any signal at that bin to DC and then taking the average of that result.

Answer 2

Whether you want to normalize or not depends on whether you want to know the level or the energy of the DFT input.

IIRC, the SciPy FFT returns energy (complies with Parseval’s relation). A signal N times as long at the same level has N times more energy. So you could divide by N to get an estimation of a level instead of energy.