2

I am using fftn in numpy to generate a discrete FT of a 20 cell 1D array and of a 20x20x20 3D array and want to suppress the high frequency terms, starting with the highest frequency and extending towards the lower frequency ones. I am more familiar with continuous FTs and struggling to identify the high frequency terms in the DFT. Where should I be looking in the array that fftn generates? (I am planning to set the values of these to zero before back-transforming.)

James
  • 387
  • 1
  • 3
  • 8
  • See: http://dsp.stackexchange.com/questions/6220/why-is-it-a-bad-idea-to-filter-by-zeroing-out-fft-bins So you might want to do a zero-padded convolution filter operation instead. – hotpaw2 Feb 24 '16 at 00:58

1 Answers1

3

According to the Notes in numpy.fft.fftn documentation

The output, analogously to fft, contains the term for zero frequency in the low-order corner of all axes, the positive frequency terms in the first half of all axes, the term for the Nyquist frequency in the middle of all axes and the negative frequency terms in the second half of all axes, in order of decreasingly negative frequency.

Note however that for odd sized arrays the Nyquist frequency is not represented. Also, assuming you are dealing with real-valued signals the Discrete Fourier Transform will have Hermitian symmetry. Whenever you process those signal in the frequency domain it is important to maintain that symmetry if you want the signal to remain real-valued after the inverse transform. While zeroing out frequency component, this means you should also zero out the frequency component at the corresponding negative frequency.

What that means specifically for your 20 cell 1D array (say array x), the L highest frequency bins (including L/2 positive frequencies and L/2 negative frequencies) are thus

lower = (len(x)-L)/2+1
upper = (len(x)+L)/2+1
x[lower:upper]

Similarly, for your 20x20x20 3D array (say array y), the L highest frequency bins along each axis are:

lower = [(d-L)/2+1 for d in y.shape]
upper = [(d+L)/2+1 for d in y.shape]
y[lower[0]:upper[0],:,:] # middle of axis 0
y[:,lower[1]:upper[1],:] # middle of axis 1
y[:,:,lower[2]:upper[2]] # middle of axis 2

Now assuming that the ringing effects described in this post by hotpaw2 is not an issue for your application, you may then zero-out those bins with:

import numpy as np;
L     = 3 # number of bins to zero out along each axis (change this to fit your needs)
          # should be odd for even length array, and even for odd length array

# Following assumes x is the 1D array
lower = (len(x)-L)/2+1
upper = (len(x)+L)/2+1
x[lower:upper] = 0 # zero-out in the middle

# Following assume y is the 3D array
lower = [(d-L)/2+1 for d in y.shape]
upper = [(d+L)/2+1 for d in y.shape]
y[lower[0]:upper[0],:,:] = 0 # zero-out in the middle of axis 0
y[:,lower[1]:upper[1],:] = 0 # zero-out in the middle of axis 1
y[:,:,lower[2]:upper[2]] = 0 # zero-out in the middle of axis 2
Community
  • 1
  • 1
SleuthEye
  • 14,379
  • 2
  • 32
  • 61
  • great - I'll try these. I saw the documentation but some things like the "term for the Nyqvist frequency" did not make sense to me - I was just expecting a set of frequencies and amplitudes. The significance of the effects of the truncation on the back-transformed function remains to be seen. Good point about maintaining symmetry – James Feb 25 '16 at 15:41
  • just tried this with the 1D array and looks to be working. All the imaginary terms are essentially zero in the back-transformed function so the symmetrizing is doing the job – James Feb 25 '16 at 15:58
  • Thanks for your answer. A few remarks : all float divisions should be changed to integer ones, e.g. `lower = (len(x)-L)//2+1` and mention that what you did is a low-pass filter in frequency. – FraSchelle Feb 05 '20 at 11:54