2016-09-30 125 views
3

我想从信号中删除一个频率(一个峰值),并在没有它的情况下绘制我的功能。在fft之后,我发现了频率和振幅,我不确定我现在需要做什么。例如,我想删除我的最高峰(​​在图上标有红点)。如何从信号中删除频率

import numpy as np 
import matplotlib.pyplot as plt 

# create data 
N = 4097 
T = 100.0 
t = np.linspace(-T/2,T/2,N) 
f = np.sin(50.0 * 2.0*np.pi*t) + 0.5*np.sin(80.0 * 2.0*np.pi*t) 

#plot function 
plt.plot(t,f,'r') 
plt.show() 

# perform FT and multiply by dt 
dt = t[1]-t[0] 
ft = np.fft.fft(f) * dt  
freq = np.fft.fftfreq(N, dt) 
freq = freq[:N/2+1] 
amplitude = np.abs(ft[:N/2+1]) 
# plot results 
plt.plot(freq, amplitude,'o-') 
plt.legend(('numpy fft * dt'), loc='upper right') 
plt.xlabel('f') 
plt.ylabel('amplitude') 
#plt.xlim([0, 1.4]) 


plt.plot(freq[np.argmax(amplitude)], max(amplitude), 'ro') 
print "Amplitude: " + str(max(amplitude)) + " Frequency: " + str(freq[np.argmax(amplitude)]) 

plt.show() 
+2

你需要对信号进行滤波。什么样的滤波器以及如何配置它将由您想要保留哪些频率以及您想要删除哪些频率来确定。你可能想要自己介绍一本DSP书,或者从这里开始:https://en.wikipedia.org/wiki/Filter_(signal_processing) – Turn

+0

这篇文章可能很有用:http://forrestbao.blogspot.rs/2014/ 07/signal-filtering-using-inverse-fft-in.html – Prefect

+0

[python中的fft带通滤波器]可能的重复(http://stackoverflow.com/questions/19122157/fft-bandpass-filter-in-python) – strpeter

回答

1

我想是这样的。我知道有更好的解决方案,但我的黑客作品。 :)第二个高峰被删除。

import numpy as np 
import matplotlib.pyplot as plt 
from scipy.fftpack import rfft, irfft, fftfreq, fft 
import scipy.fftpack 

# Number of samplepoints 
N = 500 
# sample spacing 
T = 0.1 

x = np.linspace(0.0, N*T, N) 
y = 5*np.sin(x) + np.cos(2*np.pi*x) 

yf = scipy.fftpack.fft(y) 
xf = np.linspace(0.0, 1.0/(2.0*T), N/2) 
#fft end 

f_signal = rfft(y) 
W = fftfreq(y.size, d=x[1]-x[0]) 

cut_f_signal = f_signal.copy() 
cut_f_signal[(W>0.6)] = 0 

cut_signal = irfft(cut_f_signal) 



# plot results 
f, axarr = plt.subplots(1, 3) 
axarr[0].plot(x, y) 
axarr[0].plot(x,5*np.sin(x),'g') 

axarr[1].plot(xf, 2.0/N * np.abs(yf[:N//2])) 
axarr[1].legend(('numpy fft * dt'), loc='upper right') 
axarr[1].set_xlabel("f") 
axarr[1].set_ylabel("amplitude") 


axarr[2].plot(x,cut_signal) 
axarr[2].plot(x,5*np.sin(x),'g') 

plt.show() 

enter image description here

enter image description here

+1

这很好。从花哨的角度来说,您刚刚在频域中应用了“砖墙低通滤波器”。其他答案(设计一个带阻滤波器)有什么优点和缺点与你做了什么。要看到你所做的工作的缺点,请注意频域中的两个峰都有一些宽度 - 它们不仅仅是两个脉冲。这是因为这两个频率都不在频点上,所以需要所有的频率才能生成。 –

+0

通过切断所有频率> 0.6,结果看起来像一个很好的正弦曲线,但如果您在最右边的曲线上绘制一个*实际*正弦曲线,您将看到通过用“小“的振幅。 –

+0

@AhmedFasih感谢您的注意!你能建议一些书或文章的信号处理? – Prefect

5

你可以设计一个带阻滤波器:

wc = freq[np.argmax(amplitude)]/(0.5/dt) 
wp = [wc * 0.9, wc/0.9] 
ws = [wc * 0.95, wc/0.95] 
b, a = signal.iirdesign(wp, ws, 1, 40) 
f = signal.filtfilt(b, a, f)