傅里叶变换和过滤给定的数据集

Question

总的来说，我想计算给定数据集的傅立叶变换，并过滤掉一些绝对值最大的频率。所以：

1)给定一个包含时间t的数据数组D，2)找到k个最大的傅里叶系数和3)从数据中去除这些系数，以便从原始数据中过滤掉某些信号。

在给定时间绘制过滤后的数据集时，最后出了点问题。我不确定错误在哪里。与原始数据相比，最终的 'filtered data' 图看起来甚至没有稍微平滑，并且以某种方式改变了它的位置。我的代码完全糟糕吗？

第 1 部分):

n = 1000
limit_low = 0
limit_high = 0.48
D = np.random.normal(0, 0.5, n) + np.abs(np.random.normal(0, 2, n) * np.sin(np.linspace(0, 3*np.pi, n))) + np.sin(np.linspace(0, 5*np.pi, n))**2 + np.sin(np.linspace(1, 6*np.pi, n))**2

scaling = (limit_high - limit_low) / (max(D) - min(D))
D = D * scaling
D = D + (limit_low - min(D))                       # given data

t = linspace(0,D.size-1,D.size)                    # times

第 2 部分):

from numpy import linspace                                    
import numpy as np
from scipy import fft, ifft

D_f = fft.fft(D)         # fft of D-dataset

#---extract the k biggest coefficients out of D_f ---

k = 15
I, bigvals = [], []                                        
for i in np.argsort(-D_f):                
    if D_f[i] not in bigvals:            
        bigvals.append(D_f[i])          
        I.append(i)
        if len(I) == k:
            break

bigcofs = np.zeros(len(D_f))             
bigcofs[I] = D_f[I]                      # array with only zeros in in except for the k maximal coefficients

第 3 部分):

D_filter = fft.ifft(bigcofs)
D_new = D - D_filter

p1=plt.plot(t,D,'r')
p2=plt.plot(t,D_new,'b');
plt.legend((p1[0], p2[0]), ('original data', 'filtered data'))

感谢您的帮助，在此先致谢。

Answer 1

我注意到两个问题：

您可能想要绝对值最大的分量，因此 np.argsort(-np.abs(D_f)) 而不是 np.argsort(-D_f)。
更巧妙：bigcofs = np.zeros(len(D_f)) 属于 float64 类型，并在行 bigcofs[I] = D_f[I] 处丢弃了虚部。您可以使用 bigcofs = np.zeros(len(D_f), dtype=complex)
解决此问题

我在下面稍微改进了您的代码以获得期望的结果：

import numpy as np
from scipy import fft, ifft
import matplotlib.pyplot as plt

n = 1000
limit_low = 0
limit_high = 0.48
N_THRESH = 10

D = 0.5*np.random.normal(0, 0.5, n) + 0.5*np.abs(np.random.normal(0, 2, n) * np.sin(np.linspace(0, 3*np.pi, n))) + np.sin(np.linspace(0, 5*np.pi, n))**2 + np.sin(np.linspace(1, 6*np.pi, n))**2

scaling = (limit_high - limit_low) / (max(D) - min(D))
D = D * scaling
D = D + (limit_low - min(D))                       # given data
t = np.linspace(0,D.size-1,D.size)                    # times

# transformed data
D_fft = fft.fft(D)

# Create boolean mask for N largest indices
idx_sorted = np.argsort(-np.abs(D_fft))
idx = idx_sorted[0:N_THRESH]
mask = np.zeros(D_fft.shape).astype(bool)
mask[idx] = True

# Split fft above, below N_THRESH points:
D_below = D_fft.copy()
D_below[mask] = 0
D_above = D_fft.copy() 
D_above[~mask] = 0

#inverse separated functions
D_above = fft.ifft(D_above)
D_below = fft.ifft(D_below)

# plot
plt.ion()
f, (ax1, ax2, ax3) = plt.subplots(3,1)
l1, = ax1.plot(t, D, c="r", label="original")
l2, = ax2.plot(t, D_above, c="g", label="top {} coeff. signal".format(N_THRESH))
l3, = ax3.plot(t, D_below, c="b", label="remaining signal")
f.legend(handles=[l1,l2,l3])
plt.show()

傅里叶变换和过滤给定的数据集

Fourier transform and filter given data set

python

fft

filter