python - FFT - 滤波 - 逆 FFT - 剩余偏移

Question

我正在执行以下操作:执行 FFT/削减 FFT 结果中高于 100Hz 的每个频率/执行逆 FFT

如果原始数据集没有偏移量，它效果很好!如果它有偏移，则输出结果幅度会被损坏。

示例:

Without offset

With offset (and noise)

从数学上来说，我什至不确定我是否能做我正在做的事情。我所能观察到的是，在有偏移的情况下，基频是原始频率的两倍？？？!!!

您知道为什么偏移量会改变吗？

代码:

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]

Answer 1

N 个(偶数)实值序列 f 的 DFT(f)(称为 F)具有以下属性:

F(0)(直流偏移)是一个实数

F(N/2) 是实数，即奈奎斯特频率的幅度。对于[1..N/2-1]中的i，F[N/2+i]* = F[N/2-i]，其中“*”表示复共轭。

您对 F 的操作必须保留这些属性。

仅供引用，有一些用于实值输入的特殊例程，它们使用这种对称性来计算 FFT 和 iFFT，其速度几乎是复杂数据的一般对应物的两倍。