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"
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]:
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])
In [9]:
try:
s1 = snd[:,0]
except IndexError:
s1 = snd
In [10]:
s1 # original signal
Out[10]:
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]:
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]
In [14]:
# show the impulse response snapshot
plt.plot(start_time_array, start_s1, color='k')
ylabel('Amplitude')
xlabel('Time (ms)')
Out[14]:
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)
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]:
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))
In [18]:
# print the estimated natural frequency
center_freq = freq_array[lower_tresh+natural_freq_index]
print(center_freq)
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]:
In [22]:
print('Natural Frequency:', natural_freq_mag2)
print('Natural Frequency Index:', natural_freq_index2)
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]:
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]:
In [27]:
data = np.asarray([t_envel, y_envel])
data = np.transpose(data)
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))
In [32]:
print(1/tau_est)
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)
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]: