Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

How to convert complex numbers back into "normal" numbers after performing FFT

Tags:

python

numpy

fft

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]
"""
like image 747
bluevoxel Avatar asked Apr 09 '15 16:04

bluevoxel


People also ask

What does the FFT return?

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.

How do you extract frequency from FFT?

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.

Why is FFT output complex?

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).

What is the FFT of a complex number?

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.


2 Answers

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.

like image 95
unutbu Avatar answered Oct 19 '22 19:10

unutbu


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.

like image 39
hotpaw2 Avatar answered Oct 19 '22 19:10

hotpaw2