Analyze a shock signal

The line-scanning technique has been developed and widely used to determine the position of shocks, particularly for normal shocks or those close to normal. In this example, the core method of the shock tracking algorithm is based on detecting the variation in the maximum density gradient area. The methodology is detailed in this artical. From the generated slice list (discussed in detail in the Slice list generation example), the shock will be tracked, and the corresponding signal will be generated for analysis.

Steps are as following:

Import the slice list and define important parameters:

import cv2
import numpy as np
from ShockOscillationAnalysis import SOA

if __name__ == '__main__':
   # define the slice list file
    imgdirc = r'results\Slicelist_test-results'
   imgpath = fr'{imgdirc}\2.0kHz_10mm_0.12965964343598055mm-px_tk_60px_-SliceList.png'

   f = 2000    # images sampling rate

   # from the file name or can be passed directly from:
   # SliceListGenerator.GenerateSlicesArray function
   scale = 0.12965964343598055  # mm/px

   # import the slice list
   slicelist = cv2.imread(imgpath)
   n = slicelist.shape[0]  # time

   # iniate the ShockOscillationAnalysis module
   SA = SOA(f)

   # spacify the shock region (Draw 2 vertical lines)
   newref = SA.LineDraw(slicelist, 'V', 0, Intialize=True)
   newref = SA.LineDraw(SA.clone, 'V', 1)
   # to make sure the spacified lines are correctly sorted
   newref.sort()
   # to crop the slicelist to the shock region
   shock_region = slicelist[:, newref[0]:newref[1]]
   # the width of the slice list in pixels
   xPixls = (newref[1]-newref[0])
   # the width of the slice list in mm
   shock_region_mm = xPixls*scale
   # Add info to log.txt file
   newlog = f'Shock Regions: {newref},\t Represents: {xPixls}px, \t Shock Regions in mm:{shock_region_mm}'
   SA.log(newlog, imgdirc)
   print(newlog)

Note

SA.clone is the modified image to see the first line, to keep original without modification.

$_static\img\T4\T4-1.png$

Important

The tracking domain must not contain more than one shock, as this will confuse the software and generate incorrect results.

And the console output of this step is:

Registered line: 255
Registered line: 441
Shock Regions: [255, 441],    Represents: 186px,      Shock Regions in mm:24.116693679092382

To improve the traking quality, it is good to clean optical defects by subtracting Average slice from all slices:

#%% slice list cleaning
# [subtracting the average, subtracting ambiant light frequency, improve brightness/contrast/sharpness]
ShockwaveRegion = SA.CleanSnapshots(ShockwaveRegion,'Average')

The console output of this step is:

Improving image quality ...
      - subtracting Averaging ... ✓

Note

Clean illumination disturbances by Fast Fourier Transform (FFT) also can be done by adding FFT and other parameters as follow.

ShockwaveRegion = SA.CleanSnapshots(ShockwaveRegion,
                                    'Average','FFT',
                                    filterCenter = [(0, 25)], D = 20, n = 5,
                                    ShowIm = True)

$_static\img\T4\FFT.png$ $_static\img\T4\FFT-Filtered.png$

filterCenter and other FFT parameters can be determined by enebling ShowIm to detect the defect frequency location.
The filterCenter may contain multiple locations.
The cleaning process follows the order of the argument, in the above example the Averaging will take place first then FFT.
Additional parameters such as Brightness/Contrast may also appended to the arguments if required check SOA.CleanSnapshots.

To track the shock and generate the shock signal and scale it.

import matplotlib.pyplot as plt

#%% Find shock location
shock_loc_px, uncer = SA.ShockTrakingAutomation(shock_region,
                                                method = 'integral',        # there is also 'maxGrad' and 'darkest_spot'
                                                reviewInterval = [11,14],   # to review the tracking process within this range
                                                Signalfilter = 'med-Wiener')

# Add info to log.txt file
 newlog = f'uncertainty ratio: {(len(uncer)/len(shock_loc_px))*100:0.2f}%'
 SA.log(newlog, imgdirc)
 print(newlog)

# unpack and scale the output values
shock_loc_mm= scale * np.array(shock_loc_px)     # to scale the shock location output to mm

snapshot_indx, uncertain, reason = zip(*uncer)   # unpack uncertainity columns
uncertain_mm = scale * np.array(uncertain)       # to scale the uncertain locatshock location output to mm

# plotting the output
fig1, ax1 = plt.subplots(figsize=(8,50))
# shock region image as background to review the tracked points
ax1.imshow(shock_region, extent=[0, shock_region_mm, shock_region.shape[0], 0], aspect='0.1', cmap='gray')
ax1.plot(shock_loc_mm, range(n),'x', lw = 1, color = 'g', ms = 7)        # To plot the detected shock locations
ax1.plot(uncertain_mm, snapshot_indx,'x', lw = 1, color = 'r', ms = 5)   # To plot the uncertain shock points

The tracking review:

$_static\img\T4\T4-5c.png$ $_static\img\T4\T4-2c.png$ $_static\img\T4\T4-3c.png$ $_static\img\T4\T4-4c.png$

The out put results:

$_static\img\T4\T4-6.png$

The console output of this step is:

Processing the shock location using integral method...
[====================] 100%
Appling med-Wiener filter...
Processing time: 0 Sec
uncertainty ratio: 12.00%

Note

Mostly, the tracked points follow the shock location; however, the uncertainty ratio is quite high at 14%.
The reasons for uncertainty can be reviewed from the uncertainty output. Based on this review, users may choose to change the strategy by adjusting the cleaning parameters or their order. Additionally, the selected range of the shock could be a parameter to consider.

Finally, shift the signal by the average value and use welch method to study the power spectral density (PSD).

from scipy import signal

#%% Apply welch method for PSD
avg_shock_loc = np.average(shock_loc_mm)      # find the average shock location
shock_loc_mm = shock_loc_mm - avg_shock_loc   # to shift the signal to the average

# Calculate the PSD
Freq, psd = signal.welch(x = shock_loc_mm, fs = f, window='barthann',
                           nperseg = 512, noverlap=0, nfft=None, detrend='constant',
                           return_onesided=True, scaling='density')

fig,ax = plt.subplots(figsize=(10,10))
ax.loglog(Freq, psd, lw = '2')
ax.set_ylabel(r"PSD [mm$^2$/Hz]")
ax.set_xlabel("Frequency [Hz]")

The out put results:

$_static\img\T4\T4-7.png$

The full code example:

import cv2
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
from ShockOscillationAnalysis import SOA

if __name__ == '__main__':
   # define the slice list file
   imgdirc = r'results\Slicelist_test-results'
   imgpath = fr'{imgdirc}\2.0kHz_10mm_0.12965964343598055mm-px_tk_60px_-SliceList.png'

   f = 2000    # images sampling rate

   # from the file name or can be passed directly from:
   # SliceListGenerator.GenerateSlicesArray function
   scale = 0.12965964343598055  # mm/px

   # import the slice list
   slicelist = cv2.imread(imgpath)
   n = slicelist.shape[0]  # time

   # iniate the ShockOscillationAnalysis module
   SA = SOA(f)

   # spacify the shock region (Draw 2 vertical lines)
   newref = SA.LineDraw(slicelist, 'V', 0, Intialize=True)
   newref = SA.LineDraw(SA.clone, 'V', 1)
   # to make sure the spacified lines are correctly sorted
   newref.sort()
   # to crop the slicelist to the shock region
   shock_region = slicelist[:, newref[0]:newref[1]]
   # the width of the slice list in pixels
   xPixls = (newref[1]-newref[0])
   # the width of the slice list in mm
   shock_region_mm = xPixls*scale
   # Add info to log.txt file
   newlog = f'Shock Regions: {newref},\t Represents: {xPixls}px, \t Shock Regions in mm:{shock_region_mm}'
   SA.log(newlog, imgdirc)
   print(newlog)

   # %% slice list cleaning
   # [subtracting the average, subtracting ambiant light frequency,
   # improve brightness/contrast/sharpness]
   shock_region = SA.CleanSnapshots(shock_region, 'Average')

   # %% Find shock location
   shock_loc_px, uncer = SA.ShockTrakingAutomation(shock_region,
                                                   method='integral',        # There is also 'maxGrad' and 'darkest_spot'
                                                   reviewInterval=[11, 14],  # to review the tracking process within this range
                                                   Signalfilter='med-Wiener')

   # Add info to log.txt file
   newlog = f'uncertainty ratio: {(len(uncer)/len(shock_loc_px))*100:0.2f}%'
   SA.log(newlog, imgdirc)
   print(newlog)

   # unpack and scale the output values
   # to scale the shock location output to mm
   shock_loc_mm = scale * np.array(shock_loc_px)

   # unpack uncertainity columns
   snapshot_indx, uncertain, reason = zip(*uncer)
   # to scale the uncertain locatshock location output to mm
   uncertain_mm = scale * np.array(uncertain)

   # plotting the output
   fig1, ax1 = plt.subplots(figsize=(8, 50))
   ax1.imshow(shock_region,
               extent=[0, shock_region_mm, shock_region.shape[0], 0],
               aspect='0.1', cmap='gray')
   ax1.plot(shock_loc_mm, range(n), 'x', lw=1, color='g', ms=7)
   ax1.plot(uncertain_mm, snapshot_indx, 'x', lw=1, color='r', ms=5)

   # %% Apply welch method for PSD
   # find the average shock location
   avg_shock_loc = np.average(shock_loc_mm)
   # to shift the signal to the average
   shock_loc_mm = shock_loc_mm - avg_shock_loc

   Freq, psd = signal.welch(x=shock_loc_mm, fs=f, window='barthann',
                           nperseg=512, noverlap=0, nfft=None,
                           detrend='constant', return_onesided=True,
                           scaling='density')

   fig, ax = plt.subplots(figsize=(10, 10))
   ax.loglog(Freq, psd, lw='2')
   ax.set_ylabel(r"PSD [mm$^2$/Hz]")
   ax.set_xlabel("Frequency [Hz]")