2

我正在尝试通过频谱减法来增强超声波信号。信号在时域中并且包含噪声。我将信号划分为 2 µs 的汉明窗,并计算了这些帧的傅里叶变换。然后我选择了 3 个连续的帧,我将其解释为噪声。我对这 3 帧的幅度谱进行了平均,并从每一帧的幅度谱中减去该平均值。然后我将所有负幅值谱定义为零,并通过将新幅值谱与相位谱相结合来重建增强的傅立叶变换。这给了我每帧一系列复数。现在我想通过使用傅里叶逆变换将这个系列转换回时域。但是,此操作为我提供了我不知道如何使用的复数。

我在几篇文章中读到,从傅里叶逆变换获得复数输出是正常的。然而,复数的使用是分开的。有人说要忽略虚部,因为它应该非常小(e-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)

我希望有人知道如何处理这些复数。也许我之前在处理方法上犯了一个错误,但我已经检查了很多次,对我来说似乎很可靠。这是该过程的最后一步(也是最后一步),真正希望在傅里叶逆变换后获得良好的时域信号。

4

2 回答 2

0

使用傅立叶变换转换数据时,一个重要的故障排除辅助工具是您可以执行 fft 然后获取该数据并执行逆 fft 并取回原始数据......我建议您对玩具输入时域数据进行此操作感到满意...让我们说一个 1 Khz 音频波,它是您的时域数据 ...将其发送到 fft 调用中,该调用将返回其频域表示的数组 ...不对该数据进行任何操作,将其发送到逆fft ( ifft ) ...返回的数据将是您原始的 1 Khz 音频波 ...现在就这样做以了解它的功能并在您的项目中使用这个技巧来确认您在球场上...或者如果您从频域数据开始,您也可以这样做......

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

或者

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

在此处查看更多详细信息从 FFT 获取具有最高幅度的频率

于 2019-06-07T13:51:52.807 回答
-1

这是你的问题:

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"]]))
于 2019-06-07T14:03:41.760 回答