python怎么实现eemd_使用Fourier+EMD水管中传播的单频音波进行滤波（Python）

布满荆棘的人生 2023-01-03 01:53 57阅读 0赞

### 目录 ###

*  一、背景
 *  二、原理
 *  三、实现
 *  四、代码
 *  五、总结

# 一、背景 #

对应实际问题：使用单频检测水管中阀门或其他组件。

# 二、原理 #

通过傅里叶变换计算声波的主频，将得到的主频运用于emd变换中，保留基本的主频，去除其他频率，再组合形成波形。

# 三、实现 #

傅里叶变换和emd变化python均有现有库，分别在numpy和PyEMD，值得注意的是pyemd这个库是地球移动距离(earth move distance)而非经验模态分析(Empirical Mode Decomposition)，PyEMD对应的叫emd-signal。  
因此具体实现只需对PyEMD的Visualization部分进行改编，新建类继承Visualization类，并对plot\_instant\_freq方法进行重写，使该方法能够显示频率线。![6a7902d6353ec0debe3b3af05185bbad.png][]

![fcf320f6c6a960160b96477f05a40b96.png][]  
并利用原有的私有函数\_calc\_inst\_freq计算每组频率值，对得到的频率去除离散点后取平均值与傅里叶得到的频率值比较，得到应保留的频率。再累加输出波形，得到结果。![f55844478209cd8426126828a79c073b.png][]

# 四、代码 #

代码如下  
(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 ji95:
                    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} frames\nS1 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:\\Work\\sound\\sound\\ALL\\"
        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

# 五、总结 #

对低频段处理较好，高频段因为噪音也主要集中于高频段，滤波效果差。

[6a7902d6353ec0debe3b3af05185bbad.png]: /images/20221119/0c78eef1f77c42dd83302e49bcd28b25.png
[fcf320f6c6a960160b96477f05a40b96.png]: /images/20221119/106e898ef8c84fdeb13fe8f7e5d54aa4.png
[f55844478209cd8426126828a79c073b.png]: /images/20221119/7523fa6122ee4c36bf9a1f3bfeab655d.png