6. モジュールの使い方

6.1. スクリプトとモジュール

いままでは，主にインタープリタに直接コードを書いてきましたが，一連の処理をファイルとしてまとめることなどを考えてみます．この一連の処理をまとめたものを スクリプト と呼び，再利用できるように設計されたスクリプトを モジュール と呼びます．

6.1.1. import 構文によるモジュールの読み込み

import 構文によるモジュール読み込みを考えてみましょう．モジュールは python を便利に使えるように機能を集めたもので，車輪の再発明を行わずに python を便利に使うためには重要な仕組みです．

例えば，ディレクトリにどんなファイルがあるかを知りたいといた要求があったとしましょう． python では， os モジュールとして標準的に提供されています．

さて，モジュールからオブジェクトをインポートするには以下の様な文を用います．

import os

これで読み込みは終了です． os がどのようなものかをインタープリタ上で表示させると以下のような表示が得られるます．

In [1]:  os
<module 'os' from '/usr/local/lib/python3.9/os.pyc'>

このように os は，モジュールと言われることがわかります．(バージョンやパスなどは異なる場合があります) 細かい使い方は help(os) とかやると出てきますが，ここでは os モジュール中に定義された「ディレクトリの中のファイル名を得るメソッド」:func:os.listdir() を使ってみます． import で読み込んだモジュールはオブジェクトとして取り扱われるので， . を用いてアクセスします．今， os という名前でオブジェクトを扱っていますので， os.listdir() メソッドにカレントディレクトリ ./ の内容を表示させてみます．

In[1]: os.listdir('.')

['ファイル1.txt',
 'ファイル2.txt',
 'index.rst']

などといった現在スクリプトが動いているディレクトリの ファイル名のリスト が得られます．

また，特定の機能（この場合は os.listdir() ) のみが欲しい場合は，以下のような構文で import できます．この場合はいきなりメソッド名を指定します:

In[1]: from os import listdir
In[2]: listdir('.')

（完全に英語の構文ですね．）さらにモジュール名が長いときは，略記してインポートが可能です．略記でインポートするには

In [1]: import numpy as np       # numpy を np としてインポート
In [2]: np.mean( [1, 2, 4] )     # リストの平均求めるメソッド
2.3333333333333335

こういう省略形で import するのはよくやりますが，やりすぎると訳がわからなくなりますので注意してください．

6.1.2. 名前空間に関して少し

注釈

from os import *

という書き方をする解説もあったりします，が，あまり推奨されません．アスタリスク * はワイルドカードを表し，すべての関数を表しますので， os に含まれるすべてのモジュールや関数を，現在作業しているスペースに展開してしまいます．

これをスクリプト内部で書くと以下の点でややこしい状況になります．

まず os モジュールに含まれる全ての名前が展開されてしまうので，自分で定義する関数とぶつかる（conflict）可能性が高くなります．また複数のモジュールを import した時に同じ名前の関数があると，ぶつかってしまいます．

また名前と文脈から機能を予測することが困難になります．例えば， os.name は OS の名前ですが，これがワークスペースに展開され name という名前で使えるようになっちゃうわけです．こうなると name という名前の変数は一般的すぎて何を意味するのかがわからなくなりがちです．

モジュールを import した場合，クラスの名前は名前空間と呼ばれる領域に展開されます．これは，import したときに os という名前の階層の下にメソッドを展開することを意味します．名前空間がどのようになっているのかを知るには， dir() 関数を使います． dir() 関数は名前空間に含まれるメソッド，変数名をリストにして表示してくれます．

通常の python を起動した段階で，dir() 関数を呼び出すと

In [1]: dir()
OUt[1]: ['__builtins__', '__doc__', '__name__', '__package__']

という名前が展開されていて，これ以外のものがないのがわかります．ここで，変数やリストを定義したり，モジュールをインポートしてみて，名前空間がどのように変化するのかを確認してみてください．

In [1]: a = 1
In [2]: s = 'abc'
In [3]: import os
In [4]: dir() # 名前空間を表示
Out[4]: ['__builtins__', '__doc__', '__name__', '__package__', 'a', 'os', 's']

といったように，変数名やモジュール名が登録されているのがわかります．これに

In [5]: from os import *
In [6]: dir()

とやると，名前空間に os モジュール内部のメソッド名がずらりと展開しているのがわかります．これらの名前と名前がかち合わないようにプログラムを書くのは結構大変です．

したがって，モジュールは，名前を使って階層的に管理するのが上手なやり方となります．下記では numpy を np として， scipy を sp としてます．

In [1]: import numpy as np        # data arrays
In [2]: np.linspace(0, 10, 6)     # 0 から 10 までを等間隔で6個の数値を生成
array([  0.,   2.,   4.,   6.,   8.,  10.])
In [3]: import scipy as sp
In [4]: x = np.linspace( 0, 10, 16 )
In [5]: X = sp.fft.fft( x )

このように，モジュール名.メソッド名のように指定すれば，どのモジュールに属するメソッドかを指定して呼び出すことができます．

6.2. Numpy によるデータ行列の取り扱い

Python では，配列は，リストとして扱うのですが，若干 トロい のが問題です． Python で行列や配列を扱うためのモジュール NumPy は，この弱点を解消すべく Python 上での数値演算に特化した配列を扱うためのモジュールです．

以降では， numpy, scipy, matplotlib モジュールが import されていることを仮定しています．（なのでスクリプト書くときは最初の方におまじないとして書いておいてください）

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt

ここでは， NumPy で最も重要なクラスである np.ndarray について，最低限必要と思われる知識について説明します．

np.ndarray は， N-dimensional Array すなわち， N 次元配列を扱うためのクラスです．通常 Python ではこうしたN次元配列を表現するには，多重のリストが利用されます．多重リストは，あらゆるを入れられるコンテナとしては便利なのですが，配列の処理速度が犠牲になっている面があります． np.ndarray は，この処理速度を向上させるためのものです．下記に多重リスト np.ndarray との違いを列挙します．

多重リスト
- リンクでセルを結合した形式でメモリ上に保持される．なので，多重リストは動的に変更可能．
- リスト内部で要素の型が異なっても良い
- 多重リストは，各次元の要素数が異なっても OK．なので行ごとに列が異なる2次元配列なども扱える
- 全体を操作しようとすると，若干トロい
np.ndarray
- C や Fortran の配列と同様にメモリの連続領域上に保持される．なので，np.ndarray の形状変更には全体の削除・再生成が必要．
- np.ndarray は， 全て同じ型 の要素で構成されていなければならない例えば実数は64ビットの倍精度浮動小数(np.float64)が標準となるが，単精度浮動小数(np.float32)と混ぜて使うことはできない．（Python の標準のリストはこれができる）
- np.ndarray は各次元ごとの要素数が等しくなければならない．行ごとに列数が異なる2次元配列は扱えない．リストのリストはこれができる．
- np.ndarray は，行や列を対象とした多くの高度な数学的操作を，多重リストより容易かつ高速に適用できる．

6.2.1. NumPy での配列の生成

N 次元配列 np.ndarray は，数学の概念で言えば， 1次元の場合はベクトルに，2次元の場合は行列に，そして3次元以上の場合はテンソルに該当します．ここでは np.ndarray の生成方法を説明します．

6.2.1.1. `np.array()` 関数による生成

np.ndarray にもコンストラクタはありますが，通常は， np.array() 関数などによって生成します．

np.array(object, dtype=None): numpy の配列オブジェクト np.ndarray を生成する．引数としては配列を表現する object を与える必要がある．

最初の引数 object には，配列の内容を，array_like という型で与えます．この array_like という型は，配列を np.ndarray の他，（多重）リストや（多重）タプルで表現したものです．

要素が 1, 2, 3 であるリスト [1, 2, 3] を与えて，長さ 3 のベクトルを生成するの例です:

In [1]: a = np.array([1, 2, 3])
In [2]: a
Out[2]: array([1, 2, 3])
In [3]: type(a)
Out[3]: numpy.ndarray

タプルを使った表現も可能です:

In [1]: a = np.array((10, 20, 30))
In [2]: a
Out[2]: array([10, 20, 30])

タプルを引数として与える場合は，関数の括弧と勘違いしないように注意しましょう．

次に2重にネストしたリストで表した配列の例です:

In [1]: a = np.array([[1.5, 0], [0, 3.0]])
In [2]: a
Out[2]:
array([[ 1.5,  0. ],
       [ 0. ,  3. ]])

リストの要素に np.ndarray やタプルを含むことも可能です．ただし，各次元の要素数が合っている必要があることに注意する必要があります．

In [1]: a = np.array([1.0, 2.0, 3.0])
In [2]: b = np.array([a, (10, 20, 30)])
In [3]: b
Out[3]:
array([[  1.,   2.,   3.],
       [ 10.,  20.,  30.]])

6.2.1.2. その他の numpy 関数による生成

実際には値を直接生成することはあまりしません．ここではいくつかの生成方法を例を示しながら紹介します． np.ndarray を作るための関数は， np.array() 以外にも数多くありますが，それらのうちの有名なものを記載します．

np.zeros() は，要素が全て 0 である0行列を生成

In [1]: c = np.zeros((2,1))
In [2]: c
Out[2]:
array([[ 0. ],
       [ 0. ]])

np.ones() は，要素が全て 1 である1行列を生成

In [1]: a = np.ones( (3, 3) )
In [2]: a
Out[2]:
array([[1., 1., 1.],
       [1., 1., 1.],
       [1., 1., 1.]])

np.eye() は，単位行列を生成

In [1]: e = np.eye( 2 )
In [2]: e
Out[2]:
array([[ 1., 0. ],
       [ 0., 1. ]])

np.arange() , np.linspace() は，順番にシーケンスを生成

In [1]: f = np.arange(5)
In [2]: f
Out[2]: array([0, 1, 2, 3, 4])
In [3]: g = np.linspace(0, 1, 6)  # [0, 1] を 6 個の等間隔点で
Out[3]: g
Out[3]: array([ 0. ,  0.2,  0.4,  0.6,  0.8,  1. ])

6.2.2. NumPy 配列の属性

ここでは，前節で生成した np.ndarray の属性を説明したのち，配列の要素を参照する方法について述べます．

np.ndarray には多数の属性がありますが，よく使われるものをまとめると以下のようなものになります．

class np.ndarray

固定長の多次元配列クラス．

変数:

dtype -- 各要素のデータ型
ndim -- データの次元数
shape -- 各次元の長さをタプルで表したもの

最初の属性 dtype は配列の要素の型を表します． np.ndarray は，配列の中の全要素の型は基本的に同じです．
二番目の属性 ndim は，次元数を表します．ベクトルでは 1 に，配列では 2 になります．
三番目の属性 shape は，各次元ごとの配列の大きさをまとめたタプルで指定します．例えば，長さが 5 のベクトルは (5,) となり， $2 \times 3$ 行列では (2, 3) となります．(1次元の時の表記には注意してください)

In [1]: a = np.array( [[1, 2, 3], [4, 5, 6]] )     # 2 次元 (2,3) 行列を生成
In [2]: a
Out[3]:
array([[1, 2, 3],
       [4, 5, 6]])
In [4]: a.ndim                                     # a の次元数
Out[4]: 2
In [5]: a.shape                                                        # a の形状
Out[5]: (2, 3)

これらの属性のうち dtype では，Python のビルトイン型の真理値型，整数型，浮動小数点型，複素数型に対応する np.bool ， np.int ， np.float ， np.complex が良く使われます．

配列の dtype 属性を指定するには，(1) np.array() などの配列生成関数の dtype 引数で指定する方法と，(2) np.ndarray の np.ndarray.astype() メソッドを使う方法とがあります．

ここでは，(1) の dtype 引数を指定する方法について述べます． np.array() では要素が全て整数の場合は，要素の型は整数になりますが，それを浮動小数点にするには，次のように指定します:

In [1]: a = np.array([1, 2, 3])
In [2]: a.dtype
Out[2]: dtype('int64')
In [3]: a = np.array([1, 2, 3], dtype=np.float)
In [4]: a.dtype
Out[4]: dtype('float64')

6.2.3. NumPy 配列の参照

ここで np.ndarray の要素の参照方法について述べます，が，割りとおおくの要素の参照方法があるため，ここでは基本的な方法を述べます．最も基本的な要素の参照方法とは，各次元ごとに何番目の要素を参照するかを指定します． 1次元配列であるベクトル a の要素 3 である a[3] を参照すると，次のような結果が得られます．

In [1]: a = np.array([1, 2, 3, 4, 5], dtype=float)
In [2]: a[3]
Out[2]: 4.0

添え字の範囲は， Python の規則に従って $0,\ldots,4$ であることには注意してください． a.shape[0] とすると，第1次元の要素の長さ，すなわちベクトルの長さとして 5 が得られます．

同様に， $2 \times 3$ の行列では，行は $0,\ldots,1$ の範囲で，列は $0,\ldots,2$ の範囲で指定します．

In [1]: a = np.array([[11, 12, 13], [21, 22, 23]])
In [2]: a.shape
Out[2]: (2, 3)
In [3]: a[1, 2]
Out[3]: 23

さらに，インデクスのスライスをならべて参照することが可能です．

In [1]: a = np.arange(10)
In [2]: a
Out[2]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
In [3]: a[2:9:3]   # [start:end:step]
Out[3]: array([2, 5, 8])

多次元にも対応しています．図としてまとめるとこんな感じで参照できます．

図は scipy-tutorial からの引用です．

6.2.4. Numpy 配列形状の変更

形状の変更は np.ndarray.reshape() を使うのが一般的です．例えば，0 から順番に 11 までの要素を持つ1次元配列を生成し，これを $3\times 4$ 行列に変更するには，以下のようにします． reshape() は形状を変化させるだけです．

In [1]: a = np.arange(12)
In [2]: a.shape
Out[2]: (12,)                      # <--- 1次元で12個の要素からなる
In [3]: b = a.reshape((3, 4))  # 引数がタプルなのに注意
In [3]: b                      # b を確認する
Out[3]:
array([[ 0, 1, 2, 3],
       [ 4, 5, 6, 7],
       [ 8  9,10,11]] )
In [4]: b[0,0] = 99
In [5]: b                      # b を確認再び
Out[5]:
array([[99, 1, 2, 3],
       [ 4, 5, 6, 7],
            [ 8  9,10,11]] )
In [6]: a               # a を確認
Out[6]: array([99,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11])

このように reshape は，あくまで行列の見え方 (ビューと呼びます) を変化させるだけのものです．もし，行列のコピーが必要であれば，

cp_a = np.array(a)

のように， np.array() を使って新たなオブジェクトを作成し，それの初期値に a を与えるようにして コピー を作成してください．

6.2.5. Numpy 配列の結合

Numpy どうしの配列をくっつけてみることを考えましょう．複数の np.ndarray は，いろいろな方向に延ばすことができます．これには np.hstack() や np.vstack() といった関数を用いることで実現できます．

>>> a = np.arange(12).reshape(4, 3)
>>> b = np.arange(8).reshape(4, 2)
>>> c = np.arange(4).reshape(4, 1)

といった 3 個の np.ndarray は，ともに4行であるという性質を持ちます．これを横方向に並べた配列を作るには np.hstack() 関数を使います．

>>> d = np.hstack((a, b, c))
>>> d
array([[ 0,  1,  2,  0,  1,  0],
      [ 3,  4,  5,  2,  3,  1],
      [ 6,  7,  8,  4,  5,  2],
      [ 9, 10, 11,  6,  7,  3]])

のようにくっつけることが可能になります．ただし，これは配列のくっつける方向の次元が同じでないといけないので，配列の形は気にしておいてください．

この他にも np.concatenate(), np.stack() , np.block() などがあります．

6.2.6. NumPy 配列のファイルへの読み書き

せっかく作ったデータはなるべくファイルに保存しましょう．

Numpy 配列のデータをファイルにデータ配列を書き込むには， np.save() もしくは np.savetxt() を使います．テキストファイルとして，エクセルとかで読み書きしたい場合は，後者をつかってください．

テキストファイルへの書き込み

>>> a = np.arange(25).reshape(5,5)    # 0〜24 を 5x5 の行列へ
>>> np.savetxt('dataA.txt', a)

バイナリファイルへの書き込み
1>>> b = np.arange(50).reshape(10,5) # 0〜49 を 10x5 の行列へ 2>>> np.save('dataB.npy', b)
バイナリファイルの拡張子は，どのように設定しても問題ないですが，人間が識別しやすいように .npy を指定しておくのがよいでしょう．

読み込みに関しては np.load() もしくは np.loadtxt() を使ってください．

テキストファイルからの読み込み

>>> c = np.loadtxt('dataA.txt')    # 保存したファイルから読み込み

バイナリファイルからの読み込み

>>> d = np.load('dataB.npy')     # 保存したファイルから読み込み

実験を途中で終わるときなどは，配列をファイルに保存しておき，次の実験時にそのファイルを読み込むことで，続きを行うことができます．

6.3. Numpy / Scipy / Matplotlib 連携による簡単な可視化操作

データ配列ができたら，データの操作と可視化をしてみましょう．これには matplotlib を使います．すでに matplotlib.pyplot が plt という名前で import されていることを仮定します．

注釈

jupyter では必要ないですが，コマンドラインで ipython を使っている場合は， ipython --pylab として，matplotlib を展開するのも一つの手です．
1$ ipython --pylab
2In [0]:

この状態では， matplotlib.pyplot のメソッドが，名前空間の一番の大元の名前空間に展開されるので，以下で用いているの plt. から始まる接頭語部分はつけても付けなくてもよいです．つまり，ipython を --pylab オプションつきで立ち上げると

1>>> from pylab import *
2>>> from matplotlib.pyplot import *
3>>> import numpy as np
4>>> import matplotlib.pyplot as plt
的なことをしてくれています．（他にもいろいろやってるけど）

6.3.1. 1 次元データの可視化: グラフ

たとえば，横軸を $y = x^2$ のグラフを書いてみましょう．グラフの横軸の範囲を 0 〜 19 とします．この時は，横軸に a 縦軸に a**2 の点をデータとして生成してそれをプロットすることになります．

>>> a = np.arange(20)       # a は1次元の配列で 0 〜 20
>>> plt.plot(a, a**2)       # 線の描画
[<matplotlib.lines.Line2D object at 0x95abd0c>]

全ての描画コマンドが終わったら，グラフを表示させるために

>>> plt.show()

を実行してください．ウィンドウの中にグラフが描画されるはずです．なお，Colab や Jupyter Notebook を使っている場合は， plt.show() を使わなくとも表示されるはずです．

注釈

普通の python インタープリタでは，グラフのウィンドウを消すまでは，コマンドプロンプトが戻らないことに注意してください．

グラフのウィンドウを一旦閉じて，サンプル点や， x 軸とか y 軸のラベルを貼りこむには，下のようにします．

>>> plt.plot(a, a**2, 'o' )  # 線を描画 とサンプル点描画
>>> plt.grid()              # グリッド描画
>>> plt.xlabel( 'x' )       # ｘ軸のラベル
>>> plt.ylabel( 'y' )       # y 軸のラベル
>>> plt.show()              # 画像の描画

注釈

plt.show() を忘れると描画しないことがあるので注意しましょう．

また，グラフを消去するには

>>> plt.clf()

を実行してください．

6.3.2. 2 次元データの可視化: 画像

次に画像を扱ってみましょう．画像は2次元の行列状の配列として扱えます．行列の各要素にに色の情報を載せれば（ピクセル）として扱えます．

例えば $3\times 3$ 程度の単位行列を可視化してみるとこんな感じになります．

 >>> e = np.eye(3)          # 単位行列の生成
 >>> plt.imshow(e)
 <matplotlib.image.AxesImage at 0x8616269d0>
 >>> plt.show()

とすれば，おそらくこんな画像が得られます．

ただし，pyplot のイメージ表示は，行列で記述した画像が自然に見えるようにかなりぼやかして表現してしまいます．それが嫌な人は plt.imshow() 関数の interpolation 引数に 'Nearest' もしくは 'none' を与えてください．

 >>> plt.clf()
 >>> plt.imshow( e, interpolation='none' )
 <matplotlib.image.AxesImage at 0x8616269d0>
 >>> plt.show()

こうすれば，下の画像が得られます．

さらにもう少し大きな行列 $32\times 32$ くらいを生成し，各要素をランダムな値に設定して表示させてみます．その上で，色のテーブルを，グレイスケールや，hot と呼ばれるテーブルを使って表示させると，こんな感じになります．

>>> img = np.random.rand(32,32)
>>> plt.imshow(img)
<matplotlib.image.AxesImage object at 0x9e954ac>
>>> plt.gray()
>>> plt.show()


>>> plt.imshow(img)
<matplotlib.image.AxesImage object at 0x9e954ac>
>>> plt.hot()
>>> plt.show()

6.3.3. 複数グラフの同時表示

複数のグラフを同時表示させるには subplot() を使います．例えば，上記のランダムなグラフを３つ並べて色の違いを調べたいような場合には以下のようにします．

>>> img = np.random.rand(32,32)
>>> plt.subplot( 1, 3, 1 )       # 表示区画を１行３列に設定し，その１番目に描画
>>> plt.imshow( img )

>>> plt.subplot( 1, 3, 2 )               # 表示区画の２番目を指定
>>> plt.imshow( img )
>>> plt.gray()

>>> plt.subplot( 133 )                   # 表示区画の３番目を指定
>>> plt.imshow( img )
>>> plt.hot()

>>> plt.show()

6.3.4. その他

matplotlib には他にも多くの機能があります：

色の選択
マーカーによる注釈
LaTeX のフォント
図の取り込み
頻度分布

等の機能が用意されています．

より多くの情報は：

matplotlib documentation http://matplotlib.sourceforge.net/contents.html

an example gallery with corresponding sourcecode http://matplotlib.sourceforge.net/gallery.html

を参照してみてください．

6.4. フーリエ解析をやってみる

6.4.1. 準備

さて，ここでは numpy , scipy, matplotlib を用いてデータ処理を試してみましょう．

ここでも

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt

が行われていることを仮定しています． 関数が無いよ 的なエラーが出るときは， numpy, scipy, matplotlib が正しくインポートされているか確認してください．また，自分でモジュール設計する実験になるので，展開されているソースフォルダを環境変数 PYTHONPATH に指定することを忘れないでください．

おおよそ，世の中の人が必要としている科学技術計算の基本部品は numpy, scipy の中に用意されています．なので，ポイントは，自分が処理したいデータを

どんな形に整形して入力する必要があるのか，

どんな形で出力されてくるのか

をサンプル例，文書などから把握し，適用していくことです．

題材としては， フーリエ変換 を用いたデータ処理扱います．数学的な定義は他の科目に任せることとして，ここでは，具体的なデータ処理を考えていきます．

とはいえ，数学的な定義は理解に必要なので，簡素に書いておきましょう．時間 $t$ で定義された関数 $x(t)$ の中にどんな波の成分 (周波数 $f$ の波) が含まれているかを解析するのが Fourier 変換です．この時の波の成分の方を $X(f)$ という関数で表すとすると Fourier 変換対は

\[ \begin{align}\begin{aligned}X(f) = \int dt x(t) e^{-2\pi j f t}\\x(t) = \int df X(f) e^{ 2\pi j f t}\end{aligned}\end{align} \]

と表されます．

6.4.2. scipy を用いたフーリエ変換(fft)

計算機では，データは有限個の範囲でしか扱えないので，積分の厳密な計算はできません．計算機でフーリエ変換を扱う場合，以下のような，離散フーリエ変換 (Discrete Fourier Transform: DFT) を考えます．

まず，データ $x(t)$ を区間 $[0,T]$ からサンプリング間隔 $\Delta$ で $N$ 個のデータを採ってきて $x_n = x(n\Delta)$ と表します．このときの離散フーリエ変換対は

\[ \begin{align}\begin{aligned}X_m = \sum_{n=0}^{N-1} x_n e^{-2\pi j \frac{mn}{N} }\\x_n = \sum_{m=0}^{N-1} X_m e^{ 2\pi j \frac{mn}{N} }\end{aligned}\end{align} \]

となります．（周波数データ $X_m$ を実際の（物理的な）周波数関数のサンプルとして表現するには， $X(m \tilde\Delta) = \Delta X_m$ として表されます．ただし $\tilde\Delta = \frac{1}{N\Delta}$ です．）

以上は，定義ですが，わかることは，時間データ $x_n$ をフーリエ変換すると周波数データ $X_m$ が（複素数の）数列で得られるということです．

それでは，適当に試してみましょう．最初は簡単な例から考えます

 In [1]: N = 16
 In [2]: t = np.arange(N)
 In [3]: x = np.sin( 2 * np.pi / N * t )

として，16個のサンプル点の中での正弦波を考えます．このデータを matplotlib の plot() を用いてプロット表示させてください． 16個のサンプルで波が一個含まれることがわかります．

次にこの信号を Fourier 変換してみましょう． Fourier 変換は scipy の fft というモジュールにありますので import` しておきます．なお，fft という語は，高速フーリエ変換 (Fast Fourier Transform) の略です．

from scipy.fft import fft, iff

これで，高速フーリエ変換の関数 fft() とその逆変換 ifft() とが使えるようになります．なお，実数配列しか使わない場合は，ちょっとだけ高速化された rfft() や irfft() を使うのもありです．ためしに，先ほど定義した x をフーリエ変換してみましょう

In [1]: X = fft( x )
In [2]: X
Out[2]:
array([ -9.95799250e-17 +0.00000000e+00j,
        -1.01064310e-15 -8.00000000e+00j,
        -4.27067681e-16 -6.60554663e-16j,
         9.95799250e-17 -4.44089210e-16j,
         2.33486982e-16 -5.55111512e-16j,
         9.95799250e-17 -4.44089210e-16j,
         6.71997040e-16 -4.38510058e-16j,
         3.21624530e-16 +0.00000000e+00j,
         1.22464680e-16 +0.00000000e+00j,
         3.21624530e-16 -0.00000000e+00j,
         6.71997040e-16 +4.38510058e-16j,
         9.95799250e-17 +4.44089210e-16j,
         2.33486982e-16 +5.55111512e-16j,
         9.95799250e-17 +4.44089210e-16j,
        -4.27067681e-16 +6.60554663e-16j,
        -1.01064310e-15 +8.00000000e+00j])

のような16個の複素数からなる要素が表示されると思います．これだけだとなんだかわからないので，プロット等で 可視化 します．この場合，複素数なので，絶対値や，実部，虚部などを表示してみましょう．

plt.plot(np.abs(X), 'bo-')    # 絶対値（パワースペクトル）を青点と青線で
plt.plot(np.real(X), 'go-')   # 実部を緑で
plt.plot(np.imag(X), 'ro-')   # 虚部を赤で

横軸（周波数軸）が 1 の位置と, 15 (N-1であることに注意) の位置に虚部に値を持っていて，他の部分がほぼ0であることがわかります．このことから周波数 1 の成分を取り出していることが予想されます．

さらに， fft がどういうデータを表現しているかを理解するために，もう少し例題で遊んで見ましょう．

注釈

新しいモジュールや関数を使うときは，必ず 簡単な（答えのわかっている）例題を使って実験 してみましょう． "こう動くはず" という思い込みは，バグのもとになります．

周波数を倍にしてみましょう．

f = 2
x = np.sin(2 * np.pi / N * f * t)

として同じことをすると

のようなグラフが得られます．この時は横軸 2 の位置と，14 (N-2の位置であることに注意)の位置に値をもっています．すなわち周波数が 2 （N 中に２周期分の波）の波が入っているということを示しているようです．

さらに，

f0 = 1
f1 = 3
x = 0.3 * np.sin(2 * np.pi / N * f0 * t) + 0.7 * np.sin(2 * np.pi / N * f1 * t)

のように二つの波を混ぜてみます．$x(t)$ とその Fourier 変換 $X(f)$ （のパワースペクトル）とをプロットすると下記の様になります．

周波数 1 の位置と，3の位置に，混合した比率に応じた成分が元の信号に入っていることがわかります．

6.4.3. 周波数解析の限界

それでは，周波数 f はどこまで記述できるかを考えてみましょう？ f = 7 の場合を考えてみよう．

f = 7
x = np.sin(2 * np.pi / N * f * t)
X = fft(x)

として実験してみた場合と，f=9

f = 9 x = np.sin(2 * np.pi / N * f * t) X = fft(x)

として実験してみた場合とを比較してみます．これらの場合の Power Spectrum は，以下のようになり，同一になることがわかります． (つまり f=8 を超えた波は正しくは表現できない)

なお，この限界周波数は， ナイキスト周波数 として知られます．この場合は，配列全体の長さが (N=16) で，その半分の f=8 が，このナイキスト周波数にあたります．

6.4.4. scipy を用いた逆フーリエ変換( `ifft()`)

ifft() に関しても使う分には簡単です．

N = 16
t = np.arange( N )
x = np.sin(2 * np.pi / N * t)   # 原信号を作成
X = fft(x)                # X にフーリエ変換した信号
y = ifft(X)               # y は逆フーリエ変換

とすれば，y に逆フーリエ変換が求まっているはずです．実際に y を表示してみると下記のような感じに成るはずです

In [1]: y
Out[1]:
array([  0.00000000e+00 +0.00000000e+00j,
         3.82683432e-01 +0.00000000e+00j,
         7.07106781e-01 -5.55111512e-17j,
         9.23879533e-01 -8.32667268e-17j,
         1.00000000e+00 +0.00000000e+00j,
         9.23879533e-01 -2.77555756e-17j,
         7.07106781e-01 -5.55111512e-17j,
         3.82683432e-01 +2.46519033e-32j,
         1.22464680e-16 +0.00000000e+00j,
        -3.82683432e-01 +0.00000000e+00j,
        -7.07106781e-01 +5.55111512e-17j,
        -9.23879533e-01 +8.32667268e-17j,
        -1.00000000e+00 +0.00000000e+00j,
        -9.23879533e-01 +2.77555756e-17j,
        -7.07106781e-01 +5.55111512e-17j,
             -3.82683432e-01 +2.46519033e-32j])

注意しなければならないことは， ifft() 関数は複素数の列を返してきます．虚部の値はものすごく小さい値ですが，型としては複素数です．グラフとかに表示したりするときには，実部のみを表示してみてください． (あるいは irfft() を使うか．)

6.4.5. ここまでのまとめ

計算機の上でフーリエ変換がどのように記述されるのかをまとめます．

フーリエ変換 scipy.fft.fft() は，与えられた信号 x から周波数成分を取り出して複素数列 X を吐き出す．

周波数信号 X から，原信号列 x を得るには，逆フーリエ変換 scipy.fft.ifft() を用いる

フーリエ変換 X のパワースペクトルは，周波数毎の強度に対応する値の列である．

長さ N の実数列を変換した場合，パワースペクトルは N/2 を中心として対象になる．

ナイキスト周波数 N/2 を超えた表現は得られない．

現実時間の音声信号などは，ある一定の周期のもとデータとしてとられ，データ列 x として計算機に取り込まれます（サンプリング）．この一定周期のことをサンプリング周波数と呼びこの演習では sr や， Fs といった変数名で表します．下図は，サンプリング周波数 sr [Hz] でサンプリングした N 点のデータ x と scipy.fft.fft() を用いてフーリエ変換したデータ列 X との関係を模式的に表した図です．サンプリングの間隔 $\Delta = 1/F_s$ となるので，時間信号 x の i 番目の要素 x[i] は $i\Delta$ [sec] の信号を表します．一方，周波数間隔を $\tilde\Delta$ とすると，$\frac{N}{2} \tilde\Delta$ がナイキスト周波数 $F_s/2$ 対応するので， $\tilde\Delta = F_s/N$ [Hz] となり，周波数 X の j 番目の要素 X[j] は， $j \tilde\Delta$ [Hz] の周波数信号にあたることがわかります．