如何使用Python从50Hz开始到10KHz移除鸣叫信号中的静态噪声

huangapple go评论92阅读模式
英文:

How to remove the stationary noise from the chirp signal start from 50Hz to 10KHz using python

问题

我们有一个从50Hz开始并上升到10KHz的鸣叫信号,系统输出中混入了静态噪音,希望去除与鸣叫混合的静态噪音,保留原始鸣叫信号(数据),欢迎使用Python提出任何建议。

用Python实现对鸣叫信号(50Hz到10KHz)的滤波/静态噪音降低方法。

英文:

we have chirp signal starts from 50Hz and it goes up to 10Hkz with stationary noise mixed in the system output and would like to remove the stationary noise mixed with chirp and to retain original chirp signal(data) and any suggestions would be welcome using python, thanks.

filter/stationary noise reduction approach for the chirp signal (50hz to 10Khz) using python.

答案1

得分: 1

以下是我在您的代码中建议的内容的翻译:

这里是一些代码演示了我在评论中建议的内容

import scipy.signal
import numpy as np
import matplotlib.pyplot as plt

t = np.linspace(0, 5, 44100*5)
w = scipy.signal.chirp(t, f0=50, f1=10000, t1=5, method='linear')
mu, sigma = 0, .1
s = np.random.normal(mu, sigma, len(t))
sig = w + s
n = 0
win_len = 512
hop = 64
win = np.hanning(win_len)
denoised = np.zeros([len(t)])
while n < len(t):
    if n + win_len > len(t):
        # 如果样本不足,就用零填充
        input = np.zeros([win_len])
        end_pt = len(t) - n - 1
        for i in range(0, end_pt):
            input[i] += sig[n + i]
    else:
        input = sig[n:n + win_len]
    input = input * win # 对信号进行窗口处理
    dft = np.fft.fft(input)
    mags = abs(dft)

    ft_idx = 0
    # 将仅包含噪声的频率分量置零
    zero = np.zeros(1, dtype=complex)
    for m in mags:
        if m < 40:
            dft[ft_idx] = zero[0]
        ft_idx += 1
    output = win * np.real(np.fft.ifft(dft)) / (win_len / hop) # 由于输入是实数,所以FFT的逆变换也是实数
    # 我们必须小心地对逆变换也进行窗口处理,以便重叠添加平滑

    # 注意不要溢出
    if n + win_len > len(t):
        end_pt = len(t) - n - 1
        for i in range(0, end_pt):
            denoised[n + i] += output[i]
    else:
        denoised[n:n + win_len] += output
    n += hop
plt.plot(t[0:5000], sig[0:5000])
plt.plot(t[0:5000], denoised[0:5000])
plt.show()

输出:
如何使用Python从50Hz开始到10KHz移除鸣叫信号中的静态噪声

显然,这并不完美。首先,我们保留的信号既受到衰减(由于幅度谱中的粗暴阈值处理)又具有波动峰值幅度(随机噪声的副产品,我们保留了不进行衰减的频率分量中的噪声)。需要更多工作来消除这些效应。


<details>
<summary>英文:</summary>
Here is some code that demonstrates what I suggested in my comment:
import scipy.signal
import numpy as np
import matplotlib.pyplot as plt
t = np.linspace(0,5, 44100*5)
w = scipy.signal.chirp(t, f0=50, f1=10000, t1=5, method=&#39;linear&#39;)
mu, sigma = 0, .1
s = np.random.normal(mu, sigma,len(t))
sig = w+s
n = 0
win_len = 512
hop = 64
win = np.hanning(win_len)
denoised = np.zeros([len(t)])
while n&lt;len(t):
if n+win_len &gt; len(t):
# zero pad if we run out of samples
input = np.zeros([win_len])
end_pt = len(t)-n-1
for i in range(0,end_pt):
input[i]+=sig[n+i]
else:
input = sig[n:n+win_len]
input = input * win # window the signal
dft = np.fft.fft(input)
mags = abs(dft)
ft_idx = 0
# zero out the bins that are only noise
zero = np.zeros(1, dtype=complex)
for m in mags:
if m&lt;40:
dft[ft_idx] = zero[0]
ft_idx +=1 
output = win*np.real(np.fft.ifft(dft))/(win_len/hop)# since the input is real,
# the idft of the fft will also be real
# we must take care to also window the idft so
# that the overlap add is smooth
# take care not to overflow
if n+win_len &gt; len(t):
end_pt = len(t)-n-1
for i in range(0,end_pt):
denoised[n+i]+=output[i]
else:
denoised[n:n+win_len] += output
n+=hop
plt.plot(t[0:5000], sig[0:5000])
plt.plot(t[0:5000], denoised[0:5000])
plt.show()
Output:
[![enter image description here][1]][1]
Obviously it is not perfect. For one thing, the signal we keep is both attenuated (due to the brute force thresholding in the magnitude spedtrum) and the peak amplitudes are warbly (a byproduct of the random noise that we keep in the bins we do not attenuate). More work would be needed to get rid of those effects.
[1]: https://i.stack.imgur.com/1VDAE.png
</details>

huangapple
  • 本文由 发表于 2023年6月29日 04:44:21
  • 转载请务必保留本文链接:https://go.coder-hub.com/76576597.html
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定