How to use the output of the inverse Fourier transform?

Question

I am trying to enhance an ultrasonic signal by spectral subtraction. The signal is in the time domain and contains noise. I have divided the signal into Hamming windows of 2 µs and calculated the Fourier transforms of those frames. Then I selected 3 consecutive frames which I interpreted as noise. I averaged the magnitude spectra of those 3 frames and subtracted that average from every single frame's magnitude spectrum. Then I defined all negative magnitude spectra as zero and reconstructed the enhanced Fourier transform by combining the new magnitude spectra with the phase spectra. This gives me a series of complex numbers per frame. Now I would like to transform this series back to the time domain by using the inverse Fourier transform. However, this operation provides me with complex numbers which I do not know how to use.

I have read in a couple of posts that it is normal to obtain a complex output from inverse Fourier transformation. However the use of the complex numbers is divided. Some say to neglect the imaginary part, because it is supposed to be very small (e-15) but in my case it is not negligible (0.01-0.5). To be honest I just do not know what to do with the numbers now, because I expected the inverse Fourier transform to give me real number only. And hope for very small imaginary parts, but unfortunately..

# General parameters
#
total_samples = length(time_or) # Total numbers of samples in the current series
max_time = max(time_or) # Length of the measurement in microseconds
sampling_freq = 1/(max_time/1000000)*total_samples # Sampling frequency
frame_length_t = 2 # In microseconds (time)
frame_length_s = round(frame_length_t/1000000*sampling_freq) # In samples per frame
overlap = frame_length_s/2 # Overlap in number of frames, set to 50% overlap
#
# Transform the frame to frequency domain
#
fft_frames = specgram(amp, n=frame_length_s, Fs=125, window=hamming(frame_length_s), overlap=overlap)

mag_spec=abs(fft_frames[["S"]])
phase_spec=atan(Im(fft_frames[["S"]])/Re(fft_frames[["S"]]))

#
# Determine the arrival time of noise
#
cutoff= 10 #determine the percentage of the signal that has to be cut off
dnr=us_data[(length(us_data[,1])*(cutoff/100)):length(us_data[,1]), ]
noise_arr=(length(us_data[,1])-length(dnr[,1])+min(which(dnr[,2]>0.01)))*0.008
#
# Select the frames for noise spectrum estimation
#
noise_spec=0
noise_spec=mag_spec[,noise_arr]
noise_spec=noise_spec+mag_spec[, (noise_arr+1)]
noise_spec=noise_spec+mag_spec[, (noise_arr+2)]
noise_spec_check=noise_spec/3
#
# Subtract the estimated noise spectrum from every frame
#
est_mag_spec=mag_spec-noise_spec_check
est_mag_spec[est_mag_spec < 0] = 0
#
# Transform back to frequency spectrum
#
j=complex(real=0, imaginary=1)
enh_spec = est_mag_spec*exp(j*phase_spec)
#
# Transform back to time domain
#
install.packages("pracma")
library("pracma")
enh_time=fft(enh_spec[,2], inverse=TRUE)

I hope that there is anyone with an idea on how to process these complex numbers. Maybe I have made an mistake earlier on in the processing method, but I have checked it multiple times and it seems quite solid to me. It is the (one but last) last step of the process and really hoping to obtain a nice time domain signal after inverse Fourier transforming.

Answer 1

An essential troubleshooting aid when transforming data using Fourier Transform is the notion that you can do a fft then take that data and do an inverse fft and get back your original data ... I suggest you get comfortable doing this with toy input time domain data ... lets say a 1 Khz audio wave which is your time domain data ... send it into a fft call which will return back an array of its frequency domain representation ... without doing anything with that data send it into an inverse fft ( ifft ) ... the data returned will the your original 1 Khz audio wave ... do that now to gain an appreciation of its power and use this trick on your project to confirm you are in the ball park ... Alternatively if you begin with frequency domain data you can also do this ...

freq domain data -> ifft -> time domain data -> fft -> same freq domain data

or

time domain -> fft -> freq domain -> ifft -> same time domain data

see more details here Get frequency with highest amplitude from FFT

Answer 2

This is your problem:

phase_spec=atan(Im(fft_frames[["S"]])/Re(fft_frames[["S"]]))

Here you compute an angle in a half circle, mapping the other half to the first. That is, you are loosing information.

Many languages have a function to obtain the phase of a complex value, for example in MATLAB it is angle, and in Python numpy.angle.

Alternatively, use the atan2 function, which exists in every single language I’ve ever used, except in NumPy they decided to call it arctan2. It computes the four-quadrant arctangent by taking the two components as separate values. That is, atan(y/x) is the same as atan2(y,x) if the result is in the first two quadrants.

I presume you can do

phase_spec=atan2(Im(fft_frames[["S"]]), Re(fft_frames[["S"]]))