diff --git a/gprMax/grid.py b/gprMax/grid.py index 4d10f4cf..f5a34170 100644 --- a/gprMax/grid.py +++ b/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 window can result in demanding FFT + # Time to analyse waveform - 4*pulse_width as using entire + # time window can result in demanding FFT waveform.calculate_coefficients() - time = np.arange(0, 4 * waveform.chi, G.dt) - waveformvalues = np.zeros(len(time)) - timeiter = np.nditer(time, flags=['c_index']) + iterations = round_value(4 * waveform.chi / G.dt) + if iterations > G.iterations: + iterations = G.iterations - while not timeiter.finished: - waveformvalues[timeiter.index] = waveform.calculate_value(timeiter[0], G.dt) - timeiter.iternext() + waveformvalues = np.zeros(G.iterations) + for iteration in range(G.iterations): + 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 - mag = np.abs(np.fft.fft(waveformvalues))**2 - - # 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) - + # FFT + freqs, power = fft_power(waveformvalues, G.dt) # 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 - results['maxfreq'].append(freqs[freq]) + freqthres = np.where(power[freqmaxpower:] < -G.highestfreqthres)[0][0] + results['maxfreq'].append(freqs[freqthres]) except: results['error'] = 'unable to calculate maximum power from waveform, most likely due to undersampling.' diff --git a/gprMax/utilities.py b/gprMax/utilities.py index a5f9302e..f6496cc1 100644 --- a/gprMax/utilities.py +++ b/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. diff --git a/tools/plot_Ascan.py b/tools/plot_Ascan.py index 66bb4675..921332c4 100644 --- a/tools/plot_Ascan.py +++ b/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) - with np.errstate(divide='ignore'): - power = 10 * np.log10(np.abs(np.fft.fft(outputdata))**2) - # Replace any NaNs or Infs from zero division - power[np.invert(np.isfinite(power))] = 0 + # FFT + freqs, power = fft_power(outputdata, dt) + freqmaxpower = np.where(power == 0)[0][0] - # Frequency bins - freqs = np.fft.fftfreq(power.size, d=dt) + # Set plotting range to -60dB from maximum power or 4 times + # 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 diff --git a/tools/plot_source_wave.py b/tools/plot_source_wave.py index 40fe10bc..5b20ecfc 100644 --- a/tools/plot_source_wave.py +++ b/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)) - - if w.type == 'user': - waveform = w.uservalues - w.amp = np.max(np.abs(waveform)) - - print('Maximum amplitude: {:g}'.format(w.amp)) + print('Maximum (absolute) amplitude: {:g}'.format(np.max(np.abs(waveform)))) 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) - with np.errstate(divide='ignore'): # - power = 10 * np.log10(np.abs(np.fft.fft(waveform))**2) + # FFT + freqs, power = fft_power(waveform, dt) - # 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) - - 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] + # Set plotting range to 4 times frequency at max power of waveform or + # 4 times the centre frequency + freqmaxpower = np.where(power == 0)[0][0] + if freqs[freqmaxpower] > w.freq: + pltrange = np.where(freqs > 4 * freqs[freqmaxpower])[0][0] + else: + 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')