Scipy 非常嘈杂信号的 FFT 频率分析
Scipy FFT Frequency Analysis of very noisy signal
我有噪声 data,我想为其计算频率和振幅。每 1/100 秒收集一次样本。从趋势来看,我认为频率约为 0.3
当我使用 numpy
fft
模块时,我最终得到非常高的频率(36.32 /秒),这显然是不正确的。我尝试用 pandas
rolling_mean
过滤数据以去除 fft 之前的噪声,但这也不起作用。
import pandas as pd
from numpy import fft
import numpy as np
import matplotlib.pyplot as plt
Moisture_mean_x = pd.read_excel("signal.xlsx", header = None)
Moisture_mean_x = pd.rolling_mean(Moisture_mean_x, 10) # doesn't helps
Moisture_mean_x = Moisture_mean_x.dropna()
Moisture_mean_x = Moisture_mean_x -Moisture_mean_x.mean()
frate = 100. #/sec
Hn = fft.fft(Moisture_mean_x)
freqs = fft.fftfreq(len(Hn), 1/frate)
idx = np.argmax(np.abs(Hn))
freq_in_hertz = freqs[idx]
有人可以指导我如何解决这个问题吗?
希望对您有所帮助。
https://scipy-cookbook.readthedocs.io/items/ButterworthBandpass.html
在应用 FFT 之前,您应该只过滤预期频率附近的频带并提高信噪比。
编辑:
Mark Ransom给出了一个更聪明的答案,如果你必须做FFT,你可以在转换后切断噪音。它不会产生比过滤器更差的结果。
你应该使用低通滤波器,它应该保留较大的周期性变化并首先平滑掉一些较高频率的东西。之后,就可以做 FFT 得到峰值。这是一个 recipe for FIR filter 通常用于这类事情。
FFT 是一个滤波器组。只需在 FFT 结果(而不是整个结果向量)中寻找预期频率范围内的幅度峰值 仅,其他大部分频谱将基本上被滤除。
不需要预先对信号进行滤波,因为FFT是一个滤波器。只需跳过 FFT 中与您知道包含大量噪声的频率相对应的那些部分 - 将它们归零,或者以其他方式排除它们。
你是对的,有问题。需要明确询问 pandas 第零列:
Hn = np.fft.fft(Moisture_mean_x[0])
否则会发生错误,您可以从 FFT 结果不对称这一事实看出这一点,实际输入应该是这种情况。
似乎 @tillsten 已经回答了你的问题,但这里有一些额外的确认。第一个图是您的数据(零均值,我将其更改为 csv)。第二个是功率谱密度,您可以看到峰值在 ~0.3 Hz 处的脂肪量。我'zoomed'在第三个图上看是否有第二个隐藏频率靠近主频率
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
x = pd.read_csv("signal.csv")
x = np.array(x, dtype=float)[:,0]
x = x - np.mean(x)
fs = 1e2
f, Pxx = signal.welch(x, fs, nperseg=1024)
f_res, Pxx_res = signal.welch(x, fs, nperseg=2048)
plt.subplot(3,1,1)
plt.plot(x)
plt.subplot(3,1,2)
plt.plot(f, Pxx)
plt.xlim([0, 1])
plt.xlabel('frequency [Hz]')
plt.ylabel('PSD')
plt.subplot(3,1,3)
plt.plot(f_res, Pxx_res)
plt.xlim([0, 1])
plt.xlabel('frequency [Hz]')
plt.ylabel('PSD')
plt.show()
Hn = fft.fft(x)
freqs = fft.fftfreq(len(Hn), 1/fs)
idx = np.argmax(np.abs(Hn))
freq_in_hertz = freqs[idx]
print 'Main freq:', freq_in_hertz
print 'RMS amp:', np.sqrt(Pxx.max())
这会打印:
Main freq: 0.32012805122
RMS amp: 0.0556044913489
我有噪声 data,我想为其计算频率和振幅。每 1/100 秒收集一次样本。从趋势来看,我认为频率约为 0.3
当我使用 numpy
fft
模块时,我最终得到非常高的频率(36.32 /秒),这显然是不正确的。我尝试用 pandas
rolling_mean
过滤数据以去除 fft 之前的噪声,但这也不起作用。
import pandas as pd
from numpy import fft
import numpy as np
import matplotlib.pyplot as plt
Moisture_mean_x = pd.read_excel("signal.xlsx", header = None)
Moisture_mean_x = pd.rolling_mean(Moisture_mean_x, 10) # doesn't helps
Moisture_mean_x = Moisture_mean_x.dropna()
Moisture_mean_x = Moisture_mean_x -Moisture_mean_x.mean()
frate = 100. #/sec
Hn = fft.fft(Moisture_mean_x)
freqs = fft.fftfreq(len(Hn), 1/frate)
idx = np.argmax(np.abs(Hn))
freq_in_hertz = freqs[idx]
有人可以指导我如何解决这个问题吗?
希望对您有所帮助。
https://scipy-cookbook.readthedocs.io/items/ButterworthBandpass.html
在应用 FFT 之前,您应该只过滤预期频率附近的频带并提高信噪比。
编辑:
Mark Ransom给出了一个更聪明的答案,如果你必须做FFT,你可以在转换后切断噪音。它不会产生比过滤器更差的结果。
你应该使用低通滤波器,它应该保留较大的周期性变化并首先平滑掉一些较高频率的东西。之后,就可以做 FFT 得到峰值。这是一个 recipe for FIR filter 通常用于这类事情。
FFT 是一个滤波器组。只需在 FFT 结果(而不是整个结果向量)中寻找预期频率范围内的幅度峰值 仅,其他大部分频谱将基本上被滤除。
不需要预先对信号进行滤波,因为FFT是一个滤波器。只需跳过 FFT 中与您知道包含大量噪声的频率相对应的那些部分 - 将它们归零,或者以其他方式排除它们。
你是对的,有问题。需要明确询问 pandas 第零列:
Hn = np.fft.fft(Moisture_mean_x[0])
否则会发生错误,您可以从 FFT 结果不对称这一事实看出这一点,实际输入应该是这种情况。
似乎 @tillsten 已经回答了你的问题,但这里有一些额外的确认。第一个图是您的数据(零均值,我将其更改为 csv)。第二个是功率谱密度,您可以看到峰值在 ~0.3 Hz 处的脂肪量。我'zoomed'在第三个图上看是否有第二个隐藏频率靠近主频率
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
x = pd.read_csv("signal.csv")
x = np.array(x, dtype=float)[:,0]
x = x - np.mean(x)
fs = 1e2
f, Pxx = signal.welch(x, fs, nperseg=1024)
f_res, Pxx_res = signal.welch(x, fs, nperseg=2048)
plt.subplot(3,1,1)
plt.plot(x)
plt.subplot(3,1,2)
plt.plot(f, Pxx)
plt.xlim([0, 1])
plt.xlabel('frequency [Hz]')
plt.ylabel('PSD')
plt.subplot(3,1,3)
plt.plot(f_res, Pxx_res)
plt.xlim([0, 1])
plt.xlabel('frequency [Hz]')
plt.ylabel('PSD')
plt.show()
Hn = fft.fft(x)
freqs = fft.fftfreq(len(Hn), 1/fs)
idx = np.argmax(np.abs(Hn))
freq_in_hertz = freqs[idx]
print 'Main freq:', freq_in_hertz
print 'RMS amp:', np.sqrt(Pxx.max())
这会打印:
Main freq: 0.32012805122
RMS amp: 0.0556044913489