Added FFT function (uses np.fft.fft) to utilities module, to avoid repeating same code in several modules.

gprMax/grid.py

@@ -30,6 +30,7 @@ from gprMax.constants import floattype
 from gprMax.constants import complextype
 from gprMax.materials import Material
 from gprMax.pml import PML
+from gprMax.utilities import fft_power
 from gprMax.utilities import round_value
 
 
@@ -234,45 +235,32 @@ def dispersion_analysis(G):
             else:
                 # User-defined waveform
                 if waveform.type == 'user':
-                    waveformvalues = waveform.uservalues
+                    iterations = G.iterations
 
                 # Built-in waveform
                 else:
-                    # Time to analyse waveform - 4*pulse_width as using entire
+                # Time to analyse waveform - 4*pulse_width as using entire
-                    # time window can result in demanding FFT
+                # time window can result in demanding FFT
                     waveform.calculate_coefficients()
-                    time = np.arange(0, 4 * waveform.chi, G.dt)
+                    iterations = round_value(4 * waveform.chi / G.dt)
-                    waveformvalues = np.zeros(len(time))
+                    if iterations > G.iterations:
-                    timeiter = np.nditer(time, flags=['c_index'])
+                        iterations = G.iterations
 
-                    while not timeiter.finished:
+                waveformvalues = np.zeros(G.iterations)
-                        waveformvalues[timeiter.index] = waveform.calculate_value(timeiter[0], G.dt)
+                for iteration in range(G.iterations):
-                        timeiter.iternext()
+                    waveformvalues[iteration] = waveform.calculate_value(iteration * G.dt, G.dt)
 
                 # Ensure source waveform is not being overly truncated before attempting any FFT
                 if np.abs(waveformvalues[-1]) < np.abs(np.amax(waveformvalues)) / 100:
-                    # Calculate magnitude of frequency spectra of waveform
+                    # FFT
-                    mag = np.abs(np.fft.fft(waveformvalues))**2
+                    freqs, power = fft_power(waveformvalues, G.dt)
-
-                    # Calculate power (ignore warning from taking a log of any zero values)
-                    with np.errstate(divide='ignore'):
-                        power = 10 * np.log10(mag)
-                    # Replace any NaNs or Infs from zero division
-                    power[np.invert(np.isfinite(power))] = 0
-
-                    # Frequency bins
-                    freqs = np.fft.fftfreq(power.size, d=G.dt)
-
-                    # Shift powers so that frequency with maximum power is at zero decibels
-                    power -= np.amax(power)
-
                     # Get frequency for max power
-                    freqmaxpower = np.where(np.isclose(power[1::], np.amax(power[1::])))[0][0]
+                    freqmaxpower = np.where(power == 0)[0][0]
 
                     # Set maximum frequency to a threshold drop from maximum power, ignoring DC value
                     try:
-                        freq = np.where((np.amax(power[freqmaxpower::]) - power[freqmaxpower::]) > G.highestfreqthres)[0][0] + 1
+                        freqthres = np.where(power[freqmaxpower:] < -G.highestfreqthres)[0][0]
-                        results['maxfreq'].append(freqs[freq])
+                        results['maxfreq'].append(freqs[freqthres])
                     except:
                         results['error'] = 'unable to calculate maximum power from waveform, most likely due to undersampling.'
 

gprMax/utilities.py

@@ -138,6 +138,35 @@ def round32(value):
     return int(32 * np.ceil(float(value) / 32))
 
 
+def fft_power(waveform, dt):
+    """Calculate a FFT of the given waveform of amplitude values;
+        converted to decibels and shifted so that maximum power is 0dB
+
+    Args:
+        waveform (ndarray): time domain waveform
+        dt (float): time step
+
+    Returns:
+        freqs (ndarray): frequency bins
+        power (ndarray): power
+    """
+
+    # Calculate magnitude of frequency spectra of waveform (ignore warning from taking a log of any zero values)
+    with np.errstate(divide='ignore'): #
+        power = 10 * np.log10(np.abs(np.fft.fft(waveform))**2)
+
+    # Replace any NaNs or Infs from zero division
+    power[np.invert(np.isfinite(power))] = 0
+
+    # Frequency bins
+    freqs = np.fft.fftfreq(power.size, d=dt)
+
+    # Shift powers so that frequency with maximum power is at zero decibels
+    power -= np.amax(power)
+
+    return freqs, power
+
+
 def human_size(size, a_kilobyte_is_1024_bytes=False):
     """Convert a file size to human-readable form.
 

tools/plot_Ascan.py

@@ -27,6 +27,7 @@ import matplotlib.gridspec as gridspec
 
 from gprMax.exceptions import CmdInputError
 from gprMax.receivers import Rx
+from gprMax.utilities import fft_power
 
 
 def mpl_plot(filename, outputs=Rx.defaultoutputs, fft=False):
@@ -83,22 +84,17 @@ def mpl_plot(filename, outputs=Rx.defaultoutputs, fft=False):
 
             # Plotting if FFT required
             if fft:
-                # Calculate magnitude of frequency spectra of waveform (ignore warning from taking a log of any zero values)
+                # FFT
-                with np.errstate(divide='ignore'):
+                freqs, power = fft_power(outputdata, dt)
-                    power = 10 * np.log10(np.abs(np.fft.fft(outputdata))**2)
+                freqmaxpower = np.where(power == 0)[0][0]
-                # Replace any NaNs or Infs from zero division
-                power[np.invert(np.isfinite(power))] = 0
 
-                # Frequency bins
+                # Set plotting range to -60dB from maximum power or 4 times
-                freqs = np.fft.fftfreq(power.size, d=dt)
+                # frequency at maximum power
+                try:
+                    pltrange = np.where(power[freqmaxpower:] < -60)[0][0] + freqmaxpower + 1
+                except:
+                    pltrange = freqmaxpower * 4
 
-                # Shift powers so that frequency with maximum power is at zero decibels
-                power -= np.amax(power)
-
-                # Set plotting range to -60dB from maximum power
-                pltrange = np.where((np.amax(power[1::]) - power[1::]) > 60)[0][0] + 1
-                # To a maximum frequency
-                # pltrange = np.where(freqs > 2e9)[0][0]
                 pltrange = np.s_[0:pltrange]
 
                 # Plot time history of output component

tools/plot_source_wave.py

@@ -24,6 +24,7 @@ import numpy as np
 import matplotlib.pyplot as plt
 
 from gprMax.exceptions import CmdInputError
+from gprMax.utilities import fft_power
 from gprMax.utilities import round_value
 from gprMax.waveforms import Waveform
 
@@ -83,12 +84,7 @@ def mpl_plot(w, timewindow, dt, iterations, fft=False):
 
     print('Waveform characteristics...')
     print('Type: {}'.format(w.type))
-
+    print('Maximum (absolute) amplitude: {:g}'.format(np.max(np.abs(waveform))))
-    if w.type == 'user':
-        waveform = w.uservalues
-        w.amp = np.max(np.abs(waveform))
-
-    print('Maximum amplitude: {:g}'.format(w.amp))
 
     if w.freq and not w.type == 'gaussian':
         print('Centre frequency: {:g} Hz'.format(w.freq))
@@ -100,44 +96,20 @@ def mpl_plot(w, timewindow, dt, iterations, fft=False):
         delay = np.sqrt(2) / w.freq
         print('Time to centre of pulse: {:g} s'.format(delay))
 
-    # Calculate pulse width for gaussian
-    if w.type == 'gaussian':
-        powerdrop = -3  # dB
-        with np.errstate(divide='ignore'): # Ignore warning from taking a log of any zero values
-            startpower = 10 * np.log10(waveform / np.amax(waveform))
-            stopower = 10 * np.log10(waveform[start:] / np.amax(waveform))
-
-        # Replace any NaNs or Infs from zero division
-        startpower[np.invert(np.isfinite(startpower))] = 0
-        stoppower[np.invert(np.isfinite(stoppower))] = 0
-
-        start = np.where(startpower > powerdrop)[0][0]
-        stop = np.where(stoppower < powerdrop)[0][0] + start
-        print('Pulse width at {:d}dB, i.e. full width at half maximum (FWHM): {:g} s'.format(powerdrop, time[stop] - time[start]))
-
     print('Time window: {:g} s ({} iterations)'.format(timewindow, iterations))
     print('Time step: {:g} s'.format(dt))
 
     if fft:
-        # Calculate magnitude of frequency spectra of waveform (ignore warning from taking a log of any zero values)
+        # FFT
-        with np.errstate(divide='ignore'): #
+        freqs, power = fft_power(waveform, dt)
-            power = 10 * np.log10(np.abs(np.fft.fft(waveform))**2)
 
-        # Replace any NaNs or Infs from zero division
+        # Set plotting range to 4 times frequency at max power of waveform or
-        power[np.invert(np.isfinite(power))] = 0
+        # 4 times the centre frequency
-
+        freqmaxpower = np.where(power == 0)[0][0]
-        # Frequency bins
+        if freqs[freqmaxpower] > w.freq:
-        freqs = np.fft.fftfreq(power.size, d=dt)
+            pltrange = np.where(freqs > 4 * freqs[freqmaxpower])[0][0]
-
+        else:
-        # Shift powers so that frequency with maximum power is at zero decibels
+            pltrange = np.where(freqs > 4 * w.freq)[0][0]
-        power -= np.amax(power)
-
-        if w.type == 'user':
-            freqmaxpower = np.where(np.isclose(power[1::], np.amax(power[1::])))[0][0]
-            w.freq = freqs[freqmaxpower]
-
-        # Set plotting range to 4 times centre frequency of waveform
-        pltrange = np.where(freqs > 4 * w.freq)[0][0]
         pltrange = np.s_[0:pltrange]
 
         fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, num=w.type, figsize=(20, 10), facecolor='w', edgecolor='w')