I am used to moving from the time-domain, to the frequency-domain, and then back to the time-domain by:
However I would like to skip step 1 and, instead, define an arbitrary complex (i.e., both real and imaginary component) signal in the frequency-domain and then use an ifft to get the time-domain equivalent.
To ignore the imaginary components (as in step 3 above), several Stack posts suggest the signal needs to be conjugate symmetric (e.g., https://dsp.stackexchange.com/questions/37272/how-to-get-a-real-signal-from-complex-ifft-values, https://dsp.stackexchange.com/questions/74163/what-is-imaginary-part-after-ifft-and-how-to-use-it). Following the definition of a conjugate symmetric signal here:
https://dsp.stackexchange.com/questions/9144/how-to-make-a-signal-conjugate-symmetric
I tried to force my arbitrary frequency-domain signal to be conjugate symmetric. However, when I transform the conjugate symmetric frequency-domain signal to the time-domain via an ifft, I find there are still non-zero values in the imaginary components, which makes me believe I can't simply ignore them. Interestingly, I found that the frequency domain transform of a real valued time-domain signal is not conjugate symmetric but, when transformed back to the time-domain via ifft, has insignificant (close to zero) imaginary components nonetheless.
Am I missing the definition of a conjugate symmetric signal or is there another condition required before I can ignore the imaginary components of an ifft? The code below illustrates my issues.
Finally, I decided to post in StackOverflow instead of SignalProcessing given my goal is to write a python code to define my frequency domain function and transform it to the time-domain.
import numpy as np
from scipy.fft import fft, ifft
# helper function to check if complex conjugate
def check_cc(complex_sig):
sig1 = complex_sig[::-1]
sig2 = np.conj(complex_sig)
sig_check = all(sig1 == sig2)
print(f"Signal Conjugate Symmetric: {sig_check}")
return
""" CASE 1: Damped Sine Wave in Time Domain """
# make a damped sine wave
dt = 0.01
tax = np.arange(1000) * dt + dt
f = 3
wig1 = np.exp(-1 * tax) * np.sin(2 * np.pi * f * tax)
# the frequency domain
WIG1 = fft(wig1)
# check if signal is complex conjugate symmetric
check_cc(WIG1)
# back to the time domain
wig_complex = ifft(WIG1)
wig1_imag = np.imag(wig_complex)
# print the max value of the imaginary time-domain complex signal
print(f"Max Value of wig1 imaginary component: {np.max(wig1_imag)}")
""" Case 2: Arbitraty Complex Signal in the Frequency Domain """
WIG2_real1 = np.arange(len(tax) // 2) * 0
WIG2_real1[1] = 100
WIG2_real2 = WIG2_real1[::-1]
WIG2_real = np.array([*WIG2_real1, *WIG2_real2])
WIG2_imag1 = np.arange(len(tax) // 2) * 0
WIG2_imag1[15] = 100
WIG2_imag2 = -1 * WIG2_imag1[::-1]
WIG2_imag = np.array([*WIG2_imag1, *WIG2_imag2])
WIG2 = np.array([complex(r, i) for r, i in zip(WIG2_real, WIG2_imag)])
# check if singal is complex conjugate
check_cc(WIG2)
# to the time-domain
wig2_complex = ifft(WIG2)
wig2_real = np.real(wig2_complex)
wig2_imag = np.imag(wig2_complex)
# print the max value of the imaginary time-domain complex signal
print(f"Max Value of wig2 imaginary component: {np.max(wig2_imag)}")
We need to establish the axis over which there is symmetry.
In the continuous-domain Fourier transform, where the frequency domain F(ω) has an infinite extent, the symmetry is over the origin. That is, F(ω) = F*(-ω). On the discrete Fourier transform (which is what the FFT computes), the frequency domain F[k] is defined over the integers k = 0 to N-1, for N samples. There is no F[-k]. So how do we make this signal (conjugate) symmetric?
We can take F[k] to be periodic, with F[k] = F[k+mN], for any integer m. So we can imagine a F[-k], which would be equal to F[N-k]. So, it follows that the symmetry is defined by F[k] = F*[N-k].
For this equation to hold at k = 0, we must make F[0] real-valued. Likewise, if N is even, F[N-N/2] = F[N/2], so the element at k = N/2 must also be real-valued (in the odd N case there is no such element).
In Python we can write:
N = 1000 # length of desired signal (even!)
# Create half a
F = np.random.rand(N//2 + 1) + 1j * np.random.rand(N//2 + 1)
F[0] = F[0].real
F[-1] = F[-1].real
F = np.concatenate((F, np.flip(np.conj(F[1:-1]))))
f = np.fft.ifft(F)
assert np.real_if_close(f).dtype == np.float64
To test for conjugate symmetry, we'd have to build the reverse array this way:
F_flipped = np.concatenate((F[:1], np.flip(F[1:])))
np.all(np.isclose(F, np.conj(F_flipped)))
(Note that your check all(sig1 == sig2) will fail in general, you cannot compare floating-point numbers for equality, you need to include a tolerance, which is what np.isclose() does.)