2016-09-30 123 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)