python做50Hz陷波器的应用实例

仅作为记录,大佬请跳过。

使用说明——可直接运行

将fftandshow.py、mysubplot.py和notch_implement.py文件放在同一个文件夹里,运行notch_implement.py,即可显示

直接上代码

notch_implement.py文件:

import numpy as np
import mysubplot
import fftandshow
import matplotlib.pyplot as plt
from scipy.signal import iirfilter,lfilter

def Implement_Notch_Filter(time, band, freq, order, filter_type, data):
    from scipy.signal import iirfilter
    fs   = 1/time
    nyq  = fs/2.0
    low  = freq - band/2.0
    high = freq + band/2.0
    low  = low/nyq
    high = high/nyq
    b, a = iirfilter(order, [low, high],  btype='bandstop',
                     analog=False, ftype=filter_type)
    filtered_data = lfilter(b, a, data)
    return filtered_data

ts=0.005
t=np.arange(0,1,ts)
a=np.sin(2*np.pi*80*t)+np.sin(2*np.pi*50*t)+np.sin(2*np.pi*20*t);data=a

fftandshow.fftspectrum(data,1/ts);plt.show()
data_filt=Implement_Notch_Filter(ts,6,50,2,'butter',data)
fftandshow.fftspectrum(data_filt,1/ts);plt.show()

f1,s1=fftandshow.fftspectrum_value(data,1/ts);f2,s2=fftandshow.fftspectrum_value(data_filt,1/ts)
mysubplot.mysubplot_addx([[t,data],[t,data_filt],[f1,s1],[f2,s2]])

fftandshow.py文件:

import numpy as np
import matplotlib.pyplot as plt

def fftspectrum(thesignal,fs):  # add fs
    nsel=thesignal.size
    fsel=fs*(np.arange(0,nsel/2)/nsel)  # add fs*
    ysel=np.fft.fft(thesignal)
    ysel=np.abs(ysel)
    ysel=ysel[0:len(fsel)]
    # ysel=20*np.log(ysel)
    # plt.figure()
    plt.plot(fsel,ysel)
    # plt.show()

# fs=200
# t=np.arange(0,1,1/fs)
# sig=np.cos(2*np.pi*40*t)
# plt.figure()    # function needs
# fftspectrum(sig,fs)
# plt.show() # function needs

def fftspectrum_value(thesignal,fs):  # add fs
    nsel=thesignal.size
    fsel=fs*(np.arange(0,nsel/2)/nsel)  # add fs*
    ysel=np.fft.fft(thesignal)
    ysel=np.abs(ysel)
    ysel=ysel[0:len(fsel)]
    # ysel=20*np.log(ysel)
    # plt.figure()
    # plt.plot(fsel,ysel)
    # plt.show()
    return fsel,ysel

mysubplot.py文件:

import matplotlib.pyplot as plt

def mysubplot(file):
    l=len(file)

    if l==2:
        plt.figure()
        plt.subplot(211);plt.plot(file[0])
        plt.subplot(212);plt.plot(file[1])

    if l==3:
        plt.figure()
        plt.subplot(311);plt.plot(file[0])
        plt.subplot(312);plt.plot(file[1])
        plt.subplot(313);plt.plot(file[2])

    if l==4:
        plt.figure()
        plt.subplot(411);plt.plot(file[0])
        plt.subplot(412);plt.plot(file[1])
        plt.subplot(413);plt.plot(file[2])
        plt.subplot(414);plt.plot(file[3])

    if l==6:
        plt.figure()
        plt.subplot(611);plt.plot(file[0])
        plt.subplot(612);plt.plot(file[1])
        plt.subplot(613);plt.plot(file[2])
        plt.subplot(614);plt.plot(file[3])
        plt.subplot(615);plt.plot(file[4])
        plt.subplot(616);plt.plot(file[5])

    if l==8:
        plt.figure()
        plt.subplot(811);plt.plot(file[0])
        plt.subplot(812);plt.plot(file[1])
        plt.subplot(813);plt.plot(file[2])
        plt.subplot(814);plt.plot(file[3])
        plt.subplot(815);plt.plot(file[4])
        plt.subplot(816);plt.plot(file[5])
        plt.subplot(817);plt.plot(file[6])
        plt.subplot(818);plt.plot(file[7])

    plt.show()

def mysubplot_addx(file):
    l=len(file)

    if l==2:
        plt.figure()
        plt.subplot(211);plt.plot(file[0][0],file[0][1])
        plt.subplot(212);plt.plot(file[1][0],file[1][1])

    if l==3:
        plt.figure()
        plt.subplot(311);plt.plot(file[0][0],file[0][1])
        plt.subplot(312);plt.plot(file[1][0],file[1][1])
        plt.subplot(313);plt.plot(file[2][0],file[2][1])

    if l==4:
        plt.figure()
        plt.subplot(411);plt.plot(file[0][0],file[0][1])
        plt.subplot(412);plt.plot(file[1][0],file[1][1])
        plt.subplot(413);plt.plot(file[2][0],file[2][1])
        plt.subplot(414);plt.plot(file[3][0],file[3][1])

    if l==6:
        plt.figure()
        plt.subplot(611);plt.plot(file[0][0],file[0][1])
        plt.subplot(612);plt.plot(file[1][0],file[1][1])
        plt.subplot(613);plt.plot(file[2][0],file[2][1])
        plt.subplot(614);plt.plot(file[3][0],file[3][1])
        plt.subplot(615);plt.plot(file[4][0],file[4][1])
        plt.subplot(616);plt.plot(file[5][0],file[5][1])

    if l==8:
        plt.figure()
        plt.subplot(811);plt.plot(file[0][0],file[0][1])
        plt.subplot(812);plt.plot(file[1][0],file[1][1])
        plt.subplot(813);plt.plot(file[2][0],file[2][1])
        plt.subplot(814);plt.plot(file[3][0],file[3][1])
        plt.subplot(815);plt.plot(file[4][0],file[4][1])
        plt.subplot(816);plt.plot(file[5][0],file[5][1])
        plt.subplot(817);plt.plot(file[6][0],file[6][1])
        plt.subplot(818);plt.plot(file[7][0],file[7][1])

    plt.show()

展示

在这里插入图片描述
在这里插入图片描述

在这里插入图片描述

发现滤除了50Hz的频率

参考

感谢大佬博主:传送门

猜你喜欢

转载自blog.csdn.net/weixin_41529093/article/details/112831151
hz