在频域中过滤信号振幅最高的两个频率

Question

我试图过滤振幅最高的两个频率。我想知道结果是否正确，因为过滤后的信号似乎不如原始信号平滑？ FFT 函数的输出包含基频 A0/C0 是否正确，将其包含在最大幅度的搜索中是否正确（它确实是最高的！）？

我的代码（基于我的教授和同事的代码，到目前为止我还没有理解每个细节）：

# signal
data = np.loadtxt("profil.txt")
t = data[:,0]
x = data[:,1]
x = x-np.mean(x)           # Reduce signal to mean
n = len(t)
max_ind = int(n/2-1)       
dt = (t[n-1]-t[0])/(n-1)   
T =    n*dt                  
df = 1./T                  

# Fast-Fourier-Transformation
c     = 2.*np.absolute(fft(x))/n    #get the power sprectrum c from the array of complex numbers
c[0]  = c[0]/2.                     #correction for c0 (fundamental frequency)
f = np.fft.fftfreq(n, d=dt)       
a = fft(x).real
b = fft(x).imag
n_fft = len(a)

# filter 
p = np.ones(len(c))
p[c[0:int(len(c)/2)].argsort()[int(len(c)/2-1)]] = 0         #setting the positions of p to 0 with 
p[c[0:int(len(c)/2)].argsort()[int(len(c)/2-2)]] = 0         #the indices from the argsort function
print(c[0:int(len(c)/2-1)].argsort()[int(n_fft/2-2)])        #over the first half of the c array,  
ab_filter_2 = fft(x)                                        #because the second half contains the 
ab_filter_2.real = a*p                                      #negative frequencies.  
ab_filter_2.imag = b*p
x_filter2 = ifft(ab_filter_2)*2

我不太了解 fft 返回负频率和正频率的全部内容。我知道它们只是镜像，但为什么我不能搜索整个阵列？而 ifft 函数仅适用于正频率数组？

结果图：（蓝色原件，红色被过滤）： enter image description here

Answer 1

这部分非常浪费：

a = fft(x).real
b = fft(x).imag

您无缘无故地计算了两次 FFT。您稍后第三次计算它，并且您之前已经计算过一次。您应该只计算一次，而不是 4 次。 FFT 是代码中最昂贵的部分。

然后：

ab_filter_2 = fft(x)                                        #because the second half contains the 
ab_filter_2.real = a*p                                      #negative frequencies.  
ab_filter_2.imag = b*p
x_filter2 = ifft(ab_filter_2)*2

将所有这些替换为：

out = ifft(fft(x) * p)

在这里你做了两次同样的事情：

p[c[0:int(len(c)/2)].argsort()[int(len(c)/2-1)]] = 0
p[c[0:int(len(c)/2)].argsort()[int(len(c)/2-2)]] = 0

但是你只设置了过滤器的左半部分。制作对称滤波器很重要。 abs(f) 有两个位置具有相同的值（直到舍入误差！），正频率和负频率一起出现。这两个位置应该具有相同的滤波器值（实际上是复共轭，但您有一个实值滤波器，因此在这种情况下差异无关紧要）。

我不确定索引到底做了什么。为了便于阅读，我会将声明拆分成单独行中的较短部分。

FFT 函数的输出包含基频 A0/C0 [...] 是否正确？

原则上是的，但你从信号中减去平均值，有效地将基频（直流分量）设置为 0。