<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 The instantaneous phase, amplitude, and frequency of a sawtooth wave
can changing the smoothness of a signal change the instantaneous freq dynamics? (i.e. more sharp transitions)
Step1:  1. Simulate a 20Hz sawtooth wave
Step2:  2. Filter in  13-30Hz band
Filter is very short (like steph said in email)
Step3:  3. Calculate instantaneous measures
Step4:  4. Visualize
Step5:  5. Example with real data
Step6:  6. Plot relationship between inst. amp. and inst. freq.
<ASSISTANT_TASK:>
Python Code:
%config InlineBackend.figure_format = 'retina'
%matplotlib inline
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style('white')
Explanation: The instantaneous phase, amplitude, and frequency of a sawtooth wave
can changing the smoothness of a signal change the instantaneous freq dynamics? (i.e. more sharp transitions)
End of explanation
# Define sawtooth shape in some number of samples
x1 = np.array([0,.05,.2,1,.9,.8,.7,.6,.5,.4,.3,.2,.1,.05,.01])
t1 = np.arange(len(x1))
# Interpolate sawtooth so it has 50 samples (50ms = 20Hz wave)
from scipy import interpolate
f = interpolate.interp1d(t1, x1)
t2 = np.linspace(0,len(t1)-1,50)
x2 = f(t2)
# Tile the new sawtooth to last 5 seconds
x = np.tile(x2,100)
x = x - np.mean(x)
Fs = 1000.
t = np.arange(0,5,.001)
# Plot sawtooth
plt.figure(figsize=(5,2))
plt.plot(t, x)
plt.ylim((-.7,.7))
plt.xlim((0,.5))
plt.xlabel('Time (s)')
plt.ylabel('Voltage (a.u.)')
Explanation: 1. Simulate a 20Hz sawtooth wave
End of explanation
from misshapen import nonshape
x_filt, _ = nonshape.bandpass_default(x, (13,30),Fs, w=1.5,rmv_edge=False)
Explanation: 2. Filter in  13-30Hz band
Filter is very short (like steph said in email)
End of explanation
x_amp = np.abs(sp.signal.hilbert(x_filt))
x_phase = np.angle(sp.signal.hilbert(x_filt))
# Instantaneous freq
x_freq = np.diff(x_phase)
x_freq[x_freq<0] = x_freq[x_freq<0]+2*np.pi
x_freq = x_freq*Fs/(2*np.pi)
Explanation: 3. Calculate instantaneous measures
End of explanation
samp_plt = range(1000,1200)
plt.figure(figsize=(5,5))
plt.subplot(3,1,1)
plt.plot(t[samp_plt],x[samp_plt],'k')
plt.plot(t[samp_plt],x_filt[samp_plt],'r')
plt.plot(t[samp_plt],x_amp[samp_plt],'b')
plt.xlim((t[samp_plt][0],t[samp_plt][-1]))
plt.xticks([])
plt.ylabel('raw (black)\nfiltered (red)\nInst. Amp. (blue)')
plt.subplot(3,1,2)
plt.plot(t[samp_plt],x_phase[samp_plt],'k')
plt.xlim((t[samp_plt][0],t[samp_plt][-1]))
plt.xticks([])
plt.ylabel('Inst. Phase (rad)')
plt.subplot(3,1,3)
plt.plot(t[samp_plt],x_freq[samp_plt],'k')
plt.xlim((t[samp_plt][0],t[samp_plt][-1]))
plt.xlabel('Time (s)')
plt.ylabel('Inst. Freq. (Hz)')
Explanation: 4. Visualize
End of explanation
x2 = np.load('C:/gh/bv/misshapen/exampledata.npy')
from misshapen import nonshape
x2_filt, _ = nonshape.bandpass_default(x2, (13,30),Fs, w=3,rmv_edge=False)
x_amp = np.abs(sp.signal.hilbert(x2_filt))
x_phase = np.angle(sp.signal.hilbert(x2_filt))
# Instantaneous freq
x_freq = np.diff(x_phase)
x_freq[x_freq<0] = x_freq[x_freq<0]+2*np.pi
x_freq = x_freq*Fs/(2*np.pi)
samp_plt = range(3000,5000)
plt.figure(figsize=(10,6))
plt.subplot(4,1,1)
plt.plot(t[samp_plt],x2[samp_plt],'k')
plt.plot(t[samp_plt],x2_filt[samp_plt],'r')
plt.plot(t[samp_plt],x_amp[samp_plt],'b')
plt.xlim((t[samp_plt][0],t[samp_plt][-1]))
plt.xticks([])
plt.ylabel('raw (black)\nfiltered (red)\nInst. Amp. (blue)')
plt.subplot(4,1,2)
plt.plot(t[samp_plt],x_phase[samp_plt],'k')
plt.xlim((t[samp_plt][0],t[samp_plt][-1]))
plt.xticks([])
plt.ylabel('Inst. Phase (rad)')
plt.subplot(4,1,3)
plt.plot(t[samp_plt],x_freq[samp_plt],'k')
plt.xlim((t[samp_plt][0],t[samp_plt][-1]))
plt.xticks([])
plt.ylim((0,30))
plt.ylabel('Inst. Freq. (Hz)')
plt.subplot(4,1,4)
plt.plot(t[samp_plt],x_amp[samp_plt],'k')
plt.ylabel('Inst. Amp.')
plt.xlabel('Time (s)')
ax = plt.gca()
ax2 = ax.twinx()
ax2.plot(t[samp_plt],x_freq[samp_plt],'r')
plt.xlim((t[samp_plt][0],t[samp_plt][-1]))
plt.ylim((0,30))
plt.ylabel('Inst. Freq. (Hz)',color='r')
Explanation: 5. Example with real data
End of explanation
plt.figure(figsize=(6,6))
plt.plot(x_freq, x_amp[1:],'k.',alpha=.01)
plt.xlim((0,100))
plt.xlabel('Inst. Freq. (Hz)')
plt.ylabel('Inst. Amp. (uV)')
Explanation: 6. Plot relationship between inst. amp. and inst. freq.
End of explanation
<END_TASK>