總的來(lái)說(shuō)，我想計(jì)算給定數(shù)據(jù)集的傅里葉變換并過(guò)濾掉一些具有最大絕對(duì)值的頻率。所以：1)給定一個(gè)帶有時(shí)間 t 的數(shù)據(jù)數(shù)組 D，2)找到 k 個(gè)最大的傅里葉系數(shù)，3)從數(shù)據(jù)中刪除這些系數(shù)，以便從原始數(shù)據(jù)中濾除某些信號(hào)。在繪制給定時(shí)間內(nèi)過(guò)濾后的數(shù)據(jù)集時(shí)，最終會(huì)出現(xiàn)問(wèn)題。我不太確定錯(cuò)誤在哪里。最終的“過(guò)濾數(shù)據(jù)”圖看起來(lái)甚至沒(méi)有稍微平滑，并且與原始數(shù)據(jù)相比以某種方式改變了其位置。我的代碼完全糟糕嗎？第1部分）：n = 1000limit_low = 0limit_high = 0.48D = 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))**2scaling = (limit_high - limit_low) / (max(D) - min(D))D = D * scalingD = D + (limit_low - min(D)) # given datat = linspace(0,D.size-1,D.size) # times第2部分）：from numpy import linspace import numpy as npfrom scipy import fft, ifftD_f = fft.fft(D) # fft of D-dataset#---extract the k biggest coefficients out of D_f ---k = 15I, 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: breakbigcofs = 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_filterp1=plt.plot(t,D,'r')p2=plt.plot(t,D_new,'b');plt.legend((p1[0], p2[0]), ('original data', 'filtered data'))感謝您的幫助，提前致謝。
1 回答
largeQ
TA貢獻(xiàn)2039條經(jīng)驗(yàn) 獲得超8個(gè)贊
我注意到兩個(gè)問(wèn)題：
您可能想要具有最大絕對(duì)值的分量，因此而
np.argsort(-np.abs(D_f))不是np.argsort(-D_f).更巧妙的是：
bigcofs = np.zeros(len(D_f))是類型float64并且丟棄了行 處的虛部bigcofs[I] = D_f[I]。你可以用以下方法解決這個(gè)問(wèn)題bigcofs = np.zeros(len(D_f), dtype=complex)
我在下面稍微改進(jìn)了您的代碼以獲得所需的結(jié)果：
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()
添加回答
舉報(bào)
0/150
提交
取消