我正在执行以下操作：执行 FFT / 在 FFT 结果中剪切每个高于 100Hz 的频率 / 执行逆 FFT
它运作良好......如果原始数据集没有偏移量！如果它有偏移，则输出结果幅度被破坏。
例子：
从数学上讲，我什至不确定我能做我正在做的事情。我所能观察到的是，有了偏移，基频是原来的两倍？？？！！！
你知道为什么偏移量会改变吗？
编码：
def FFT(data,time_step):
"""
Perform FFT on raw data and return result of FFT (yf) and frequency axis (xf).
"""
"""
#Test code for the manual frequency magnitude plot
import numpy as np
import matplotlib.pyplot as plt
#Generate sinus waves
x = np.linspace(0,2*np.pi,50000) #You need enough points to be able to capture information (Shannon theorem)
data = np.sin(x*2*np.pi*50) + 0.5*np.sin(x*2*np.pi*200)
time_step = (x[-1]-x[0])/x.size
list_data = FFT(data,time_step)
x = list_data[0]
y = list_data[1]
plt.figure()
plt.xlim(0,300)
plt.plot(x,np.abs(y)/max(np.abs(y)),'k-+')
"""
N_points = data.size
#FFT
yf_original=np.fft.fft(data*time_step) #*dt really necessary? Better for units, probably
#Post-pro
#We keep only the positive part
yf=yf_original[0:N_points/2]
#fundamental frequency
f1=1/(N_points*time_step)
#Generate the frequency axis - n*f1
xf=np.linspace(0,N_points/2*f1,N_points/2)
return [xf, yf, yf_original]
def Inverse_FFT(data,time_step,freq_cut):
list_data = FFT(data,time_step)
N_points = data.size
#FFT data
x = list_data[0]
yf_full = list_data[2]
#Look where the frequency is above freq_cut
index = np.where(x > freq_cut)
x_new_halfpos = x[0:index[0][0]-1] #Contains N_points/2
yf_new = np.concatenate((yf_full[0:index[0][0]-1], yf_full[N_points/2 +1:index[0][0]-1]))
#Apply inverse-fft
y_complex = np.fft.ifft(yf_new)
#Calculate new time_step ??!!
N_points_new = x_new_halfpos.size *2
f1 = x_new_halfpos[1]
time_step_new = 1/(N_points_new*f1)
#Create back the x-axis for plotting. The original data were distributed every time_step. Now, it is every time_step_new
""" WARNING - It assumes that the new x_new axis is equally distributed - True ?!? """
x_new = np.linspace(0,N_points_new*time_step_new,N_points_new/2)
y_new = y_complex.real / time_step_new
return [x_new,y_new]
生成示例的示例代码：
import numpy as np
import matplotlib.pyplot as plt
#Generate sinus waves
x = np.linspace(0,2*np.pi,50000) #You need enough points to be able to capture information (Shannon theorem)
data = np.sin(x*2*np.pi*5) + 0.5*np.sin(x*2*np.pi*20) + 0.2*np.random.normal(x)
plt.figure()
plt.xlim(0,np.pi/4)
plt.plot(x,data)
time_step = (x[-1]-x[0])/x.size
list_data = FFT(data,time_step)
x = list_data[0]
y = list_data[1]
plt.figure()
plt.xlim(0,30)
plt.xlabel("Frequency [Hz]")
plt.ylabel("Normalized magnitude")
plt.plot(x,np.abs(y)/max(np.abs(y)),'k-+')
#Anti-FFT
data_new = Inverse_FFT(data,time_step,100)
x_new = data_new[0]
y_new = data_new[1]
time_step_new = (x_new[-1]-x_new[0])/x_new.size
plt.figure()
plt.xlim(0,np.pi/4)
plt.plot(x_new,y_new)
list_data_new = FFT(y_new,time_step_new)
x_newfft = list_data_new[0]
y_newfft = list_data_new[1]
plt.figure()
plt.xlim(0,30)
plt.xlabel("Frequency [Hz]")
plt.ylabel("Normalized magnitude")
plt.plot(x_newfft,np.abs(y_newfft)/max(np.abs(y_newfft)),'k-+')
谢谢 ！
亲切的问候，
编辑：更正的功能：
def Anti_FFT(data,time_step,freq_cut):
list_data = FFT(data,time_step)
N_points = data.size
#FFT data
x = list_data[0]
yf_full = list_data[2]
#Look where the frequency is above freq_cut
index = np.where(x > freq_cut)
x_new_halfpos = x[0:index[0][0]-1] #Contains N_points/2
#Fill with zeros
yf_new = yf_full
yf_new[index[0][0]:N_points/2 +1]=0
yf_new[N_points/2 +1 :-index[0][0]]=0 #The negative part is symmetric. The last term is the 1st term of the positive part
#Apply anti-fft
y_complex = np.fft.ifft(yf_new)
#Calculate the """new""" x_axis
x_new = np.linspace(0,N_points*time_step,N_points)
#Divide by the time_step to get the right units
y_new = y_complex.real / time_step
return [x_new,y_new]