第2章

フーリエ変換による周波数分析

2. 1 実験目的
本実験では，フーリエ変換という基本的な周波数分析の方法について MATLAB のシミュレーション実験を通して理解することを目的としている。

2. 2 実験原理と内容
2. 2. 1 フーリエ級数展開
任意の周期信号はいろいろな周波数の三角関数の無限の和として表すことができる。これがフランスの有名な数学者 Fourier(1768～1830)が考え出した，フーリエ級数展開である。
まず周期 T （秒）で繰り返される信号を x(t ) とする。フーリエ級数の考え方は，信号 x(t ) は T 秒で 1周期。即ち周波数 1 / T （Hz）の成分と，その整数倍の周波数の成分[周波数が 2 / T ， 3 / T ， 4 / T !
（Hz）]の無限和として表される，となる。これを式で表すと

2π t
4π t
6π t+ a2 cos+ a3 cos+!!
TTT
2π mt
2π t
4π t!+ am cos+ b1 sin+ b2 sin
TTT
6π t
2π mt+ b3 sin+!!+ bm sin
TT∞∞
2π mt
2π mt
= a0 + ∑ am cos+ ∑ bm sinTT
m=1
m=1x(t) = a0 + a1 cos（1）と な る 。 式 （ 1 ） に お い て ， a0 は 信 号 x(t ) に お け る 直 流 成 分 の 大 き さ を 表 す 。 cos 2p t / T ，sin 2p t / T を基本波成分と呼ぶ。これは周期 T （秒）即ち周波数 1 / T （Hz）の cos 波（余弦波）とsin 波（正弦波）の成分である。 cos 2p mt / T ， sin 2p mt / T を高調波成分とよばれ，周期が T / m（秒），言い換えると，周波数が m / T （Hz）の cos 波， sin 波成分である。係数 a m ， bm はフーリエ係数とよばれ，それぞれ cos ， sin 波成分の振幅を表す。
フーリエ係数を原信号 x(t ) から求める場合には次のようになる。

2-1

直流成分の振幅 a 0 =cos 成分 a m =1 Tx(t ) dt
T ò0

また， 周波数

2p mt
2 Tx(t ) cosdt
ò
T 0Tm（Hz）の成分の振幅T

sin 成分 bm =2p mt
2 Tx(t ) sindtò
T 0T

例題1： 三角関数の角度，ラジアン及び周期波形

sin と cos は 360 ! （度）を1周期とする三角関数である。度のかわりにラジアン（rad）を用いると 1! = p / 180 (rad) であるから，下表に対応関係を求められる。但し， p = 3.14159 !である。

0 ・・・ 90 ・・・ 180 ・・・ 270 ・・・ 360

度p ・・・ 3p / 2 ・・・ 2p
ここで， m = 1, 2 の場合，時間 t に対応する sin 2p mt / T と cos 2p mt / T の値，及び波形は次のよう

ラジアン

0 ・・・p / 2 ・・・に求められる。

m = 1のとき
t0 ・・・ T / 4 ・・・ T / 2 ・・・ 3T / 4 ・・・ T2p mt / T

0 ・・・ p / 2 ・・・p・・・ 3p / 2 ・・・ 2p

sin 2p mt / T

0 ・・・

1・・・0・・・cos 2p mt / T

1 ・・・0・・・- 1 ・・・- 1 ・・・ 00・・・ 1m = 2 のとき
0 ・ T / 8 ・T / 4 ・・・ ・・・ ・・・・・・・

T2p mt / T

0 ・ p /2 ・p・・・ ・・・ ・・・・・・・

4p

sin 2p mt / T

0 ・ 1・0・・・ ・・・ ・・・・・・・

0cos 2p mt / T

1 ・ 0・- 1・・・ ・・・ ・・・・・・・

1t

上記の例題から sin 2p t / T （ cos 2p t / T ）の周期は T （秒）で， sin 4p t / T （ cos 4p t / T ）の周期は

T / 2 （秒）であることが直感的に分かる。同様に sin 2p mt / T （ cos 2p mt / T ）の周期は T / m （秒）2-2で，周波数は m / T （Hz）となる。

2．2．2 振幅，パワー，位相及びそれらのスペクトル同じ周波数の cos と sin は位相に 90 度のずれがあることを利用して，式（1）を書き換えると，次式が得られる。

⎛ 2π mt⎞x(t) = a0 + ∑ Am cos ⎜− θm ⎟⎝ T⎠
m=1∞（2）

こ こ で ， Am = a m2 + bm2 を 振 幅 と 言 い ， そ の 2 乗 ( | Am |2 ) を 振 幅 の パ ワ ー と 言 う 。 ま たq m = tan -1 (bm / a m )を位相と呼ばれている。
周波数の関数としてのフーリエ係数，すなわち周波数の変化に対するフーリエ係数の変化をスペクトルとよんでいる。周波数対振幅は振幅スペクトルで，周波数対パワーはパワースペクトルで，周波数対位相は位相スペクトルである。
各周波数成分を纏めた振幅，パワー，位相及びスペクトルは同じように複素フーリエ級数からも求められる。

例題 2： 周期信号のスペクトル
合成信号の周期は T = 100ms であり，次式のように表す。

2pt
2 ´ 5pt p
2 ´ 7pt p- 0) + 10 × cos(- ) + 5 × cos(- )TT
2T
42 ´ 12pt p
2 ´ 15pt+ 2 × cos(- ) + 1.5 × cos(- 0)T
6Tx(t ) = 8 × cos((3)

周期 T （秒）の信号に対し，振幅スペクトル，パワースペクトル及び位相スペクトルを求めたもの図(a)～図（ｄ）に示せ。
合成信号の周期は T = 100ms であるため，その周波数の間隔は 1 / T = 10 Hz である。番号m = 1 から m = 20 までに設定すると，周波数の範囲は 0～200Hz となる。ここで，式（4.2）を展開した式（4.3）と比較すれば，番号 m = 1 のとき，すなわち周波数 10Hz であるときに，式（4.3）の第一項から
振幅スペクトル | Am | =8 で，パワースペクトル | Am |2 =64 及び位相スペクトル q m = 0 であることが直感的に分かる。同様に m = 5，7，12，15，周波数が 50，70，120，150 Hz の各スペクトルを求め（この部分

2-3

は課題として解いてもらう。），図(b)～図（ｄ）に示すようなスペクトルを描くことができる。
例題 2 においては，信号の周期性があるため，図(b)～図（ｄ）ように，求められたスペクトルは連続でなく，ある周波数間隔で出現するので，これは線スペクトルと呼ばれている。信号が周期性のない場合（周期 T ® ¥ と看做される）の場合には，フーリエ変換（積分）で求められたスペクトルが線スペクトルでなく，連続的なスペクトルとなる。

2．2．3 離散フーリエ変換(DFT)と高速フーリエ変換(FFT)
フーリエ変換は時間領域の信号を周波数領域におけるスペクトル(spectrum)として観察するための手段としてよく用いられる。離散フーリエ変換は有限長の離散時間信号から，その周波数によって表現される数列に変換する一つの計算方法である。離散時間フーリエ変換は次式のように定義される。

N -1

X ( k ) = å x ( n )e

-j

2pkn
N,k = 0, 1, ! , N - 1

(4)

n =0

ここで，x(n)は時間領域の信号である。N は離散時間信号のサンプル数であり，周波数 0～ 2p を等分割するための数である。また，インデックス k は周波数を表す整数，フーリエ変換された信号の絶対値|X(k)|は信号 x(n)の振幅スペクトル(magnitude spectrum)，また角度 ÐX (k ) は位相スペクトル(phase spectrum)である。
離散時間フーリエの逆変換は次式のようになる。

2-4

2px ( n) =jkn
1 N -1X ( k )e N ,åN k =0n = 0, 1, ! , N - 1(5)

上述した離散時間フーリエ変換とその逆変換は，それぞれ周波数領域における N 個の複素数，時間領域における N 個の数列を作成するための計算方法と考えてもよい。これらの計算法は Nが大きくなると演算量の問題が生じてくる。そのため，高速フーリエ変換(FFT)が考案されている。
Matlab においては，高速フーリエ変換の計算は関数 fft()で行い，その逆変換の計算は関数 ifft()で行う。
例題 3：周波数 f1=697[Hz]，f2=1209[Hz]の正弦波信号 x(t ) = sin( 2p f 1t ) + sin( 2pf 2 t ) をサンプリ

ング周波数 fs=8192[Hz]で標本化し，得られた離散信号に対して高速フーリエ変換を行いその振幅特性を示してみよう。また，そのフーリエ変換された信号に対して逆高速フーリエ変換してみよう。
本例題の Matlab 計算プログラムを次表のようになる。
% Investigation of the fft() and ifft()% Filename: Samp6_3.m% Jianting Cao, Feb. 28, 2003
clear all
f1=697;f2=1209; % 周波数の定義
fs=8192;% サンプリング周波数の定義

n=[0:1/fs:0.1]; % サンプル数の定義
x=sin(2*pi*f1*n)+sin(2*pi*f2*n);%合成信号の生成
X=fft(x);m=abs(X);% 信号のフーリエ変換の計算% 振幅特性の計算

xx=real(ifft(X));% 信号の逆フーリエ変換の計算
t=linspace(0,0.01,82);% 元信号のプロット
subplot(3,1,1);plot(t,x(1:length(t)));grid;
xlabel('Time [sec]','FontSize',8);ylabel('Amplitude','FontSize',8);title('Original Signal, FFT signal and IFFT signal','FontSize',10)
subplot(3,1,2); %フーリエ変換した信号のプロット

2-5

f=(0:length(X)-1)*fs/length(X);plot(f,m);axis([0 1500 0 450]);grid;
xlabel('Frequency [Hz]','FontSize',8);ylabel('Magnitude','FontSize',8);
subplot(3,1,3); %逆フーリエ変換した信号のプロット
plot(t,xx(1:length(t)),'r-');grid;
xlabel('Time [sec]','FontSize',8);ylabel('Amplitude','FontSize',8);

上記のプログラムを実行すると，次の結果が得られる。

上段の波形は，サンプリング周波数 fs=8192 [Hz]で標本化し得られたものである。中段の波形は，上段の離散信号に対して求められた振幅特性を示している。この波形に現れている二つのピークで示される周波数はそれぞれ f1，f2 と一致していることが分かる。更に下段の波形は，その逆フーリエ変換されたものであり，上段の離散信号と一致していることが示されている。

2-6

2．3 実験課題
１．例題1において， m = 3 のとき，下表に「？」マックが付けられた所の値を求めよ。また，時間範囲

t = [0 T ]の正弦波と余弦波をそれぞれで描いてください。

0 ・ T / 12 ・T / 6 ・・・ ・・・ ・・・・・・・ Tt
2p mt / T? ・?

sin 2p mt / T? ・cos 2p mt / T?

・?

・・・ ・・・ ・・・・・・・ ?

?

・?

・・・ ・・・ ・・・・・・・ ?

?

・ ?

・・・ ・・・ ・・・・・・・・?

２．例題２の続きとして，周波数が10，50，70，120，150 Hzであるとき（番号 m = 1，5，7，12，15）の振幅スペクトル，パワースペクトル及び位相スペクトルをそれぞれ求め，下表に入れてください。
周波数 (Hz)10振幅スペクトル | Am |701201508パワースペクトル | Am |264q m (rad)0位相スペクトル

50p /2

3. プッシュ回線の電話機で発信の際にダイヤルボタンを押す毎に発信される信号音が伴っている。これは Dual Tone Multi Frequency(DTMF)と呼ばれ，一つの標準となっており，1 つのボタンにつき 2 種類の周波数が割り当てられる。例えば，電話機番号のボタン「1」の音は周
波数 f1=697[Hz]，f2=1209[Hz]の正弦波信号 x(t ) = sin( 2p f 1t ) + sin( 2pf 2 t ) を用いて発生させる。各番号のボタンと周波数の対応関係を次表のように示す。

2-71209 Hz 1336 Hz 1477 Hz１）

697 Hz123770 Hz

456852 Hz

789941 Hz*0#

この表を利用して，各自の学籍番号の最後 1 桁の数字による合成正弦波を作成せよ。但しサンプリング周波数は fs=8192[Hz]とする。

２）

合成された正弦波を音声関数 sound()で聞いてみて，自分が持っている携帯電話でその数字を押すときの音と一致するかどうかを比べてみよ。

３）

次に，その作成した正弦波に対して振幅特性と位相特性[関数 unwrap(angle(fft(・)))を利用すること]を求めよ。

2-8