БПФ — фильтрация — обратное БПФ — оставшееся смещение

Я делаю следующее: выполняю БПФ / вырезаю каждую частоту выше 100 Гц в результатах БПФ / выполняю обратное БПФ.

Это работает хорошо... если исходный набор данных не имеет смещения! Если он имеет смещение, величина выходного результата повреждена.

Примеры:

Без смещения

Со смещением (и шумом)

Я даже не уверен, что могу делать то, что делаю, с математической точки зрения. Все, что я могу заметить, это то, что со смещением основная частота в два раза превышает исходную ???!!!

Есть ли у вас идеи, почему смещение изменено?

Код:

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]