FFT - фильтрация - Inverse FFT - Offset остальные

Question

я делаю следующее: Выполните FFT/Вырезать каждую частоту выше 100 Гц в результатах FFT/Выполнить Inverse-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 (0), (смещение по постоянному току), является действительным числом

F (N/2) - действительное число, амплитуда частоты Найквиста. Для i в [1..N/2-1] это F [N/2 + i] * = F [N/2-i], где '*' указывает на комплексное сопряжение.

Ваши манипуляции с F должны сохранять эти свойства.

FYI. Существуют специальные подпрограммы для вещественного ввода, которые используют эту симметрию для вычисления FFT и iFFT почти в два раза быстрее, чем их общие аналоги для сложных данных.