Blast and Impact Dynamics

Department of Civil & Structural Engineering

dispersion.m - A MatLab script for phase angle and amplitude correction of pressure bar signals


In processing the signals from split Hopkinson pressure bar (SHPB) experiments it is often assumed that longitudinal stress waves in the pressure bars propagate one-dimensionally at a common velocity c_0, and so measurements taken at the strain gauges are often simply translated to the end of the bar using a suitable time delay. In reality, stress waves propagate at a specific phase velocity, c_p, which is a function of frequency and the bar's diameter, one-dimensional wave speed and Poisson's ratio. Phase velocity decreases as the frequency of a wave increases, leading to dispersion of a signal as it propagates down the bar. Dispersion of the stress pulse is accompanied by a frequency-dependent variation in stress and strain across the bar cross-section, so that a signal recorded on the surface of the bar at some distance from the specimen will not accurately describe the stresses the specimen was subjected to, and hence cannot be used to accurately determine the specimen response.

This script uses an implementation of the dispersion-correction method described by Tyas and Pope (2005) to ensure that the inferred measurements of axial stress and strain accurately represent the specimen behaviour. In this method:

  1. The time-domain strain signal is converted into the frequency domain using the fast Fourier transform (FFT)
  2. A correction is applied to the phase angle of each frequency component to account for the dispersion over the distance between the strain gauge and the bar end, using Bancroft's equation.
  3. A correction is applied to the amplitude of each frequency component using the factors M_1 and M_2, which account for the variation of strain and Young's modulus across the bar cross section, respectively. These are derived from Davies' analysis of the radial effects in a cylindrical pressure bar.
  4. The signal is transformed back into the time domain using the inverse FFT.

Dispersion.m uses a pre-calculated, normalised look-up table of phase velocity, M1 and M2 to improve calculation time. A lookup table for a Poisson's ratio of 0.29 has been provided in the .zip file, and other tables can be constructed using the relationships defined in Tyas and Pope (2005).

People involved

Research project(s)

Strain-rate effects in quartz sand

Soil-filled wire and geotextile gabions are commonly used to construct perimeter walls and other defensive infrastructure in military bases. They serve to protect personnel and key assets from the effects of blast and fragmentation. The attenuating properties of soil make it a highly effective defence against such threats, and as...


% dispersion() - First-mode dispersion correction of a finite  
% arbitrary signal in a cylindrical bar
% - Finds FFT of the signal
% - Corrects phase velocity and amplitude of each frequency
%   using method described by Tyas & Pope (2005)
% - Reconstructs signal using IFFT
% - (Frequencies above fa/c0 = 0.2619 stripped (d/L = 0.6) 
%   due to limitations of M1 correction)
% x        Zero-padded strain signal in time domain 
%          (1xN numeric)
% fs       Sampling frequency, Hz
% a        Bar radius, m
% c0       One-dimensional wave velocity of the bar, m/s
% E        Young's modulus of the bar, GPa
% z        Distance to correct over, m 
%          (+ve in direction of propagation)
% xStrain  Dispersion-corrected strain signal
% xStress  Dispersion-corrected stress signal, MPa
function [xStrain xStress] = dispersion(x,fs,a,c0,E,z)
    % Input signal
    N = length(x); % Number of elements in signal
    dt = 1/fs; % Time step, s 
    t = 0:dt:dt*(N-1); % Time base, s
    f = (0:N-1)*(fs/N); % FFT frequencies, Hz
    fMax = 0.2619*c0/a; % Max correctable frequency due to factor M1 limitations, Hz
    % FFT the signal
    X = fft(x);
    XStrain = X; % Create a copy for strain correction
    XStress = X; % Create a copy for stress correction
    % Phase shift, adjust magnitude of frequency components
    numberOfBins = length(X);
    DCbin = 1; % DC frequency bin
        % N is even
        positiveBins = 2:numberOfBins/2; % Positive frequency bins
        nyquistBin = numberOfBins/2+1; % Nyquist frequency bin
        binsToEdit = [positiveBins nyquistBin]; % Total bins to edit individually
        negativeBins = numberOfBins/2+2:numberOfBins; % Negative frequency bins (populated through conjugation of positive bins)
        % N is odd
        positiveBins = 2:(numberOfBins+1)/2; % Positive frequency bins
        binsToEdit = positiveBins; % Total bins to edit individually
        negativeBins = (numberOfBins+1)/2+1:numberOfBins; % Negative frequency bins (populated through conjugation of positive bins)
    for b = binsToEdit
        if f(b) <= fMax
            % Edit phase and amplitude of positive bins
            [angleMod M1 M2] = dispersionFactors(f(b),a,c0,z); % Find phase shift and factors M1 and M2 for current freq
            XStrain(b) = M1*abs(X(b)) * exp(1i * (angle(X(b))-angleMod)); % Apply phase shift and factors M1 to obtain corrected strain
            XStress(b) = M1*M2*abs(X(b)) * exp(1i * (angle(X(b))-angleMod)); % Apply phase shift and factors M1 and M2 to obtain corrected stress(/E)
            % Above fMax zero X data (apply perfect low-pass filter)
            XStrain(b) = 0;
            XStress(b) = 0;
    XStrain(negativeBins) = conj(XStrain(positiveBins(end:-1:1))); % Correct negative bins by taking complex conjugate of positive bins
    XStress(negativeBins) = conj(XStress(positiveBins(end:-1:1))); % Correct negative bins by taking complex conjugate of positive bins
    % Convert the corrected frequency components back into the time domain
    xStrain = ifft(XStrain); % Corrected strain
    xStress = ifft(XStress)*E*1000; % Corrected stress, MPa
% dispersionFactors() - Calculate corrections to amplitude 
% and phase angle to account for dispersion at a particular 
% frequency
% <dataTable>  A .mat file containing pre-calculated vectors 
%              for a particular Poisson's ratio
%              - normFreqs: normalised frequencies (f*a/c0)
%              - vRatios: normalised velocities (c/c0)
%              - M1: amplitude factor M1
%              - M2: normalised amplitude factor M2 (M2/E)
% f            Frequency, Hz
% a            Bar radius, m
% c0           One-dimensional wave velocity of the bar, m/s
% z            Distance to apply correction over, m
% angleMod  Phase angle correction, rad
% M1        Correction for variation in response across bar 
%           cross-section
% M2        Correction for variation in ratio of axial stress
%           and axial strain (dynamic Young's modulus)
function [angleMod M1 M2] = dispersionFactors(f,a,c0,z)
    load('dispersionFactors_PR29.mat'); % File containing phase velocity data
    normFreq = f*a/c0; % Normalised frequency
	% Find change in phase angle
	phaseVelocity = interp1(normFreqs,vRatios,normFreq)*c0; % Interpolated phase velocity value
	angleMod = 2 * pi * f * z / phaseVelocity; % Change in phase angle at normFreq
    % Find amplitude factors M1 and M2
	M1 = interp1(normFreqs,M1,normFreq); % Interpolated value of M1 at normFreq
	M2 = interp1(normFreqs,M2,normFreq); % Interpolated value of M2/E at normFreq