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?