FFT：在信號中找到並切割50Hz的噪音

Question

我有噪音data（峯值週期爲1.8s，每個週期爲2048個bin），我想計算頻率並刪除50Hz。我很確定，我期待的頻率是50Hz，因爲我通過使用originlab找到它。

當我嘗試在python中做同樣的操作時，平均峯值是〜47Hz。我正在看教程和例子，但結果仍然是一樣的。

import numpy as np 
from scipy.fftpack import fft 
from scipy.fftpack import fftfreq 
import matplotlib.pyplot as plt 
data = np.loadtxt('3.dat', comments="#") 
t = data[:, 0] 
y = data[:, 2] 
len_data = len(data) 
bins = 2048 
plt.figure(figsize=(7, 9)) 
plt.subplot(211) 
plt.plot(t, y, 'b-') 
plt.xlabel("time[sec]") 
plt.ylabel("original signal") 
plt.subplot(212) 
F = fft(y) 
freq = fftfreq(len(t), (t[1] - t[0])) 
ipos = np.where(freq > 0) 
freqs = freq[ipos] 
mags = np.abs(F[ipos]) 
plt.plot(freqs, mags, 'b-') 
plt.xlabel("freq") 
plt.ylabel("POWER") 
plt.savefig('stoc.png') 
plt.show() 

有人可以幫助我如何解決？

我必須恢復關於隔絕噪音的問題。當我減去頻率時，信號幅度顯着下降。它是否正確？

data = np.loadtxt('3.dat', comments="#") 
t = data[:, 0] 
phase = data[:, 1] 
y = data[:, 2] 
pulse_no = data[:, 3] 
len_data = len(data) 
bins = 2048 
ti = np.linspace(t[0], t[-1], len_data) 
yi = np.interp(ti, t, y) 
t, y = ti, yi 

plt.figure(figsize=(10, 10)) 
plt.subplot(511) 
plt.plot(t, y, 'b-') 
plt.xlabel("time[sec]") 
plt.ylabel("original signal") 
plt.subplot(512) 
F = fft(y) 
N = len(t) 
w = fftfreq(N, (t[1] - t[0])) 
ipos = np.where(w > 0) 
freq = w[ipos] 
mags = abs(F[ipos]) 
plt.plot(freq, mags) 
ip = np.where(F > 0)[0] 
Fs = np.copy(F) 
yf = ifft(Fs) 
ip = np.where(F > 0)[0] 
Ff = np.copy(F) 
Ff[ip > ip[[(181)]]] = 0 
Ff[ip < ip[[(175)]]] = 0 
magsf = abs(Ff[ipos]) 
plt.plot(freq, magsf, 'r-') 
plt.subplot(513) 
Fr = mags - magsf 
plt.plot(freq, Fr) 
plt.subplot(514) 
yf = ifft(Ff) 
yr = ifft(Fr) 
plt.plot(t, yf) 
plt.subplot(515) 
flux = y - np.real(yf) 
plt.plot(t, flux) 
plt.plot(t, y) 
plt.show() 

Answer 1

您的問題似乎是，你的時間網格並不均勻間隔：

In [83]: d = np.diff(data[:,0]) 

In [84]: d 
Out[84]: 
array([ 0.0006144 , 0.0006144 , 0.00049152, ..., 0.0006144 , 
     0.0006144 , 0.00049152]) 

如果我插的值，以恆定的間隔時間：

data = np.loadtxt('3.dat', comments="#") 
t = data[:, 0] 
y = data[:, 2] 
len_data = len(data) 

ti = np.linspace(t[0], t[-1], len_data) 
yi = np.interp(ti, t, y) 
t, y = ti, yi 

峯值爲50Hz：