2016-11-18 8 views
3

У меня есть сложный сигнал, который я могу видеть в пространстве Фурье и хотел бы фильтровать некоторые частоты, которых я не желаю. Я читал в Интернете, что я должен применить окно Ханнинга, прежде чем принимать преобразование Фурье, чтобы избежать утечки.Восстанавливающий сигнал после фильтра Фурье с окном Ханнинга

Следовательно, то, что я делаю, как видно из приведенного ниже кода, заключается в том, чтобы применить окно Ханнинга к моим данным, а затем принять его преобразование. В качестве теста я хотел посмотреть, не фильтрую ли я что-нибудь, могу ли я вернуть исходный сигнал. Однако сигнал по краям достиг нулевого уровня.

Теперь я понимаю, что это происходит из-за того, что оконный фильтр Ханнинга также заканчивается на нуле. В этом случае, как я могу применить окно Ханнинга, перейдите в частотный домен и вернитесь в мой временной домен с восстановленным сигналом? Если мой сигнал достигает нуля на концах, когда я пытаюсь фильтровать частоты, которые я желаю, результат во временной области будет по-прежнему идти к нулю по краям.

Что мне не хватает/не работает в моем методе? Спасибо за любую помощь!

Вот пример кода, что я делаю:

import sys 
import matplotlib 
def fourier(time,array): 

    fft = np.fft.fft(array*np.hanning(len(array))) 

    Npts = len(array) 
    spacing_array = time[::-1][:-1][::-1] - time[:-1] 

    if np.mean(spacing_array) - spacing_array[0] > 1.e-16: 
      print "time axis not equally separated. cannot compute fft" 
      sys.exit() 
    spacing = spacing_array[0] 

    freq = np.fft.fftfreq(Npts, spacing) 

    return freq,fft 

if __name__ == "__main__": 
    import numpy as np 
    import matplotlib.pyplot as plt 

    # generate a sample signal 
    sample_rate = 100.0 
    nsamples = int(2.e3) 
    t = np.arange(nsamples)/sample_rate 
    x = np.cos(2*np.pi*0.5*t) + 0.2*np.sin(2*np.pi*2.5*t+0.1) + \ 
      0.2*np.sin(2*np.pi*15.3*t) + 0.1*np.sin(2*np.pi*16.7*t + 0.1) + \ 
       0.1*np.sin(2*np.pi*23.45*t+.8) 
    y = np.cos(7*np.pi*1.1*t) + 3*np.sin(0.2*np.pi*2.8*t+1) + \ 
      0.2*np.sin(8*np.pi*2.7*t) + 1*np.sin(2*np.pi*t + 2.1) + \ 
       0.1*np.sin(0.2*np.pi*0.45*t+1.4) 
    z = x*np.exp(y*1j) 

    z_freq,z_fft = fourier(t,z) 

    plt.clf() 
    plt.figure(figsize=(8,12)) 

    plt.subplot(4,1,1) # original signal 
    plt.plot(t,np.absolute(z)) 

    plt.subplot(4,1,2) # fourier transform 
    plt.semilogy(sorted(z_freq),[b for (a,b) in sorted(zip(z_freq,np.absolute(z_fft)/nsamples))]) 

    # filtering 
    plt.subplot(4,1,3) 
    idx = np.where(np.abs(z_freq)>2.0) 
    z_fft[idx]=0 
    z_filter = np.fft.ifft(z_fft) 
    plt.plot(t,np.real(z_filter)) 

    z_freq,z_fft = fourier(t,z_filter) 
    plt.subplot(4,1,4) 
    plt.semilogy(sorted(z_freq),[b for (a,b) in sorted(zip(z_freq,np.absolute(z_fft)/nsamples))]) 
    plt.show() 

И образ, который выходит из него следующий: Real space and Frequency space of signal before and after filtering with Hanning window

+0

Обычно фильтрация в частотной области осуществляется с использованием перекрывающихся блоков и либо [сложение с перекрытием] (https://en.wikipedia.org/wiki/Overlap-add_method) или [перекрытие сохранить] (HTTPS://en.wikipedia.org/wiki/Overlap-save_method) для линейной свертки. –

ответ

1

ножевое на некоторые из ваших вопросов.

В этом случае, как применить окно Ханнинга, перейти в частотную область и вернуться в мою временную область с восстановленным сигналом?

До цифровой точности реализации fft в numpy, применяя преобразование Фурье, а затем обратное преобразование Фурье должно возвращать исходный сигнал. Однако, если вы умножаете сигнал на окно hanning перед fft, результатом fft'ing и ifft'ing будет сигнал раз в окне hanning. Поэтому третий сюжет, который вы предоставили, имеет смысл. Это исходный сигнал, умноженный на окно сглаживания (и сглаживается путем исключения частот выше 2.0).

Я прочитал онлайн, что перед применением преобразования Фурье я должен применить окно Ханнинга, чтобы избежать утечки.

Окно ганнинга используется для минимизации спектральной утечки. В идеальном случае вы можете использовать свой конечный клип данных в виде плитки, чтобы воссоздать всю функцию во всем пространстве. Когда ваш конечный клип данных не может точно служить этой цели, вы понесете некоторое количество спектральной утечки.

Окно укрощения ослабляет начало и конец сигнала, чтобы помочь удовлетворить ограничение «повторяющейся плитки». При этом я не думаю, что спектральная утечка является проблемой на данных, которые вы предоставили, но, возможно, это в другом наборе данных, который вы анализируете.

Для получения дополнительной информации о спектральной утечке, перейдите к why does spectral leakage arise in a fft.

Если мой сигнал заканчивается на концах, когда я пытаюсь фильтровать частоты, которые я желаю, результат во временной области будет по-прежнему идти к нулю по краям.

Это не обязательно так. Возможно, что частоты, которые вы устраняете в частотной области, фактически были режимами, которые привели функцию к нулю во временной области.

Четвертый участок

Ваш четвертый участок интересен. В нем вы рисуете результат преобразования Фурье высокочастотных исключенных временных рядов. Функция быстро распадается на +/- 2, потому что вы устранили эти частоты. Тот факт, что fft не идет ровно в 0, состоит в том, что дискретное преобразование Фурье не является точным.

Быстрое кодирование комментарий

Одна строка кода имеет интересное применение массива нарезания. Это заняло у меня секунду, но я хотел прокомментировать то, что вы сделали.

time[::-1][:-1][::-1] # Reverse array, disregard "new" last element, ... 
         # ... then reverse again 
time[1:] # Don't consider the first element. 
time[1:] == time[::-1][:-1][::-1] # These are exactly the same. 
# Result: True 

 Смежные вопросы

  • Нет связанных вопросов^_^