Как использовать выходные данные обратного преобразования Фурье?

Question

Я пытаюсь усилить ультразвуковой сигнал путем спектрального вычитания. Сигнал находится во временной области и содержит шум. Я разделил сигнал на окна Хэмминга длительностью 2 мкс и рассчитал преобразования Фурье этих кадров. Затем я выбрал 3 последовательных кадра, которые я интерпретировал как шум. Я усреднил спектры величин этих трех кадров и вычел это среднее из спектра величин каждого отдельного кадра. Затем я определил все спектры отрицательной величины как ноль и восстановил усиленное преобразование Фурье, объединив новые спектры амплитуды с фазовыми. Это дает мне серию комплексных чисел на кадр. Теперь я хотел бы преобразовать этот ряд обратно во временную область, используя обратное преобразование Фурье. Однако эта операция предоставляет мне комплексные числа, которые я не знаю, как использовать.

Я прочитал в нескольких сообщениях, что это нормально, чтобы получить комплексный результат обратного преобразования Фурье. Однако использование комплексных чисел разделено. Некоторые говорят, что пренебрегать мнимой частью, потому что она должна быть очень маленькой (15), но в моем случае она не пренебрежимо мала (0,01-0,5). Честно говоря, я просто не знаю, что теперь делать с числами, потому что ожидал, что обратное преобразование Фурье даст мне только действительное число. И надеюсь на очень маленькие мнимые части, но, к сожалению..

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

Я надеюсь, что есть кто-нибудь с идеей о том, как обрабатывать эти сложные числа. Возможно, я ранее допустил ошибку в методе обработки, но я проверял ее несколько раз, и она кажется мне довольно серьезной. Это (последний, но последний) последний шаг процесса, и он действительно надеется получить хороший сигнал во временной области после обратного преобразования Фурье.

Answer 1

Важным средством устранения неполадок при преобразовании данных с использованием преобразования Фурье является представление о том, что вы можете выполнить fft, затем взять эти данные и выполнить обратное fft и вернуть ваши исходные данные... Я предлагаю вам освоиться с игрушечными входными данными во временной области... допустим, звуковая волна с частотой 1 кГц, которая является данными вашей временной области... отправить ее в вызов fft, который вернет обратно массив его представления в частотной области... ничего не делая с этими данными, отправьте их в обратную fft ( ifft) ... возвращенные данные будут вашей исходной звуковой волной 1 кГц... сделайте это сейчас, чтобы оценить ее мощь и используйте этот трюк в своем проекте, чтобы подтвердить, что вы находитесь в парке мячей... В качестве альтернативы если вы начнете с данными в частотной области, вы также можете сделать это...

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

или же

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

подробности смотрите здесь. Получите частоту с наибольшей амплитудой из FFT

Answer 2

Это ваша проблема:

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

Здесь вы вычисляете угол в полукруге, сопоставляя другую половину с первой. То есть вы теряете информацию.

Многие языки имеют функцию для получения фазы комплексного значения, например, в MATLAB это angle и в Python numpy.angle,

В качестве альтернативы используйте atan2 функция, которая существует на каждом языке, который я когда-либо использовал, за исключением NumPy, который они решили назвать arctan2, Он вычисляет арктангенс четырех квадрантов, принимая два компонента как отдельные значения. Это, atan(y/x) такой же как atan2(y,x) если результат в первых двух квадрантах.

Я полагаю, вы можете сделать

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