Scipy巴特沃思带通滤波器未返回正确的振幅

hiz5n14c 于 2023-08-05 发布在其他

关注(0)|答案(1)|浏览(94)

我试图使用SciPy巴特沃思创建一个带通滤波器，但它没有返回正确的振幅（尽管它似乎得到了正确的周期）。首先，我创建一个由三个波组成的波：

x = np.arange(1,365,1)
y1 = 1*np.sin( (1/10) * x)
y2 = 5*np.sin( (1/45) * x)
y3 = 10*np.sin( (1/120) * x)

y_total = y1 + y2 + y3

字符串
然后我使用过滤器：

import math
from scipy.signal import butter, lfilter

def firstNonNan(listfloats):
    for item in listfloats:
        if math.isnan(item) == False:
            #print(item)
            return item

def butter_bandpass(highcut, lowcut, fs, order=2):
    nyq = 0.5 * fs
    high = highcut / nyq

    low = lowcut / nyq

    b, a = butter(order, [high, low], btype='band')
    return b, a

def butter_bandpass_filter(data, cap, cap1, order=2):
    highcut = 1/cap

    lowcut = 1/cap1

    fs = 1
    b, a = butter_bandpass(highcut,lowcut, fs, order=order)

    y = lfilter(b, a, data)
    return y

def butter_hpf(data, cap, cap1):
    indc = firstNonNan(list(data))
    year1 = list(np.random.uniform(low=indc-0.1, high=indc+0.1, size=(cap,))) + list(data)
    return(butter_bandpass_filter(year1, cap, cap1)[cap:])

型
（我在两端使用随机生成的数据来尝试和处理边缘效应）
过滤器的应用如下：

plt.scatter(x,y_total)

plt.scatter(x,butter_hpf(y_total,96,30))

型
其中butter_hpf函数采用时间序列以及高和低截止周期。我希望它只返回一个振幅为5的正弦波，但正如你在下面看到的，振幅太小，波不完全只是正弦波（蓝色是输入波，橙子是滤波器输出）。如何修复我的过滤器？

的数据

scipy

来源：https://stackoverflow.com/questions/76639767/scipy-butterworth-bandpass-filter-does-not-return-correct-amplitudes

1条答案

按热度按时间

6tdlim6h1#

np.sin和大多数（可能是所有）其他Scipy基三角函数以弧度工作;你在“循环”中工作;你的频率是（我认为）“每天”，但你需要将它们缩放到“每天弧度”，以生成y1等。
结果仍然不是很好，尽管它可能已经足够好了。使用filtfilt会有一个小的改进。

import matplotlib.pyplot as plt
import numpy as np
import math
from scipy.signal import butter, lfilter, filtfilt, iirdesign, sosfilt

def firstNonNan(listfloats):
    for item in listfloats:
        if math.isnan(item) == False:
            #print(item)
            return item

def butter_bandpass(highcut, lowcut, fs, order=2):
    nyq = 0.5 * fs
    high = highcut / nyq

    low = lowcut / nyq

    b, a = butter(order, [high, low], btype='band')
    return b, a

def butter_bandpass_filter(data, cap, cap1, order=2):
    highcut = 1/cap

    lowcut = 1/cap1

    fs = 1
    b, a = butter_bandpass(highcut,lowcut, fs, order=order)

    y = lfilter(b, a, data)
    return y

def butter_hpf(data, cap, cap1):
    indc = firstNonNan(list(data))
    year1 = list(np.random.uniform(low=indc-0.1, high=indc+0.1, size=(cap,))) + list(data)
    return(butter_bandpass_filter(year1, cap, cap1)[cap:])

x = np.arange(1,365,1)
y1 = 1*np.sin( 2*np.pi*(1/10) * x)
y2 = 5*np.sin( 2*np.pi*(1/45) * x)
y3 = 10*np.sin( 2*np.pi*(1/120) * x)

y_total = y1 + y2 + y3

fig, ax = plt.subplots()
ax.plot(x,y_total)
ax.plot(x,butter_hpf(y_total,96,30))
ax.plot(x,y2)
ax.legend(['noisy','filtered','desired'])
ax.set_title('using rad/day')

# filtfilt
b, a = butter(2, [1/96, 1/30], btype='band', fs=1)
yff = filtfilt(b, a, y_total, method='gust')

fig, ax = plt.subplots()
ax.plot(x,y_total)
ax.plot(x,yff)
ax.plot(x,y2)
ax.legend(['noisy','filtered','desired'])
ax.set_title('using rad/day and filtfilt')

字符串
x1c 0d1x的数据

的

赞(0）回复(0）举报 2023-08-05

我来回答

Scipy巴特沃思带通滤波器未返回正确的振幅

1条答案

相关问题

热门标签

最新问答