4

I want to remove one frequency (one peak) from signal and plot my function without it. After fft I found frequency and amplitude and I am not sure what I need to do now. For example I want to remove my highest peak (marked with red dot on plot).

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()
Prefect
  • 153
  • 1
  • 1
  • 7
  • 2
    You need to filter the signal. What kind of filter and how you configure it is going to be determined by both which frequencies you want to keep and which you want to remove. You might want to get yourself an intro to DSP book or start here: https://en.wikipedia.org/wiki/Filter_(signal_processing) – Turn Sep 30 '16 at 22:41
  • This post may be useful: http://forrestbao.blogspot.rs/2014/07/signal-filtering-using-inverse-fft-in.html – Prefect Oct 01 '16 at 01:31
  • Possible duplicate of [fft bandpass filter in python](http://stackoverflow.com/questions/19122157/fft-bandpass-filter-in-python) – strpeter May 09 '17 at 19:50

2 Answers2

8

You can design a bandstop filter:

from scipy import signal
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)
MD004
  • 581
  • 1
  • 7
  • 19
HYRY
  • 94,853
  • 25
  • 187
  • 187
8

One option is to transform the signal to the frequency domain then remove the selected frequency.

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

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

x = np.linspace(0.0, (N-1)*T, N)
# x = np.arange(0.0, N*T, T)  # alternate way to define x
y = 5*np.sin(x) + np.cos(2*np.pi*x) 

yf = 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  # filter all frequencies above 0.6

cut_signal = irfft(cut_f_signal)

# plot results
f, axarr = plt.subplots(1, 3, figsize=(9, 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

Steven C. Howell
  • 16,902
  • 15
  • 72
  • 97
Prefect
  • 153
  • 1
  • 1
  • 7
  • 3
    This is totally fine. In fancy terms, you’ve just applied a “brick-wall low-pass filter in the frequency-domain”. What the other answer (designing a bandstop filter) has pros and cons versus what you did. To see the drawback of what you did, notice how both peaks in the frequency domain have some width to them—they’re not just two impulses. That’s because neither frequency lies exactly on a frequency bin, so it takes *all* frequencies to make them. – Ahmed Fasih Oct 01 '16 at 02:09
  • 1
    By cutting off all frequencies >0.6, the result looks like a nice sinusoid, but if you plot an *actual* sinusoid on top of your right-most plot, you’ll see the subtle distortions introduced by filtering out the frequencies with “small” amplitudes. – Ahmed Fasih Oct 01 '16 at 02:10
  • @AhmedFasih thanks for notices! Can you suggest some book or article for signal processing? – Prefect Oct 01 '16 at 13:13
  • I think http://www.dspguide.com/ might be helpful to start even for those without math background (i.e., if you’re not comfortable with integrals and complex exponentials). I haven’t read this but http://dsp-nbsphinx.readthedocs.io/en/nbsphinx-experiment/ teaches a DSP course using IPython Notebooks, so hopefully it’s very interactive and hands-on, and has lots of demos. Happy learning—please write a book based on what you learn! (And ask more FFT/DSP questions. There’s a twitter bot that tweets FFT questions: https://twitter.com/StackOverFFT !) – Ahmed Fasih Oct 01 '16 at 15:02
  • Instead of zeroing out, `use np.min(, cut_f_signal)` after the first peak to retain the small contributions. – percusse Oct 01 '16 at 23:56