目标：通过接近人耳朵听力极限的高频部分进行数据的传输
基本思路：用音频靠近20k左右的9个频率来对数据进行编码，该频率出现代表1，没有则为0
流程：
：
编码验证：
#coding=utf-8
import struct
import numpy as np
import bitarray as ba
import getopt
import sys
import wave
import codetable as CT
import cfg as CFG
#encode the serial number into sound wave, author:shuaiwen 405243523@qq.com
sample_rate=44100
sample_width=2 #16bits
def usage():
print('Usage: -s serialnumber(8bytes,avoid use sn:xxxc88cxxx or 8cxxxxxc8) -f output.pcm16')
def print_help_exit():
usage()
sys.exit(-1)
#A:为信号幅值
#fi:为信号频率
#time_s:为时间长度（s）
#sample:为信号采样频率
def signal_xHz(A, fi, time_s, sample,p):
phase_offset=(time_s*fi-int(time_s*fi))/2.0 #由于时间并不能保证完整周期，与其让相位都积累到最后，不如平均一下让相位对称一点可能有助于减少频谱泄漏
#print(phase_offset)
return A*np.sin(np.pi*2*(np.linspace(0,sample*time_s,sample*time_s)*fi/sample+phase_offset)+p)
def signal_xHz2(A, fi, time_s, sample,p):
phase_offset=time_s*fi-int(time_s*fi) #由于时间并不能保证完整周期
pads=int(phase_offset*sample/fi)
signal_parts=A*np.sin(np.pi*2*(np.linspace(0,sample*time_s-pads,int(sample*time_s-pads))*fi/sample)+p)
pad_parts=np.zeros(pads)
return np.append(signal_parts,pad_parts)
def choose_windows(name='Hamming', N=20):
# Rect/Hanning/Hamming
if name == 'Hamming':
window = np.array([0.54 - 0.46 * np.cos(2 * np.pi * n / (N - 1)) for n in range(N)])
elif name == 'Hanning':
window = np.array([0.5 - 0.5 * np.cos(2 * np.pi * n / (N - 1)) for n in range(N)])
elif name == 'Rect':
window = np.ones(N)
elif name == 'Blackman':
window=np.array([0.42 - 0.5 * np.cos(2 * np.pi * n / (N - 1))+0.08*np.cos(4*np.pi*n/(N-1)) for n in range(N)])
elif name == 'Blackman-harris':
window=np.array([0.35875-0.48829*np.cos(2*np.pi*n/(N-1))+0.14128*np.cos(4*np.pi*n/(N-1))-0.01168*np.cos(6*np.pi*n/(N-1))for n in range(N)])
return window
def code2pcm(code):
A=CFG.encode_amplifier #amplitude,volume
duration=CFG.encode_duration_time#720/(408*44.1) #0.040016006402561026 second
#f=[ [400,0.5] ,[412,0.5],[420,0.5],[424,0.5],[432,0.5],[436,0.5],[444,1],[448,1]]
f=CFG.carrier_frq
global sample_rate
data=[]
if len(code)!=len(f):
print("error code length")
exit(1)
data.append(signal_xHz2(A*0.5,CFG.key_freq,duration,sample_rate,0)) #seperate frequency
for i in range(0,len(f)):
if code[i]==True:
data.append(signal_xHz2(A*f[i][1],f[i][0]*44.1,duration,sample_rate,f[i][2]))
[h,w]=np.shape(data)
dataum=np.sum(data,axis=0)/h
#experiment purpose only
# dataum=signal_xHz2(A*0.5,428*sample_rate/1000,duration,sample_rate,0)
# return dataum
#exp end
# win=choose_windows('Blackman',int(sample_rate*duration))
win=choose_windows('Blackman-harris',int(sample_rate*duration))
pcm=np.multiply(dataum,win)
return pcm
def char2pcm(ch):
code=CT.code_table.get(ch.lower(),ba.bitarray('11111111'))
if code.tobytes()[0]==0xff:
return None
return code2pcm(code)
def main():
#guard byte:c(0x5a) w(0xa5)------wccw ,avoid use sn:xxxwccwxxx or wcxxxxxcw
guardbytes='wccw'
filename=''
sn=''
global sample_rate
global sample_width
try:
opts, args = getopt.getopt(sys.argv[1:], "hf:s:")
#print(opts)
for name, value in opts:
if name in ("-h"):
print_help_exit()
elif name in ("-f"):
if (value != None):
filename=value
else:
print('filename wanted')
print_help_exit()
elif name in ("-s"):
sn=value
if len(sn)!=8:
print_help_exit()
else:
print_help_exit()
except getopt.GetoptError:
print_help_exit()
if sn.find(guardbytes)!=-1:
print_help_exit()
dataum=np.zeros(int(CFG.encode_lead_time*sample_rate),dtype=np.int16)
inteval=CFG.encode_interval#0.001
pads=np.zeros(int(inteval*sample_rate),dtype=np.int16)
for i in guardbytes:
print(i)
dataum=np.append(dataum,char2pcm(i))
dataum=np.append(dataum,pads)
for i in sn:
print(i)
dataum=np.append(dataum,char2pcm(i))
dataum=np.append(dataum,pads)
dataum=np.append(dataum,np.zeros(int(0.027*sample_rate),dtype=np.int16))
wf = wave.open(filename, 'wb')
wf.setnchannels(1)
wf.setframerate( sample_rate)
wf.setsampwidth(sample_width)
for x in dataum:
data = struct.pack('H', np.uint16(x))
wf.writeframesraw(data)
wf.close()
print('generating finished')
if __name__ == '__main__':
main()
解码程序：
#coding=utf-8
import numpy as np
import matplotlib.pyplot as plt
import os
import wave
import math
import cv2
import bitarray as ba
import getopt
import sys
import codetable as CT
import cfg as CFG
import difflib
# import jieba
import Levenshtein
import operator
from functools import reduce
#author shuaiwen
#
def usage():
'''
@summary: the usage
@param : None
@return: None
'''
print('Usage: -s serialnumber(8bytes,avoid use sn:xxxc88cxxx or 8cxxxxxc8) -f full_path_to_file.wav')
def print_help_exit():
'''
@summary: print the help
@param : None
@return: None
'''
usage()
sys.exit(-1)
def draw_wav_shape(waveData,nframes,framerate):
#'''绘制语音波形'''
print("plotting signal wave...")
time = np.arange(0,nframes) * (1.0 / framerate)#计算时间
time= np.reshape(time,[nframes,1]).T
plt.plot(time[0,:nframes],waveData[0,:nframes],c="b")
plt.xlabel("time")
plt.ylabel("amplitude")
plt.title("Original wave")
plt.show()
def draw_spectrogram(waveData,framerate,framesize):
print("plotting spectrogram...")
#framelength = 0.02 #帧长20~30ms
#framesize = framelength*framerate #每帧点数 N = t*fs,通常情况下值为256或512,要与NFFT相等\
#而NFFT最好取2的整数次方,即framesize最好取的整数次方
def get_key (dict, value):
return [k for k, v in dict.items() if v == value]
def encodeId2fftId(encode_id):
return int(round(encode_id*CFG.encode_frq_interval/CFG.fft_frq_interval))
def get_auto_corr(timeSeries,k):
'''
Descr:输入：时间序列timeSeries，滞后阶数k
输出：时间序列timeSeries的k阶自相关系数
l：序列timeSeries的长度
timeSeries1，timeSeries2:拆分序列1，拆分序列2
timeSeries_mean:序列timeSeries的均值
timeSeries_var:序列timeSeries的每一项减去均值的平方的和
'''
l = len(timeSeries)
#取出要计算的两个数组
timeSeries1 = timeSeries[0:l-k]
timeSeries2 = timeSeries[k:]
timeSeries_mean = timeSeries.mean()
timeSeries_var = np.array([i**2 for i in timeSeries-timeSeries_mean]).sum()
auto_corr = 0
for i in range(l-k):
temp = (timeSeries1[i]-timeSeries_mean)*(timeSeries2[i]-timeSeries_mean)/timeSeries_var
auto_corr = auto_corr + temp
return auto_corr
def get_up_limit():
return int(math.floor(CFG.freq_up_limit/CFG.fft_frq_interval))
def get_bottom_limit():
return int(math.floor(CFG.freq_down_floor/CFG.fft_frq_interval))#18k 以上的频率
def get_key_frq_id():
id=int(math.floor(CFG.key_freq/CFG.fft_frq_interval)-get_bottom_limit())
return id
#返bottom->up 频率范围内的频谱能量，过滤掉不关心的频率
def get_spectrum(waveData,framerate):
if int(max(abs(waveData)))==0:
return None
waveData = waveData * 1.0/max(abs(waveData))# 归一化
framesize=20*framerate/1000 #how long we will analysis per time
#找到与当前framesize最接近的2的正整数次方
nfftdict = {}
lists = [32,64,128,256,512,1024,2048]
for i in lists:
nfftdict[i] = abs(framesize - i)
sortlist = sorted(nfftdict.items(), key=lambda x: x[1])#按与当前framesize差值升序排列
framesize = int(sortlist[0][0])#取最接近当前framesize的那个2的正整数次方值为新的framesize
NFFT = framesize #NFFT必须与时域的点数framsize相等，即不补零的FFT
overlapSize=int(round(framesize-CFG.step_time*framerate)) #步进距离
print("帧长为{},帧叠为{},傅里叶变换点数为{}".format(framesize,overlapSize,NFFT))
figure,(ax0,ax1) = plt.subplots(nrows=2,figsize=(9,6))
spectrum,freqs,ts,fig = ax0.specgram(waveData,NFFT = NFFT,Fs =framerate,window=np.hanning(M = NFFT),noverlap=overlapSize,mode='default',scale_by_freq=True,sides='default',scale='dB',xextent=None)
upper_limit=get_up_limit()
low_limit=get_bottom_limit()
key_frq_id=get_key_frq_id()
print("display range:{}~{},key frq:{}".format(low_limit,upper_limit,key_frq_id))
#for i in sorted(spectrum[:,0]):#np.argsort(spectrum[:,0]):
# print(i)
#先强制隔离各个频率
#for i in CFG.guard_frq:
# print('clear interval:{}'.format(encodeId2fftId(i)))
# spectrum[encodeId2fftId(i),:]=0
#只截取我们关心的频率范围
spectrum=spectrum[low_limit:upper_limit,:]
#如果关键频率的最大值和平均值的差值太小，则直接丢弃该数据
if (max(spectrum[key_frq_id,:])-np.mean(spectrum[key_frq_id,:]))/np.mean(spectrum[key_frq_id,:]) <1: #值需要摸索
print("unqualifer data,can't handle it")
spectrum=null
return spectrum
def main():
path = "F:\sw\locker"
name = 'rec2.wav'#'menjin8.wav'
filename=''
sn=''
try:
opts, args = getopt.getopt(sys.argv[1:], 'hf:s:')
for name, value in opts:
if name in ("-h"):
print_help_exit()
elif name in ("-f"):
if (value != None):
filename=value
else:
print('filename needed')
print_help_exit()
elif name in ("-s"):
sn=value
else:
print_help_exit()
except getopt.GetoptError:
print_help_exit()
if(filename==''):
filename = os.path.join(path, name)
f = wave.open(filename,'rb')
# 得到语音参数
params = f.getparams()
nchannels, sampwidth, framerate,nframes = params[:4]
audio_data=f.readframes(nframes)
f.close()#关闭文件
wData = np.fromstring(audio_data,dtype=np.short)
#将音频信号规整成每行一路通道信号的格式，即该矩阵一行为一个通道的采样点，共nchannels行
wData = np.reshape(wData,[nframes,nchannels]).T # .T 表示转置,actually,we only use monophony
#analysis 1 second data per time
analysis_len=int(framerate*1.0) #maximum nframes
if analysis_len>nframes:
analysis_len=nframes
key_frq_id=get_key_frq_id()
#载频处理，用于后面的索引
carriers=np.array(CFG.carrier_frq)[:,0]
for i in range(len(carriers)):
carriers[i]=encodeId2fftId(carriers[i])-get_bottom_limit()
carriers=np.append(carriers,key_frq_id)
carriers=carriers.astype(int)
carrier_with_off=[] #载波和其上下两个频率，频率漂移后可能会跑到上下两个位置
for i in carriers:
carrier_with_off.append(i-1)
carrier_with_off.append(i)
carrier_with_off.append(i+1)
#对齐数据下
supplement=analysis_len-len(wData[0,:])%analysis_len
wData=np.c_[wData,np.zeros([len(wData[:,0]),supplement],dtype=np.short)]
offset=analysis_len
start_point=0
waveData=wData[0,start_point:offset]
while len(waveData)>0:
print("start point:{}".format(start_point))
spectrum=get_spectrum(waveData,framerate)
start_point=start_point+offset
waveData=wData[0,start_point:start_point+offset]
if spectrum is None:
continue
#直接把非关键频率的行设置为零，这些值对我们没有用
for i in range(len(spectrum[:,0])):
if i not in carrier_with_off:
spectrum[i,:]=0
#整体放大一下,因为各个频率的强弱不同，按照统一的放大倍数会导致弱的频率会被过滤掉，目前采用以各个频率的最大值来进行放大
for i in carriers:
max_val=max(max(spectrum[i-1,:]),max(spectrum[i,:]),max(spectrum[i+1,:]))
order=1.5#offset,10^order*max大于1
while max_val<1:
max_val=max_val*10
order=order+1
print('order used in frq{} is:{}'.format(i,order))
spectrum[i-1:i+1,:]=np.trunc(spectrum[i-1:i+1,:]*10**order)#过滤太小的数值，这种方法无法抵抗大的波动
#过滤掉少于5个点的列和中间缺失必须出现的频率的列，因为至少需要5个可以观察到的频率（其中包含一个必须出现的频率）需要出现才有意义
#cv2.imshow("spectrum1",np.flipud(spectrum.astype(float)))
#cv2.waitKey(1000000)
#关键频率缺失的列将被设置为0
filter_val=np.mean(spectrum[key_frq_id,:])#np.median(spectrum[key_frq_id,:])#
for idx,val in enumerate(spectrum[key_frq_id,:]):
if val<=filter_val:
#print('set col {} to zero'.format(idx))
spectrum[:,idx]=0
cv2.imshow("split col",np.flipud(spectrum.astype(float)))
cv2.waitKey(1000000)
#过滤掉只有一个值的关键频率点(将该列设置为0)
#row_key=spectrum[key_frq_id,:].astype(int)
#for idx,val in enumerate(row_key):
# if idx!=0 and idx<(len(row_key)-1):
# if row_key[idx-1]==0 and row_key[idx+1]==0:
# spectrum[:,idx]=0
#cv2.imshow("split col2",np.flipud(spectrum.astype(float)))
#cv2.waitKey(1000000)
#以载波频率为中心，将上下两行的值合并到中间（防止频率飘动），对中间值进行局部峰值过滤(低于最大值的一半的值将被设置为0)，找出频率间隔点
for i in carriers:
spectrum[i,:]=spectrum[i-1,:]+spectrum[i,:]+spectrum[i+1,:]
filter_val=int(max(spectrum[i,:])/20)
spectrum[i-1,:]=0
spectrum[i+1,:]=0
for j in range(len(spectrum[i,:])):
if spectrum[i,j]<filter_val:
spectrum[i,j]=0
cv2.imshow("split col3",np.flipud(spectrum.astype(float)))
cv2.waitKey(1000000)
#横向除了载波频率外的值全部设置为零
[h,w]=np.shape(spectrum)
for i in range(h):
if i in carriers:
print("skip:{}".format(i))
else:
spectrum[i,:]=0
cv2.imshow("split col4",np.flipud(spectrum.astype(float)))
cv2.waitKey(1000000)
#移除前面无效的数据列（关键频率不存在的列）
while spectrum[key_frq_id,0]==0:
spectrum=np.delete(spectrum,0,axis=1)
cv2.imshow("split col5",np.flipud(spectrum.astype(float)))
cv2.waitKey(1000000)
#以编码串间隔为依据进行数据分割,有必要移除前导间隔，否则会影响解码分析
minum_blank=int(CFG.encode_lead_time/CFG.step_time) #至少连续这么长的频谱范围内不应该发现主载波即认为是进入了一个字符串分析的前导部分
heads=[]
count=0
for i in range(len(spectrum[key_frq_id,:])):
if spectrum[key_frq_id,i]>0:
if count>=minum_blank:
heads.append(i)
count=0
else:
count=count+1
print("split into {} unit".format(len(heads)))
s=0
for i in heads:
text="analysis {}".format(i)
print(text)
keys=[]
# d=spectrum[:,s:i]
# cv2.imshow(text,np.flipud(d.astype(float)))
# cv2.waitKey(1000000)
raw=analysis_units(spectrum[:,s:i],key_frq_id,carriers,s)
str=raw2str(raw)
#print("raw code:{},parsed code:{}".format(i.tobytes().hex(),get_key(CT.code_table,i)),end=' ')
ratio=get_sim(str,sn)
print("parse result:{},Levenshtein similarity::{}".format(str,ratio))
s=i
raw=analysis_units(spectrum[:,s:],key_frq_id,carriers,s)
for j in raw:
print("raw str:{}".format(get_bits(j)))
str=raw2str(raw)
ratio=get_sim(str,'wccw'+sn)
print ('difflib similarity1:{} target:{} parsed:{}'.format(ratio,sn,str))
#尾巴数据会遗漏，要分析一下
def get_bits(bitarr):
s=" "
for m in range(0,bitarr.length()):
if bitarr[m]== True:
s=s+"1"
else:
s=s+"0"
return s
def get_sim(str1,str2):
seq = difflib.SequenceMatcher(None, str1,str2)
ratio = seq.ratio()
return ratio
def raw2str(raw):
keys=[]
s=''
for n in raw:
keys.append(get_key(CT.code_table,n))
s=''.join(reduce(operator.add, keys)) #ToDo, 如果keys全部为数字，会变成数字相加的合，需要注意
return s
#对分割的单元数据进行解析
def analysis_units(spectrum,key_frq_id,carriers,seq=0):
#移除前后无效的数据列（关键频率不存在的列）
while spectrum[key_frq_id,0]==0:
spectrum=np.delete(spectrum,0,axis=1)
while spectrum[key_frq_id,-1]==0:
spectrum=np.delete(spectrum,-1,axis=1)
cv2.imshow("a{}".format(seq),np.flipud(spectrum.astype(float)))
cv2.waitKey(1000000)
#纵向以关键频率为线索进行纵向数值合并，同时找到间距点（周期）,通过自相关进行查找,
[h,w]=np.shape(spectrum)
if w <8:
print("skip unit analysis because data isn't enough")
return -1
corr=[]
for i in range(6): #自相关间距不能设置太大，不然可能会重复大的周期，现在这个只是经验值
corr.append(get_auto_corr(spectrum[key_frq_id,:],i+1))
max_index = corr.index(max(corr))+1
#纵向可以以max_index(间距点) 进行合并，然后将剩下的值设置为0
#凑个整数列，数据不够0来凑
supplement=max_index-len(spectrum[key_frq_id,:])%max_index
spectrum=np.c_[spectrum,np.zeros([len(spectrum[:,0]),supplement],dtype=np.int)]
key_pos_col=[] #用来保存有效数据列位置，和carriers搭配形成索引点下标
key_pos_col.append(0)
idx=0
while idx< len(spectrum[key_frq_id,:])-1:
i=1
while int(spectrum[key_frq_id,idx+i])>0: #合并连续的数
spectrum[:,idx]=spectrum[:,idx]+spectrum[:,idx+i]
spectrum[:,idx+i]=0
i=i+1
idx=idx+i
while int(spectrum[key_frq_id,idx]) ==0: #跳过0
idx=idx+1
if idx>=len(spectrum[key_frq_id,:]):
break
if idx<len(spectrum[key_frq_id,:]):
key_pos_col.append(idx)
# for i in range(int(len(spectrum[key_frq_id,:])/max_index)):
# key_pos_col.append(i*max_index)
# spectrum[:,i*max_index]=spectrum[:,i*max_index]+spectrum[:,i*max_index+1]+spectrum[:,i*max_index+2]+spectrum[:,i*max_index+3]
# spectrum[:,i*max_index+1]=spectrum[:,i*max_index+2]=spectrum[:,i*max_index+3]=0
cv2.imshow("a1{}".format(seq),np.flipud(spectrum.astype(float)))
cv2.waitKey(1000000)
raw=[]
for i in key_pos_col:
dataum=ba.bitarray()
for j in carriers:
if j!=key_frq_id: #这个频率点不是有效数据，要剔除
if int(spectrum[j][i])>0:
print(spectrum[j][i])
dataum.append(True)
else:
dataum.append(False)
raw.append(dataum)
return raw
if sys.version_info<(3,0):
print ('<----parse result:',)
else:
print('<----parse result:',end=' ')
for i in raw:
print("raw code:{},parsed code:{}".format(i.tobytes().hex(),get_key(CT.code_table,i)),end=' ')
#print(dataum.tobytes().hex(),end=' ')
#print(get_key(CT.code_table,dataum),end=' ')
print('---->')
# cv2.imshow("spectrum6",spectrum.astype(float))
# cv2.waitKey(1000000)
if __name__ == '__main__':
main()
2352
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
