本文介绍了一种基于傅里叶变换和经验模态分解（EMD）的单频检测技术，该技术能有效识别和提取特定频率信号。通过计算声波主频并应用于EMD变换中，保留目标频率成分的同时去除其他干扰频率。实验证明，该方法在低频段信号处理方面表现良好。
目录
- 一、背景
- 二、原理
- 三、实现
- 四、代码
- 五、总结
一、背景
对应实际问题：使用单频检测水管中阀门或其他组件。
二、原理
通过傅里叶变换计算声波的主频，将得到的主频运用于emd变换中，保留基本的主频，去除其他频率，再组合形成波形。
三、实现
傅里叶变换和emd变化python均有现有库，分别在numpy和PyEMD，值得注意的是pyemd这个库是地球移动距离(earth move distance)而非经验模态分析(Empirical Mode Decomposition)，PyEMD对应的叫emd-signal。
因此具体实现只需对PyEMD的Visualization部分进行改编，新建类继承Visualization类，并对plot_instant_freq方法进行重写，使该方法能够显示频率线。
并利用原有的私有函数_calc_inst_freq计算每组频率值，对得到的频率去除离散点后取平均值与傅里叶得到的频率值比较，得到应保留的频率。再累加输出波形，得到结果。
四、代码
代码如下
(1)主文件
import wave
import pyaudio
import numpy
import matplotlib.pylab as pylab
from matplotlib import pyplot
from pyhht.emd import EMD
from pyhht.visualization import plot_imfs
from PyEMD import EMD,Visualisation,EEMD
from EmdVisFreq import EmdVisFreq
from statsmodels.tsa.stattools import grangercausalitytests
import pandas
def ReadWav(filename):
wf = wave.open(filename, "rb")
#创建PyAudio对象
p = pyaudio.PyAudio()
stream = p.open(format=p.get_format_from_width(wf.getsampwidth()),
channels=wf.getnchannels(),
rate=wf.getframerate(),
output=True)
nframes = wf.getnframes()
framerate = wf.getframerate()
#读取完整的帧数据到str_data中，这是一个string类型的数据
str_data = wf.readframes(nframes)
wf.close()
#将波形数据转换为数组
# A new 1-D array initialized from raw binary or text data in a string.
wave_data = numpy.fromstring(str_data, dtype=numpy.short)
#将wave_data数组改为2列，行数自动匹配。在修改shape的属性时，需使得数组的总长度不变。
#单声道时为-1,1
wave_data.shape = -1,1
#将数组转置
wave_data = wave_data.T
#time 也是一个数组，与wave_data[0]或wave_data[1]配对形成系列点坐标
time = numpy.arange(0,nframes)*(1.0/framerate)
return wave_data,time,framerate
def liner_correlation(wave0,wave1,templen=44100):
pylab.plot(wave0[0][:templen],wave1[0][:templen])
pylab.show()
def drawWave(time,wave,templen=44100,figname=""):
time=time[:templen]
wave_temp=wave[0][:templen]
pylab.figure(figsize=(16,9),dpi=600)
pylab.plot(time, wave_temp,label="Wave")
pylab.xlabel("Time/seconds)")
pylab.title(figname[20:32]+"FourierTransform")
pylab.ylabel('Range')
#pylab.set(title=figname+"FourierTransform", xlabel='Spectrum/Hz',ylabel='Range')
#ax.set_xticklabels([int(item-timerange) for item in ax.get_xticks()])
pylab.legend()
pylab.savefig(figname+"_"+"Wave.png",format="png")
pylab.close()
#pylab.show()
def fourierTransform(wave_data,framerate,figname="",originalSpectrum=""):
N=44100#采样率
start=0 #开始采样位置
df = framerate/(N-1) # 分辨率
freq = [df*n for n in range(0,N)] #N个元素
wave_data2=wave_data[0][start:start+N]
#wave_data2=wave_data[start:start+N]
spectrum=numpy.fft.fft(wave_data2)*2/N
#常规显示采样频率一半的频谱
d=int(len(spectrum)/2)
#print(d)
#仅显示频率在2000以下的频谱
while freq[d]>3000:
d-=10
#d=halfd(spectrum,freq)
#print(d)
pylab.figure(figsize=(16,9),dpi=600)
pylab.plot(freq[:d-1],abs(spectrum[:d-1]),'r',label='Fourier Spectrum')
pylab.axvline(originalSpectrum,linestyle='--',label='Original Spectrum')
pylab.title(figname[20:32]+"FourierTransform")
pylab.xlabel('Spectrum/Hz')
pylab.xticks(fontsize=10.5)
pylab.ylabel('Range')
#pylab.set(title=figname+"FourierTransform", xlabel='Spectrum/Hz',ylabel='Range')
#ax.set_xticklabels([int(item-timerange) for item in ax.get_xticks()])
pylab.legend()
pylab.savefig(figname+"_"+"Fourier.png",format="png")
pylab.close()
#pylab.show()
return spectrum,freq
def halfd(spectrum,freq):
d=int(len(spectrum)/2)
#仅显示频率在2000以下的频谱
while freq[d]>2000:
d-=10
return d
def range_filter(spectrum,freq):
for index in range(len(spectrum)):
#print(n)
if abs(spectrum[index])<0:
spectrum[index]=0
#print(n)
#print(c)
d=halfd(spectrum,freq)
pylab.plot(freq[:d-1],abs(spectrum[:d-1]),'r')
pylab.show()
return spectrum
def freq_filter(spectrum,freq):
for index in range(len(spectrum)):
#print(n)
if index<400 or index>700:
spectrum[index]=0
#print(c)
d=halfd(spectrum,freq)
pylab.plot(freq[:d-1],abs(spectrum[:d-1]),'r')
pylab.show()
return spectrum
def ifftTransform(spectrum,time):
N=44100
f_wave=numpy.fft.ifft(spectrum*N/2)
pylab.plot(time[:44100], f_wave)
pylab.xlabel("time (seconds)")
pylab.show()
return f_wave
def emdTransform(wave,time,figname="",originalSpectrum=""):
emd = EMD()
emd.emd(wave)
imfs, res = emd.get_imfs_and_residue()
#print(imfs.shape)
#print(time.shape)
# In general:
#components = EEMD()(S)
#imfs, res = components[:-1], components[-1]
#print(imfs[0][0])
#print(type(imfs[0][0]))
'''
vis = Visualisation()
vis.plot_imfs(imfs=imfs, residue=res, t=time, include_residue=True)
vis.plot_instant_freq(time, imfs=imfs)
vis.show()
'''
evf=EmdVisFreq()
#freqs=evf.getfreq(imfs)
print('evf')
evf.plot_imfs(imfs=imfs, residue=res, t=time, include_residue=True)
evf.savefig(figname+"_"+"EMDRange.png")
evf.myplot_instant_freq(originalSpectrum,time, imfs=imfs)
evf.savefig(figname+"_"+"EMDSpectrum.png")
return imfs
def emdTransformNoPic(wave,time):
emd = EMD()
emd.emd(wave)
imfs, res = emd.get_imfs_and_residue()
return imfs
def eemdTransform(wave,time):
eemd=EEMD()
eemd.trials=50
eemd.noise_seed(12345)
eemd.eemd(wave,time,10)
imfs,res=eemd.get_imfs_and_residue()
vis = Visualisation()
vis.plot_imfs(imfs=imfs, residue=res, t=time, include_residue=True)
vis.plot_instant_freq(time, imfs=imfs)
vis.show()
return imfs
def emdSelect(imfs,start=2,stop=4):
result=numpy.zeros([len(imfs[0])],float)
for i in range(start,stop):
result=result+imfs[i]
return result
def emdSelectFreq(imfs,time,originalSpectrum=""):
evf=EmdVisFreq()
freqs=evf.getfreq(imfs,time)
result=numpy.zeros([len(imfs[0])],float)
x=0
y=0
for i in freqs:
i05=numpy.percentile(i,5)
i95=numpy.percentile(i,95)
meani=numpy.mean(i)
for j in i:
if j<i05 or j>i95:
j=meani
#print(len(i))
if abs(meani-originalSpectrum)<500:
result=result+imfs[x]
#print(x)
x=x+1
return result
def emdFilter(filename,templen=500,start=2,stop=4):
wave,time,framerate=ReadWav(filename)
_wave=wave[0][:templen]
_time=time[:templen]
imfs=emdTransformNoPic(_wave,_time)
_imfs=emdSelect(imfs)
return _imfs
def emdFilterWave(wave,time,templen=500,originalSpectrum="",figname=""):
#wave,time,framerate=ReadWav(filename)
#_wave=wave[0][:templen]
#_time=time[:templen]
_wave=wave
_time=time
imfs=emdTransformNoPic(_wave,_time)
_imfs=emdSelectFreq(imfs,_time,originalSpectrum)
#橙线
pylab.figure(figsize=(16,9),dpi=600)
pylab.title(figname[20:32]+"EMDFilter")
pylab.xlabel('Time/s')
pylab.ylabel('Range')
pylab.plot(time,_wave,label="Original Wave")
pylab.plot(time,_imfs,label="EMDFIlter Wave")
pylab.legend()
pylab.savefig(figname+"_"+"EMDFilter.png",format="png")
pylab.close()
return _imfs
def timeGrangeCorrelation(filename1,filename2):
wave1,time1,framerate1=ReadWav(filename1)
wave2,time2,framerate2=ReadWav(filename2)
templen=500
_wave1=wave1[0][:templen]
_time1=time1[:templen]
_wave2=wave2[0][:templen]
_time2=time2[:templen]
imfs1=emdTransformNoPic(_wave1,_time1)
imfs2=emdTransformNoPic(_wave2,_time2)
start=1
stop=3
_imfs1=emdSelect(imfs1,start,stop)
_imfs2=emdSelect(imfs2,start,stop)
Temp=numpy.empty([2,len(_imfs1)],float)
Temp[0]=_imfs1
Temp[1]=_imfs2
Temp=Temp.T
grangercausalitytests(Temp,maxlag=100)
def crosscorr(datax, datay, lag=0, wrap=False):
""" Lag-N cross correlation.
Shifted data filled with NaNs
Parameters
----------
lag : int, default 0
datax, datay : pandas.Series objects of equal length
Returns
----------
crosscorr : float
"""
if wrap:
shiftedy = datay.shift(lag)
shiftedy.iloc[:lag] = datay.iloc[-lag:].values
return datax.corr(shiftedy)
else:
shiftedy = datay.shift(lag)
shiftedy=shiftedy.fillna(0)
return datax.corr(shiftedy)
def timeLateCorrelation(filename1,filename2,timerange=500):
_wave1=emdFilter(filename1)
_wave2=emdFilter(filename2)
_wave1df=pandas.DataFrame(_wave1)
_wave2df=pandas.DataFrame(_wave2)
_wave1df=_wave1df[0]
_wave2df=_wave2df[0]
rs = [crosscorr(_wave1df,_wave2df, lag) for lag in range(-int(timerange-1),int(timerange))]
print(rs)
offset = numpy.ceil(len(rs)/2)-numpy.argmax(rs)
f,ax=pyplot.subplots(figsize=(16,9),dpi=600)
ax.plot(rs)
ax.axvline(numpy.ceil(len(rs)/2),color='k',linestyle='--',label='Center')
ax.axvline(numpy.argmax(rs),color='r',linestyle='--',label='Peak synchrony')
ax.set(title=f'Offset = {offset} framesnS1 leads <> S2 leads',ylim=[0,.31],xlim=[0,timerange*2], xlabel='Offset',ylabel='Pearson r')
ax.set_xticklabels([int(item-timerange) for item in ax.get_xticks()])
pyplot.legend()
#pyplot.figure(figsize=(16,9),dpi=600)
pyplot.savefig(filename1+"To"+"TLCC.png",format="png")
#pyplot.show()
rsmax=numpy.argmax(rs)
pylab.figure(figsize=(16,9),dpi=600)
pylab.plot(range(0,timerange),_wave2df.shift(-timerange+rsmax).fillna(0))
pylab.plot(range(0,timerange),_wave1df)
pylab.savefig(filename1+"To"+"MaxCorr.png",format="png")
#pylab.show()
if __name__ == "__main__":
path="D:WorksoundsoundALL"
spectures=numpy.array([0,100,200,300,400,500,600,1000,2000,3000])
plots=range(1,5)
openrates=["20%"]
filetype=".wav"
filename=""
templen=500
_damping=range(0,len(plots))
daming=numpy.array(_damping)
for specture in spectures:
_tempdampi=0
damping0=0
for plot in plots:
for openrate in openrates:
filename=path+str(specture)+"hz"+"-"+str(plot)+"-"+openrate+filetype
_wave,_time,_framerate=ReadWav(filename)
drawWave(_time,_wave,templen=templen,figname=filename)
fourierTransform(_wave,_framerate,filename,specture)
wave0=_wave[0][:templen]
time0=_time[:templen]
emdTransform(wave0,time0,filename,specture)
tempimfs=emdFilterWave(wave0,time0,originalSpectrum=specture,figname=filename)
filename=path+str(specture)+"hz"+"-"+str(plot)+"-"+openrate+filetype
#_wave,_time,_framerate=ReadWav(filename)
#fourierTransform(_wave,_framerate,filename)
if _tempdampi==0:
daming[_tempdampi]=numpy.mean(numpy.abs(tempimfs))
damping0=numpy.mean(numpy.abs(tempimfs))
print("mean0:")
print(numpy.mean(numpy.abs(tempimfs)))
daming[_tempdampi]=numpy.mean(numpy.abs(tempimfs))/damping0
print("mean:")
print(numpy.mean(numpy.abs(tempimfs)))
print("danming")
print(daming[_tempdampi])
_tempdampi=_tempdampi+1
print("tempdmi:")
print(_tempdampi)
print("1")
pyplot.plot(range(1,len(plots)+1),daming)
print("dam:")
print(daming)
pyplot.savefig(path+str(specture)+"hz"+"-"+"-")
print("2")（2）Visualization子类
from PyEMD import Visualisation
import pylab as plt
import numpy as np
class EmdVisFreq(Visualisation):
PLOT_WIDTH = 18
PLOT_HEIGHT_PER_IMF = 4.5
def myplot_instant_freq(self,mark, t, imfs=None, order=False, alpha=None):
if alpha is not None:
assert -0.5 < alpha < 0.5 , '`alpha` must be in between -0.5 and 0.5'
imfs, _ = self._check_imfs(imfs, None, False)
num_rows = imfs.shape[0]
imfs_inst_freqs = self._calc_inst_freq(imfs, t, order=order, alpha=alpha)
fig, axes = plt.subplots(num_rows, 1, figsize=(self.PLOT_WIDTH, num_rows*self.PLOT_HEIGHT_PER_IMF),dpi=600)
if num_rows == 1:
axes = fig.axes
axes[0].set_title("Instantaneous frequency")
for num, imf_inst_freq in enumerate(imfs_inst_freqs):
ax = axes[num]
ax.plot(t, imf_inst_freq)
marking= np.ones([len(imf_inst_freq),], dtype = int)
marking=marking*mark
ax.plot(t,marking)
ax.set_ylabel("IMF {} [Hz]".format(num+1))
# Making the layout a bit more pleasant to the eye
plt.tight_layout()
def plot_imfs(self, imfs=None, residue=None, t=None, include_residue=True):
"""Plots and shows all IMFs.
All parameters are optional since the `emd` object could have been passed when instantiating this object.
The residual is an optional and can be excluded by setting `include_residue=False`.
"""
imfs, residue = self._check_imfs(imfs, residue, include_residue)
num_rows, t_length = imfs.shape
num_rows += include_residue is True
t = t if t is not None else range(t_length)
fig, axes = plt.subplots(num_rows, 1, figsize=(self.PLOT_WIDTH, num_rows*self.PLOT_HEIGHT_PER_IMF))
if num_rows == 1:
axes = list(axes)
axes[0].set_title("Time series")
for num, imf in enumerate(imfs):
ax = axes[num]
ax.plot(t, imf)
ax.set_ylabel("IMF " + str(num+1))
if include_residue:
ax = axes[-1]
ax.plot(t, residue)
ax.set_ylabel("Res")
# Making the layout a bit more pleasant to the eye
plt.tight_layout()
def savefig(self,figname="",):
plt.savefig(figname+"_"+"Fourier.png",format="png")
plt.close()
def getfreq(self,imfs,time):
t=time
freqs = self._calc_inst_freq(imfs, t, order=False, alpha=None)
return freqs五、总结
对低频段处理较好，高频段因为噪音也主要集中于高频段，滤波效果差。
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
