Envelope Fitting Tool for Longitudinal Stress Wave of Wood


In [1]:
# put your file name below
# example:
wood_file = "./wav_data/SO1-1"
wood_name = "Sonokeling 1"
# wood_name = "put your file target here"

Acquiring the signal


In [2]:
from pylab import *
from scipy.io import wavfile
from scipy.signal import butter, lfilter
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import freqz
from ipykernel import kernelapp as app
import scipy as sp
import scipy.optimize

In [3]:
%matplotlib inline

In [4]:
# bandpass filter
def butter_bandpass(lowcut, highcut, fs, order=5):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a
def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    y = lfilter(b, a, data)
    return y

In [5]:
samp_freq, snd = wavfile.read(wood_file + '.wav')

In [6]:
snd.dtype


Out[6]:
dtype('int16')

In [7]:
snd = snd / (2.**15)

In [8]:
print(snd.shape)
print(snd.shape[0])
print(snd.shape[0] / samp_freq)
print(snd.shape[1])


(213444, 2)
213444
4.84
2

In [9]:
try:
    s1 = snd[:,0]
except IndexError:
    s1 = snd

In [10]:
s1 # original signal


Out[10]:
array([ -7.32421875e-04,  -7.62939453e-04,  -7.01904297e-04, ...,
         6.10351562e-05,   9.15527344e-05,   1.22070312e-04])

Calculating the envelope fitting analysis

Time domain view


In [11]:
# acquiring signal in time domain
time_array = arange(0, snd.shape[0], 1)
time_array = time_array / samp_freq
time_array = time_array * 1000  #scale to milliseconds

In [12]:
# plot the signal in time domain
plt.plot(time_array, s1, color='k')
ylabel('Amplitude')
xlabel('Time (ms)')


Out[12]:
<matplotlib.text.Text at 0x8a2c630>

In [13]:
# get the peak point as the starting point of impulse response
start = s1.argmax()
print(start)
# get the stop point of the impulse response snapshot
stop = start+samp_freq//25
start_s1 = s1[start:stop]
start_time_array = time_array[start:stop]


71264

In [14]:
# show the impulse response snapshot
plt.plot(start_time_array, start_s1, color='k')
ylabel('Amplitude')
xlabel('Time (ms)')


Out[14]:
<matplotlib.text.Text at 0x8a2cc18>

Frequency domain view

Original Signal


In [15]:
n = len(s1) 
p = fft(s1) # take the fourier transform 
n_unique_pts = int(ceil((n+1)/2.0))
p = p[0:n_unique_pts]
p = abs(p)
print(p)


[  1.87500000e+00   7.11737011e-01   3.31320289e+00 ...,   4.05069728e-03
   9.49487996e-03   2.68554688e-03]

In [16]:
# plot the FFT result
p = p / float(n) # scale by the number of points so that
# the magnitude does not depend on the length
# of the signal or on its sampling frequency
p = p**2 # square it to get the power
# multiply by two (see technical document for details)
# odd nfft excludes Nyquist point
if n % 2 > 0: # we've got odd number of points fft
    p[1:len(p)] = p[1:len(p)] * 2
else:
    p[1:len(p) -1] = p[1:len(p) - 1] * 2 # we've got even number of points fft
freq_array = arange(0, n_unique_pts, 1.0) * (samp_freq / n);
freq_db = 10*log10(p)
lower_tresh = int(len(freq_db)/(samp_freq/2)*800)
higher_tresh = int(len(freq_db)/(samp_freq/2)*4100)
natural_freq_mag = max(freq_db[lower_tresh:higher_tresh])
natural_freq_index = freq_db[lower_tresh:higher_tresh].argmax()
gca().set_position((.1, 0.18, 1, 1))
plt.plot(freq_array/1000, freq_db, '-k')
xlabel('Frequency (kHz)')
ylabel('Power (dB)')
title(wood_name + ' Wood Longitudinal Wave. Feliana and Ayutyastuti')
figtext(.02, .02,
        'Estimated Natural Frequency = %0.2f kHz %0.2f dB ' % (freq_array[lower_tresh+natural_freq_index]/1000, natural_freq_mag))


Out[16]:
<matplotlib.text.Text at 0x95f0e48>

Filtered signal near estimated natural frequency


In [17]:
# validation
rms_val = sqrt(mean(s1**2))
print(rms_val)
rms_val_equivalent = sqrt(sum(p))
print(rms_val_equivalent)
print(math.isclose(rms_val,rms_val_equivalent))


0.00140477440669
0.00140477440669
True

In [18]:
# print the estimated natural frequency
center_freq = freq_array[lower_tresh+natural_freq_index]
print(center_freq)


1590.70247934

In [19]:
lowcut = center_freq - 100
highcut = center_freq + 100

In [20]:
#filter the signal to s2
s2 = butter_bandpass_filter(s1, lowcut, highcut, samp_freq, order=3)
n2 = len(s2) 
p2 = fft(s2) # take the fourier transform 
n_unique_pts2 = int(ceil((n2+1)/2.0))
p2 = p2[0:n_unique_pts2]
p2 = abs(p2)

In [21]:
# plot the FFT result
p2 = p2 / float(n2) # scale by the number of points so that
# the magnitude does not depend on the length
# of the signal or on its sampling frequency
p2 = p2**2 # square it to get the power
# multiply by two (see technical document for details)
# odd nfft excludes Nyquist point
if n2 % 2 > 0: # we've got odd number of points fft
    p2[1:len(p2)] = p2[1:len(p2)] * 2
else:
    p2[1:len(p2) -1] = p2[1:len(p2) - 1] * 2 # we've got even number of points fft
freq_array2 = arange(0, n_unique_pts2, 1.0) * (samp_freq / n2);
freq_db2 = 10*log10(p2)
lower_tresh2 = int(len(freq_db2)/(samp_freq/2)*800)
higher_tresh2 = int(len(freq_db2)/(samp_freq/2)*4100)
natural_freq_mag2 = max(freq_db2[lower_tresh2:higher_tresh2])
natural_freq_index2 = freq_db2[lower_tresh2:higher_tresh2].argmax()
gca().set_position((.1, 0.18, 1, 1))
plt.plot(freq_array2/1000, freq_db2, '-k')
xlabel('Frequency (kHz)')
ylabel('Power (dB)')
# title(wood_name + ' Wood Longitudinal Wave. Feliana and Ayutyastuti')
title(wood_name + ' Wood Longitudinal Wave. Feliana and Ayutyastuti')
figtext(.02, .02,
        'Estimated Natural Frequency = %0.2f kHz %0.2f dB ' % (freq_array2[lower_tresh2+natural_freq_index2]/1000, natural_freq_mag2))


Out[21]:
<matplotlib.text.Text at 0x9bbf550>

In [22]:
print('Natural Frequency:', natural_freq_mag2)
print('Natural Frequency Index:', natural_freq_index2)


Natural Frequency: -89.9074729927
Natural Frequency Index: 3827

In [23]:
start_point = s2.argmax()
time_array2 = time_array[start_point:start_point+1764]
s2_decay = s2[start_point:start_point+1764]

In [24]:
# plot of the impulse response snapshot
plt.plot(time_array2, s2_decay, color='k')
ylabel('Amplitude')
xlabel('Time (ms)')


Out[24]:
<matplotlib.text.Text at 0x9557f98>

Extract the envelope


In [25]:
# y is the filtered signal with band pass filter 
# y_envel is y-axis value of envelope
# t_envel is x-axis value of envelope
# len is function to calculate the length of a list
y_envel = [s2_decay[data] for data in range(len(s2_decay)-1) if (s2_decay[data-1]<s2_decay[data] and s2_decay[data+1]<s2_decay[data])]
t_envel= [(data+start_point)/44.1 for data in range(len(s2_decay)-1) if (s2_decay[data-1]<s2_decay[data] and s2_decay[data+1]<s2_decay[data])]

In [26]:
y_envel = np.asarray(y_envel)
t_envel = np.asarray(t_envel)
plt.plot(time_array2, s2_decay, color='blue', label="Filtered Spatial-Domain Signal")
plt.plot(t_envel, y_envel, color='k', label="Envelope of the signal")
ylabel('Amplitude')
xlabel('Time (ms)')
plt.legend()


Out[26]:
<matplotlib.legend.Legend at 0xb6d4860>

Performing envelope fitting


In [27]:
data = np.asarray([t_envel, y_envel])
data = np.transpose(data)

Save to file


In [28]:
np.savetxt(wood_name + ".csv", data, delimiter=",")

In [29]:
# function to fit exponent
from math import log
 
def fit_exponent(tList,yList,ySS=0):
   '''
   This function finds a 
       tList in sec
       yList - measurements
       ySS - the steady state value of y
   returns
       amplitude of exponent
       tau - the time constant
   '''
   bList = [log(max(y-ySS,1e-6)) for y in yList]
   b = matrix(bList).T
   rows = [ [1,t] for t in tList]
   A = matrix(rows)
   #w = (pinv(A)*b)
   (w,residuals,rank,sing_vals) = lstsq(A,b)
   tau = -1.0/w[1,0]
   amplitude = exp(w[0,0])
   return (amplitude,tau)

In [30]:
import random

(amplitude_est,tau_est) = fit_exponent(t_envel, y_envel)

In [31]:
print ('Amplitude estimate = %f, tau estimate = %f'
       % (amplitude_est,tau_est))
        
y_est = amplitude_est*(exp(-t_envel/tau_est))


Amplitude estimate = 5555606135631311813568574827606219744608256.000000, tau estimate = 15.758976

In [32]:
print(1/tau_est)


0.0634558999183

Calculating the R-Squared of envelope fitting


In [33]:
residuals = y_est - y_envel
ss_res = np.sum(residuals**2)

In [34]:
ss_tot = np.sum((y_est - np.mean(y_est))**2)

In [35]:
ss_tot2= np.sum((y_envel - np.mean(y_envel))**2)
r_squared2 = 1 - (ss_res / ss_tot2)
print("R-Squared:", r_squared2)


R-Squared: 0.776046645864

In [36]:
plt.plot(t_envel,y_envel, color='k', label="Original envelope")
plt.plot(t_envel, y_est, color='r', label="Envelope fittin result")

# plt.plot(t_envel,y_envel,'k-',t_envel,yEst,'r-')
title(wood_name + ' Wood Envelope Fitting R-square ' + str(r_squared2))
plt.legend()


Out[36]:
<matplotlib.legend.Legend at 0xb7ec2b0>