.. index:: filter

===================
フィルター
===================


.. toctree::
   :maxdepth: 1
   
　前ページの周波数解析で，雑音がある場合，雑音が大きくなりすぎると周期的成分が雑音に埋もれて，
あるはずのピークが周波数解析により見いだせなくなった．
サンプリング周波数を十分にあげれば何とか見つけることもできたが，
元の信号の成分を知らない時にどうやって十分と判断すればよいのか．
実験環境によってはサンプリング周波数を上げられないこともある．


フィルター
==============

フィルター（濾過器）によっていらない信号を取り除こう．
まずは，二つの周波数からなる信号の一部を取り除く．
前ページと同じく :math:`f(x)=1.0\sin(0.1*2\pi*x)+1.0\cos(1.0*2\pi*x)` より得られる信号があった時に，
たとえば後半の1Hzの信号を取り除きたいとする．

フーリエ変換と逆フーリエ変換を利用する．

+ フーリエ変換（時間→周波数）を行う．
+ 得られた周波数信号のうちいらない周波数成分をゼロにする．
+ 逆フーリエ変換（周波数→時間）を行う．いらない周波数成分が取り除かれた信号が得られる．

.. code-block:: python
    :linenos:
    :emphasize-lines: 9,28 

    import numpy as np
    import matplotlib.pyplot as plt 
    from scipy import fftpack

    def f(x):
        pi2 = 2.*np.pi
        return 1.0*np.sin(0.1*pi2*x) + 1.0*np.cos(1.*pi2*x)

    low_pass_freq = 0.5
    num_of_sample = 1024
    sampling_time = 100.0
    sampling_freq = num_of_sample /sampling_time
    print 'sampling_freq is',sampling_freq,'[Hz]'

    #Periodic data with/without random noise
    xdata = np.linspace(0, sampling_time, num=num_of_sample)
    np.random.seed(1234)
    ydata = f(xdata) #+ 1.0*np.random.randn(xdata.size)

    #FFT
    time_step = xdata[1]-xdata[0]
    sample_freq = fftpack.fftfreq(ydata[:].size, d=time_step)
    y_fft = fftpack.fft(ydata[:])
    pidxs = np.where(sample_freq > 0)
    freqs, power = sample_freq[pidxs], np.abs(y_fft)[pidxs]
    freq = freqs[power.argmax()]
    #iFFT
    y_fft[np.abs(sample_freq) > low_pass_freq] = 0
    filtered_ydata = fftpack.ifft(y_fft) 
    filtered_ydata = filtered_ydata.real
    #FFT2
    y_fft2 = fftpack.fft(filtered_ydata[:])
    pidxs2 = np.where(sample_freq > 0)
    freqs2, power2 = sample_freq[pidxs2], np.abs(y_fft2)[pidxs2]
    freq2 = freqs[power2.argmax()]

    #PLot
    plt.figure(figsize=(8,10))
    plt.subplot(211)
    plt.plot(xdata,ydata,'b-', linewidth=0.5)
    plt.plot(xdata,filtered_ydata,'r-',linewidth=1, label='filtered') 
    plt.xlim(0.,20.)
    plt.xlabel('Time')
    plt.ylabel('Ydata')
    plt.grid(True)

    plt.subplot(212)
    #plt.semilogx(freqs, power,'b.-',lw=1)
    plt.loglog(freqs, power,'b.-',lw=1)
    plt.loglog(freqs2, power2,'r.-',lw=1)
    plt.ylim(0.1,10.**4)
    plt.xlabel('Frequency')
    plt.ylabel('Power')
    plt.grid(True)

    plt.show()

※出だしの書き方が前ページから少し変更．

９行目：ここで指定する値以下の周波数成分だけを残すことにする．0.5Hzとした．

２８行目：low_pass_freq（0.5Hz）より大きい周波数成分をゼロにする．
２９，３０行目：逆フーリエ変換を行い，時間領域データを得る．
３２行目～：さらに得られたデータをフーリエ変換し，周波数領域でフィルターの効果を確認．

実行結果は以下の通り，

.. figure:: ../img/230_filter.png
   :scale: 70 %
   :alt: filter

上の図の青い線が元の信号で，赤い線が0.5Hz以下の成分だけを残した信号である．
（ある値以下の信号だけを残すフィルターをLow pass filter（低い成分だけが通過できるフィルター）という）
0.1Hzのきれいなsinカーブである．
下の図は，周波数領域での表示．フィルター通過後の赤い線では，カットした0.5Hz以上の成分がほとんどなくなっているのが分かる．

しかし，やはり前述のエイリアス現象には注意しなくてはいけない．
サンプリング数（num_of_sample = 1024）が不十分なまま，cosの周波数を10Hzにすると，

.. code-block:: python

    return 1.0*np.sin(0.1*pi2*x) + 1.0*np.cos(10.*pi2*x)

.. figure:: ../img/231_filter_with_alias.png
   :scale: 70 %
   :alt: filter with alias

やはりうまくいかないのである．おそるべしエイリアス．
もちろん，サンプリング数を倍にしてサンプリング周波数が10Hzの2倍になるようにすればうまくいく．

アンチエイリアスフィルター
=================================

そのうち．