首页 > 解决方案 > Python - 执行 FFT 忽略来自 MEMS 麦克风的直流偏移

问题描述

我正在尝试对 wav 文件执行 FFT,效果很好,但在我的图中,我看到 0 Hz 处的振幅很大。我认为这是一个直流偏移。我的目的是在我的绘图中或直接在代码中忽略这个 DC 偏移,因为它阻止我看到实际的噪音。在我的示例中,我记录了大约 6.1kHz 的噪声,我可以清楚地看到它,如果我在该点放大,但在一般情况下,由于 0Hz 的幅度,它是不可观察的。如果您告诉我如何忽略 0Hz(或 DC 偏置),我会很高兴。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from __future__ import print_function
import scipy.io.wavfile as wavfile
import scipy
import scipy.fftpack
import numpy as np
from matplotlib import pyplot as plt

fs_rate, signal = wavfile.read("file.wav")
print ("Frequency sampling", fs_rate)
l_audio = len(signal.shape)
print ("Channels", l_audio)
if l_audio == 2:
    signal = signal.sum(axis=1) / 2
N = signal.shape[0]
print ("Complete Samplings N", N)
secs = N / float(fs_rate)
print ("secs", secs)
Ts = 1.0/fs_rate # sampling interval in time
print ("Timestep between samples Ts", Ts)
t = scipy.arange(0, secs, Ts) # time vector as scipy arange field / numpy.ndarray
FFT = abs(scipy.fft(signal))
FFT_side = FFT[range(N/4)] # one side FFT range
freqs = scipy.fftpack.fftfreq(signal.size, t[1]-t[0])
fft_freqs = np.array(freqs)
freqs_side = freqs[range(N/4)] # one side frequency range
fft_freqs_side = np.array(freqs_side)

print (abs(FFT_side))

plt.subplot(211)
p1 = plt.plot(t, signal, "g") # plotting the signal
plt.xlabel('Time')
plt.ylabel('Amplitude')

plt.subplot(212)
p2 = plt.plot(freqs_side, abs(FFT_side), "b") # plotting the positive fft spectrum
plt.xlabel('Frequency (Hz)')
plt.ylabel('Count single-sided')
plt.show()

放大图 正常情节

标签: pythonmatplotlibscipyfftfrequency

解决方案


较大的偏移量通常表示信号预处理不正确。常见的方法包括贬低数据和使用线性回归对数据进行漂移。这是一个例子

from matplotlib.pyplot import *

from numpy import *

dt = 1/1000
T  = 1
t  = arange(0, T, dt)

n = t.size
y = sin(pi * t * 3)   + 39  + 3 * t + random.rand(n)

from scipy import optimize
# subtract drift
lin = lambda x, a, b : a * x + b 
coeff, _ = optimize.curve_fit(lin,t, y)
dmy= y- coeff[0] * t + coeff[0]

# compute power
fy = abs(fft.fft(y))[:n//2] ** 2
fyn= abs(fft.fft(dmy - dmy.mean()))[:n//2] ** 2 # NB demeaned

freq= linspace(0, T / dt, n//2) # get freqs
fig, ax = subplots(2, sharex = 'all')
for axi, data, label in zip(ax, [fy,fyn], 'raw processed'.split()):
    axi.plot(freq, data)
    axi.set(xlim = (0, 10), title = label)
axi.set_xlabel('freq')
subplots_adjust(hspace = .5)

在此处输入图像描述


推荐阅读