<jupyter_start><jupyter_text>Analysing ultrasound beacons in the wild<jupyter_code>from scipy.io.wavfile import read as wavread
import matplotlib.pyplot as plt
import numpy as np

from scipy.signal import hilbert, chirp
from scipy.signal import argrelextrema<jupyter_output><empty_output><jupyter_text>Let's open the signal<jupyter_code>audiofilename = "starbuckmontpar.wav"
[samplerate, x] = wavread(audiofilename) # x is a numpy array of integer, representing the samples 
# scale to -1.0 -- 1.0
if x.dtype == 'int16':
    nb_bits = 16 # -> 16-bit wav files
elif x.dtype == 'int32':
    nb_bits = 32 # -> 32-bit wav files
max_nb_bit = float(2 ** (nb_bits - 1))
samples = x / (max_nb_bit + 1.0) # samples is a numpy array of float representing the samples<jupyter_output><empty_output><jupyter_text>Let's plot the signal<jupyter_code>t = range(len(x))
for k in range(len(x)):
    t[k] = 1.0*t[k] / samplerate
plt.plot(t,samples)  # let's plot the sound
plt.xlabel("Time in s")
plt.show()<jupyter_output><empty_output><jupyter_text>FFT analysis<jupyter_code>L = len(x)

f1 = int(19500*L/samplerate)
f2 = int(20500*L/samplerate)

f = range(len(x))
for k in range(len(x)):
    f[k] = 1.0e-3*f[k]*samplerate/len(x)
# Let's get the signal around 20khz:
A = np.fft.fft(samples)


plt.plot(f[f1:f2],np.abs(A)[f1:f2])  
plt.xlabel("Freq in kHz")
plt.show()<jupyter_output><empty_output><jupyter_text>We observe a peak at 20khz, with symetrical peaks around - right what we were looking for !<jupyter_code>T = 4*1024
d = 200000
nS = int(L/T)

d1 = int(19500*T/samplerate)
d2 = int(20500*T/samplerate)

B = np.fft.fft(samples[d:d+T])
plt.plot(np.abs(B[d1:d2]))   
plt.show()
SP = []
for i in range(nS):
    B = np.fft.fft(samples[i*T:(i+1)*T-1])
    SP.append(np.abs(B[d1:d2]))  
SPP = np.asarray(SP)
m1 = 0
m2 = L
xl = range(d2-d1)
yl = range(m2-m1)
plt.figure(figsize=(15,5))

plt.imshow(np.transpose(SPP[m1:m2]), cmap='hot',vmin=SPP.min(),vmax=SPP.max()) 
plt.title("Plot of the 19.5khz to 20.5khz ")
plt.show()<jupyter_output><empty_output><jupyter_text>* Seems like peaks at 20khz, with ripples at +- 250hz* See https://www.sec.cs.tu-bs.de/pubs/2016-batmobile.pdf<jupyter_code>C = A
for k in range(len(A)/2+1):
    if k < int(0.9*f1) or k > int(1.1*f2):
        C[k] = 0
        C[-k] = 0

Cleaned = np.real( np.fft.ifft(C) )
plt.plot(t,samples) 
plt.plot(t,Cleaned,"r")   
plt.show()
analytic_signal = hilbert(Cleaned[0:4096/2])
amplitude_envelope = np.abs(analytic_signal)

plt.plot(t[0:4096/10],Cleaned[0:4096/10],"r")   
plt.plot(t[0:4096/10],amplitude_envelope[0:4096/10],"b")   
plt.show()
print 1.0/250/2

D = argrelextrema(amplitude_envelope[0:4096/10], np.less)
O = []
for k in D[0]:
    O.append(t[k])
for k in range(len(O)-1):
    tmp = O[k+1]-O[k]
    if tmp > 1e-4: 
        print tmp<jupyter_output>0.000541666666667
0.000520833333333
0.001375
0.0004375
0.000604166666667
0.000541666666667
0.000104166666667
0.0014375
0.00127083333333<jupyter_text>Seing a 250hz modulation ?<jupyter_code>m = max(amplitude_envelope)
 
plt.plot(t[0:4096/2],(8*amplitude_envelope/m).astype(int),"b")   
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1662412252712189, 'fraction_numerical': 0.09285258455647734, 'mean_word_length': 4.773480662983426, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 32, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 6, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '194176', 'license': 'MIT', 'n_tokens_mistral': 1614, 'n_tokens_neox': 1413, 'n_words': 263, 'path': 'mobile/20180110-SMontparnasseBeacon.ipynb', 'repo_name': 'kelu124/ultrasoundhunter', 'size': '3134'}