フーリエ級数展開

参考文献及びWEB

参考文献

(1)「ディジタル信号処理技術」玉井徳みち、長島厚、藤田泰弘、若井修造著日経BP社
(2)「ディジタル信号処理の基礎」三上直樹著 CQ出版
(3)「C言語によるディジタル信号処理入門」三上直樹著 CQ出版
(4)「アナログ&ディジタルフィルタ入門」小野浩司著日刊工業
(5)「フーリエの冒険」ヒッポファミリークラブ

フーリエ級数展開

これから説明するフーリエ級数展開の式をまず示す。

f( t+T) = f(t)

$f(t) = a_0+ \sum_{n=1}^{\infty}(a_n \cos n\omega t + b_n \sin n\omega t)$

ここで、

$a_{0} = \frac{1}{T}\int_0^{T} f(t) dt$

$a_{n} = \frac{2}{T}\int_0^{T} f(t)\cos n\omega t dt$ $n=1,2,\cdots$

$b_{n} = \frac{2}{T}\int_0^{T} f(t)\sin n\omega t dt$ $n=1,2,\cdots$

奇関数

-f(t) = f(-t) 奇関数

a_0 = a_n = 0
$b_{n} = \frac{4}{T}\int_0^{\frac{T}{2}} f(t)\sin n\omega t dt$ $n=1,2,\cdots$

偶関数

f(t) = f(-t) 偶関数

$a_0 = \frac{2}{T}\int_0^{\frac{T}{2}} f(t) dt$
$a_{1} = \frac{4}{T}\int_0^{\frac{T}{2}} f(t)\cos n\omega t dt$ $n=1,2,\cdots$

野菜ジュースを考える

フーリエ級数展開を説明する前に、
突然ですが、例として「野菜ジュース」を考えてみます。
いま、トマト、セロリ、みかんがそれぞれ次図のように配合された
野菜ジュースがあるとします。

トマト	50%
セロリ	30%
みかん	20%
計	100%

■時間軸上

    |
50g |-------------------------- <---(トマト、セロリ、みかん入りジュース)
    |
0g  |------------------------- 時間

このように野菜ジュールはどの時間からみても容量50gのジュースで
時間軸上では、トマトとセロリとみかんがどのくらいの割りあいで入っているかを
調べることはできません。

■成分軸上

                                        *
                                *       *
  |                             *       *
  |    *                        *       *
  |    *                        *       *
  -----|------|-------|---------|-------|---  成分
     みかん  ぶどう  きゃべつ  セロリ  トマト

もし上のような成分軸に世界を変えられるならば、このように各野菜の成分が分かります。

同じように、波に対して上記のような思考で考えてみます。

波の基本

結論から言いますと波というのは、
基本周波数の整数倍の周波数の波が足し合わさってできている
と考えることができます。

例として、のこぎり波で説明します。

のこぎり波は、以下の式

$f(x) = \sum_{n=1}^{\infty}(\frac{2}{n} \cdot (-1)^{(n+1)} \sin n x)$

で表すことができます。

このときを

n=1の場合
n=1～10までを足し合わせた場合
n=1～100までを足し合わせた場合

をグラフ化してみます。

図のようにに波を足し合わせればするほど三角波に近づいていきます。

$\sin$ 波と $\cos$ 波

円と正弦波のアプリ

さて、波には $\sin$ 波と $\cos$ 波があります。
円を一周したときの
yの値をプロットすると $\sin$ 波を
xの値をプロットすると $\cos$ 波を表すことができます。

$\sin$ 波も $\cos$ 波も、円を一周すると(つまり360度になると)にもとの位置へ戻ります。

正弦波を説明する

さて、正弦波は一周 360 度( $2\pi$ ラジアン)すると元に戻ります。
$\sin(2\pi) = 0$
となるわけですね。

ここで横軸を時間として、時間が経過したときの正弦波を式で表してみます。
一周した時の時間を周期&mimetex(T)で表すと、

$\sin(\frac{2\pi}{T}\cdot t)$ (1-1)

と表すことができます。
つまり t=T[s] になると一周して元の位置に戻るわけですね。

さて、この一周期の逆数は
1秒間に何回一周したかを表す
ことが分かります。

これを周波数 $f=\frac{1}{T}$ [Hz](1-2)といいます。

つまり
$y(t) = \sin(\frac{2\pi}{T} \cdot t) = \sin(2\pi f t)$ (1-3)

です。

$\cos$ 波(余弦波)というのは正弦波が $\frac{\pi}{2}$ ずれた波形です。
確かめてください。
余弦波が $\frac{\pi}{2}$ ずれた波形は正弦波です。
確かめてください。

波の性質と量

先ほど、円を回転させて
ある位置の
縦軸をプロットすると $\sin$ 波に
横軸をプロットすると $\cos$ 波が
書けました。

円の半径が各々の波形の振幅になります。

そして、波の位置は、
$\frac{2\pi}{T} \cdot t$ (1-4)
または
$2 \pi f \cdot t$ (1-5)
から求められます。

距離を求めるときには、速度に時間をかけましたよね。
つまり速度というのは、その物体のかけっこの性質を表していました。
かけっこが速い、遅いのを速度で定義するとみんなが納得できるわけです。

同じように波の性質は(1-4),(1-5)で定義できるわけです。
この
$2 \pi f = \omega$ (1-6)

と定義した $\omega$ を角速度と定義してます。

つまり、波の性質は角速度の違いで表現できるのです。

$y= \sin(\omega t)$ (1-7)

$\omega$ が違うと波の性質が違うというわけですね。

性質というのは味覚で言えば「甘い」「しょっぱい」「辛い」と同じようなものと考えてかまいません。
しかしながら甘いものとしょっぱいものを容器に入れたとして
甘いものの量が多ければ全体の味覚は甘く感じるし、しょっぱいものを入れればしょっぱく感じる。

それを表すのが先ほどの円の半径

$y=A \sin(\omega t)$ (1-8)

のAつまり振幅が量を表しているのです。

ある任意の波は基本周波数の整数倍で表される

ある任意の波は、性質を現す角速度 $\omega$ の整数倍の集まりで表現できます。
その性質がない波形は、振幅0として消えます。
ですので、各々の波の振幅(量)を求めることで、以下のような形

$f(t) = a_{0}cos(0) + b_{0}sin(0) + a_{1}\cos(\omega t)+b_{1}\sin(\omega t) + a_{2}\cos(2\omega t)+b_{2}\sin(2\omega t) + \cdots + a_{n}\cos(n \omega t)+b_{n}\sin(n \omega t)$
$= a_{0} + \sum^{\infty}_{n=1}(a_{n}\cos(n \omega t)+b_{n}\sin(n \omega t) )$

で複雑な波を表すことができます。

では複雑な波を構成する各々の波の振幅はどう求めるのでしょうか？
次からそれを説明します。

フーリエ級数展開

を求める

次に複雑な波から一つ一つの波の振幅(量)を取り出す操作をやってみます。
つまり複雑な波を個々の波に分ける作業です。

まず a_0 を考えてみる。
a_0 は波が0からどの位上にあるか　
つまりバイアス電圧のようなものと考えてもらってかまわない。

その他の波はバイアスを中心にして上下に揺れている。つまり1周期分の面積は0になる。

つまり複雑な波の1周期分の面積を求めれば図のように赤の面積部位だけが残る。

後は、それを1周期で割れば a_0 を求めることができます。
つまり

$a_0 = \frac{1}{T}\int^{T}_{0} f(t) dt$ (2-1)

となります。
この式は、
「a0は複雑な波f(t)の一周期(0からT)までの面積を求めてTで割る」
という意味です。

を求める

波は上下に揺れているので同じ周波数の波を掛けてやればマイナスの部分がプラスになる。

図の青線のように。青線をよく眺めれば
青線が0の基準線との間に作る面積は
0～T/2までの長方形の面積に等しくなることが波形の対称性より明らかである。

つまりanを求めたいならば求めたい正弦波と同じ周波数の正弦波を掛けて
マイナス成分を消しそれの１周期分の面積を求めると

横 0～T/2の長さ,
縦振幅anの長さ

を掛け合わせた長方形となる。
求めたい周波数の信号以外は、マイナスとプラスの成分が残り打ち消しあって面積0となる。

つまり

$a_n = \frac{2}{T}\int^{T}_{0} f(t)\cdot cos(n\omega t) dt$ $n=1,2,\cdots$ (2-2)
$b_n = \frac{2}{T}\int^{T}_{0} f(t)\cdot sin(n\omega t) dt$ $n=1,2,\cdots$ (2-3)

となる。
この式の意味は、

a_n は、複雑な波f(t)に求めたい波cos(nwt)をかけて、0からTまでの面積を求めてT/2で割る。
b_n は、複雑な波f(t)に求めたい波sin(nwt)をかけて、0からTまでの面積を求めてT/2で割る。

という意味です。

つまり、複雑な波を構成する各々の波の量を求めるわけです。

まとめると

(1)複雑な波の表現

$f(t) = a_{0} + \sum^{\infty}_{n=1}(a_{n}\cos(n \omega t)+b_{n}\sin(n \omega t) )$

(2)複雑な波の各々の波の成分の量(振幅)を求める式

$a_0 = \frac{1}{T}\int^{T}_{0} f(t) dt$ (2-1)

$a_n = \frac{2}{T}\int^{T}_{0} f(t)\cdot cos(n\omega t) dt$ $n=1,2,\cdots$ (2-2)
$b_n = \frac{2}{T}\int^{T}_{0} f(t)\cdot sin(n\omega t) dt$ $n=1,2,\cdots$ (2-3)

もう一度　のこぎり波を考える

のこぎり波は、
$f(x) = \sum_{n=1}^{\infty}(\frac{2}{n} \cdot (-1)^{(n+1)} \sin n x)$
のように表される。

この式を求める。

フーリエ展開

■問い
f(x) = x を $-\pi \lt x \lt \pi$ 　なる区間でフーリエ展開する。

■解答

$a_{0} = \frac{1}{T}\int^{T}_{0} f(x) dx$
$= \frac{1}{2\pi}\int^{\pi}_{-\pi} x dx$
$= \frac{1}{2\pi}\Big[\frac{x^2}{2}\Big]^{\pi}_{-\pi}$
= 0

は

$a_n = \frac{2}{T}\int^{T}_{0} f(x)\cdot cos(n x) dx$
$= \frac{1}{\pi}\int^{\pi}_{-\pi} x\cdot cos(n x) dx$
$= \frac{1}{n\pi}\int^{\pi}_{-\pi} x\cdot (sin(n x))' dx$
$= \frac{1}{n\pi}\Big{ \[ x\sin(n x) \]^{\pi}_{-\pi} -\int^{\pi}_{-\pi} \sin(nx) dx \Big}$
$= \frac{1}{n\pi}\[ \frac{1}{n}\cos(n x) \]^{\pi}_{-\pi}$
= 0

は

$b_n = \frac{2}{T}\int^{T}_{0} f(x)\cdot \sin(n x) dx$
$= \frac{1}{\pi}\int^{\pi}_{-\pi} x\cdot \sin(n x) dx$
$= -\frac{1}{n\pi}\int^{\pi}_{-\pi} x\cdot (\cos(n x))' dx$
$= -\frac{1}{n\pi}\Big{ \[ x\cos(n x) \]^{\pi}_{-\pi} -\int^{\pi}_{-\pi} \cos(nx) dx \Big}$
$= -\frac{1}{n\pi}\[ x\cos(n x) \]^{\pi}_{-\pi} -\[\frac{1}{n}\sin(nx)\]^{\pi}_{-\pi}$
$= -\frac{1}{n\pi}\[ \pi\cos(n \pi) - (-\pi)\cos(n(-\pi)) \]$
$= -\frac{1}{n\pi}\[ 2\pi\cos(n \pi) \]$
$= -\frac{2}{n}\[ \cos(n \pi) \]$
$= \frac{2}{n}(-1)^{(n+1)}$

従って複雑な波は

$f(t) = a_{0} + \sum^{\infty}_{n=1}(a_{n}\cos(n x)+b_{n}\sin(n x) )$

$f(x) = \sum_{n=1}^{\infty}( \frac{2}{n}(-1)^{(n+1)} \sin n x)$

となる。

Matlabスクリプト

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 三角波形を表示
%             programming by embedded.samurai
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
echo off
clear all
close all

%サンプリング周期 1[ms]
dx=1e-3;
% xの範囲は -pi/2から0.001単位で3piまで
x=[-2*pi:dx:2*pi]; 
xx = x/pi;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% n=1のとき
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1);

y1=0;
for m=1:1:1
 y=(2/m)*(-1)^(m+1)*sin(m.*x);
 y1=y1+y;
end

y1
plot(xx,y1/pi,'r-','linewidth',3);
hold on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% n=10のとき
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y2=0;
for m=1:1:10
 y=(2/m)*(-1)^(m+1)*sin(m.*x);
 y2=y2+y;
end

plot(xx,y2/pi,'g-','linewidth',3);
hold on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% n=100のとき
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y3=0;
for m=1:1:100
 y=(2/m)*(-1)^(m+1)*sin(m.*x);
 y3=y3+y;
end

plot(xx,y3/pi,'m-','linewidth',3);
hold on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 理想パルス
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
data=[];

for n=x
 if ( n >= -2*pi) && ( n < -pi)
   y=n  +2*pi;
   data=[data,y];
 
 elseif (n >= -pi) && (n < pi)
   y=n;
   data=[data,y];
 elseif (n >= pi) && (n <= 2*pi) 
   y=n -2*pi;
   data=[data,y];
 
 end
end

plot(xx,data/pi,'b-','linewidth',3);
hold on;

% xラベルとその文字の大きさ、線の太さの設定
str1={'-3p','-2p','-p','0','p','2p','3p'}
set(gca,'FontName','symbol','xtick',[-3:1:3],'xticklabel',str1)

str2={'-3p','-2p','-p','0','p','2p','3p'}
set(gca,'FontName','symbol','ytick',[-3:1:3],'yticklabel',str2)

set(gca,'LineWidth',2,'FontSize',18)

% x-y範囲
axis([-Inf Inf -Inf Inf]);
grid on

% xラベル、yラベル
xlabel('Positon t','Fontsize',20,'FontName','Century');
ylabel('pulsu wave y(t)','Fontsize',20,'FontName','Century')

% キャプション
h=legend('n=1',...,
'n=10',...,
'n=100',...,
'ideal triangle wave',...,
4);
set(h,'FontSize',15,'FontName','Century');

print -djpeg sankakuwave_t.jpg

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2);

%サンプリング周波数 1kHz
df=1/dx
%xの数 7854
[gyo retu]=size(x)

%基本周波数 0.1273[Hz]
bf=df/retu

%FFTの結果は複素数である
fftdata1 = fft(y1)';
fftdata2 = fft(y2)';
fftdata3 = fft(y3)';
fftdata4 = fft(data)';

%FFTの絶対値を求める
power1=abs(fftdata1).^2;
power2=abs(fftdata2).^2;
power3=abs(fftdata3).^2;
power4=abs(fftdata4).^2;

%離散周波数を作る
freq = [0:retu-1]'*bf

plot(freq,power1,'r-','linewidth',3);
hold on
plot(freq,power2,'g-','linewidth',3);
hold on
plot(freq,power3,'m+--','linewidth',2);
hold on
plot(freq,power4,'bo-','linewidth',1);
hold on

axis([0 1  0 Inf]);

set(gca,'LineWidth',2,...,
'FontSize',15);

xlabel('Frequency[Hz]','Fontsize',20,'FontName','Century');
ylabel('Normalize','Fontsize',20,'FontName','Century');

% キャプション
h=legend('n=1',...,
'n=10',...,
'n=100',...,
'ideal triangle wave',...,
1);
set(h,'FontSize',20,'FontName','Century');
grid on

print -djpeg sankakuwave_f.jpg