Say, I have an example signal consists of three cosines where each of which represents 4, 6 and 8 frequency band. Now, I throw this signal into frequency domain with the use of FFT, and in frequency domain I cut off unwantend 6 Hz band. Finally, I want to inverse the signal from the frequency domain back into time domain. But when I simply use numpy.fft.ifft
I get array of complex numbers, which is not the best result for further analysis of the signal. How can I inverse FFT after performing bandpassing so I'll get whole information carried by real and imaginary part as one number? I looked into z = sqrt(real^2 + imaginary^2)
thing, but it's not the "thing".
Below I provide a working example. I'll be grateful for your help.
import numpy as np
from scipy.fftpack import fftfreq
# Define signal.
Fs = 128 # Sampling rate.
Ts = 1 / Fs # Sampling interval.
Time = np.arange(0, 10, Ts) # Time vector.
signal = np.cos(4*np.pi*Time) + np.cos(6*np.pi*Time) + np.cos(8*np.pi*Time)
def spectrum(sig, t):
"""
Represent given signal in frequency domain.
:param sig: signal.
:param t: time scale.
:return:
"""
f = fftfreq(sig.size, d=t[1]-t[0])
y = np.fft.fft(sig)
return f, np.abs(y)
def bandpass(f, sig, min_freq, max_freq):
"""
Bandpass signal in a specified by min_freq and max_freq frequency range.
:param f: frequency.
:param sig: signal.
:param min_freq: minimum frequency.
:param max_freq: maximum frequency.
:return:
"""
return np.where(np.logical_or(f < min_freq, f > max_freq), 0, sig)
freq, spec = spectrum(signal, Time)
signal_filtered = np.fft.ifft(bandpass(freq, spec, 5, 7))
print(signal_filtered)
"""
print(signal_filtered) result:
[ 2.22833798e-15 +0.00000000e+00j 2.13212081e-15 +6.44480810e-16j
1.85209996e-15 +1.23225456e-15j ..., 1.41336488e-15 -1.71179288e-15j
1.85209996e-15 -1.23225456e-15j 2.13212081e-15 -6.44480810e-16j]
"""
The FFT function returns a result equal to the complex, discrete Fourier transform of Array. The result of this function is a single- or double-precision complex array. The FFT function calls the MKL_FFT function unless it is performing an 8D transform.
We can obtain the magnitude of frequency from a set of complex numbers obtained after performing FFT i.e Fast Fourier Transform in Python. The frequency can be obtained by calculating the magnitude of the complex number. So simple ab(x) on each of those complex numbers should return the frequency.
Those complex numbers in the FFT result are simply just 2 real numbers, which are both required to give you the 2D coordinates of a result vector that has both a length and a direction angle (or magnitude and a phase).
FFT(x, T, Freq) «T» is the time index and «Freq» is the Frequency index, and these two indexes should have the same length. The time series, «x», is indexed by «T» and may be real or complex, and the result is indexed by «Freq» and contains complex numbers in general.
If you want to cut out the frequencies between 5 and 7, then you'll want to keep frequencies where
(f < min_freq) | (f > max_freq)
which is equivalent to
np.logical_or(f < min_freq, f > max_freq)
Therefore, use
return np.where(np.logical_or(f < min_freq, f > max_freq), sig, 0)
instead of
return np.where(np.logical_or(f < min_freq, f > max_freq), 0, sig)
since the second argument to np.where
contains the value returned by np.where
when the condition is True.
With this one change, your code yields
[ 3.00000000 +0.00000000e+00j 2.96514652 +1.24442385e-15j
2.86160515 +2.08976636e-15j ..., 2.69239924 +4.71763845e-15j
2.86160515 +5.88163496e-15j 2.96514652 +6.82134642e-15j]
Note that if your signal is real, you could use rfft
to take the discrete Fourier transform of a real sequence, and irfft
to take its inverse, and rfftfreq
to generate the frequencies.
For example,
from __future__ import division
import numpy as np
import scipy.fftpack as fftpack
# Define signal.
Fs = 128 # Sampling rate.
Ts = 1 / Fs # Sampling interval.
Time = np.arange(0, 10, Ts) # Time vector.
signal = np.cos(4*np.pi*Time) + np.cos(6*np.pi*Time) + np.cos(8*np.pi*Time)
def spectrum(sig, t):
"""
Represent given signal in frequency domain.
:param sig: signal.
:param t: time scale.
:return:
"""
f = fftpack.rfftfreq(sig.size, d=t[1]-t[0])
y = fftpack.rfft(sig)
return f, np.abs(y)
def bandpass(f, sig, min_freq, max_freq):
"""
Bandpass signal in a specified by min_freq and max_freq frequency range.
:param f: frequency.
:param sig: signal.
:param min_freq: minimum frequency.
:param max_freq: maximum frequency.
:return:
"""
return np.where(np.logical_or(f < min_freq, f > max_freq), sig, 0)
freq, spec = spectrum(signal, Time)
signal_filtered = fftpack.irfft(bandpass(freq, spec, 5, 7))
print(signal_filtered)
yields
[ 3. 2.96514652 2.86160515 ..., 2.69239924 2.86160515
2.96514652]
Note you must use scipy
's fftpack
here; do not mix SciPy's implementation with NumPy's.
If you want a strictly real result (minus rounding error noise), the input to your IFFT needs to be hermician symmetric (e.g. you need to make sure the second half of the complex array is the complex conjugate mirror of the first half). Look at your initial FFT of real data and you'll see the symmetry.
But it looks like you didn't filter the negative frequencies, and thus sent an unsymmetrical input to the IFFT which then outputs a complex result.
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With