「諸概念の迷宮(Things got frantic)」用語集

本編で頻繁に使うロジックと関連用語のまとめ。

【高校数学】フーリエ解析(Fourier analysis)による波形合成について

涌井良幸「高校生からわかるフーリエ解析2019年)」第1章の内容がコンピューター言語による検証向けだったので、早速統計言語Rでのプログラミングにチャレンジしてみました。

f:id:ochimusha01:20190815124456j:plain

 そもそも「フーリエ解析Fourier analysis)」は、全ての波形について(Sin波やCos波の様な)正弦波に適正な数字を掛けて足し合わせたら再現可能と考えます。 

正弦波(sinusoidal wave) - Wikipedia

正弦関数として観測可能な周期的変化を示す波動。その波形は正弦曲線(せいげんきょくせん、sine curve)もしくはシヌソイド (Sinusoid) と呼ばれ、数学、信号処理、電気工学およびその他の分野において重要な働きをする。

自然現象における正弦波

単一の周波数成分のみを持つ波動であり、厳密な意味では自然界には存在しないが、一般の物理学や電磁気学、音響学などにおいて観測されるべき波動の振幅が、付随される雑音に比べて十分に大きい場合、これを正弦波と見なすことが多い。 この広義の意味での正弦波は自然界でも海の波、音波、光波などで発生する。また、日々の平均気温を年間を通してプロットした際などにも荒いシヌソイドパターンが現れる。

商用電源など発電機から得られる交流電圧の波形も一般に正弦波形をとる。

観測対象としての波動

バネによって吊り下げられた重りの振動は、平衡点まわりでは正弦波として近似できる。

f:id:ochimusha01:20190815193046g:plain

基本形…固定された観測位置における正弦波は次のような関数として記述することができる。

y=A*sir(ωt-φ)

  • ただしtは時刻(Times)、A は振幅(Amplitude、波の中心からの最大偏差)、ω(オメガ)は角周波数(Angular frequency、角振動数、円振動数、国際単位系では単位はラジアン毎秒)、−φ(ファイ)は初期位相(IP=Initial phase、t = 0 における位相)という。
  • −φは位相シフトとも関係がある。例えば、初期位相−φが負の値であれば、波形全体が未来の時間へシフトされる、すなわち波の到達が遅れる。シフトされる時間は、φ / ω である。

一般形…これにさらに波動の発生源からの距離xや波数k 、直流成分(振幅の中心となる値) D などを含めると以下の様な関数となる。
y=A*sir(kx-ωt-φ)+D

  • この関数の形で波形を記述できるものを正弦波と総称する。
  • 波数は角周波数と「k=ω/c=2πf/c=2π/λ」のような関係にある。ここで、λ は波長(独:Wellenlänge、英:wavelength、空間を伝わる波(波動)の持つ周期的な長さ)、f は周波数(frequency、工学、特に電気工学・電波工学や音響工学などにおいて、電気振動(電磁波や振動電流)などの現象が、単位時間(ヘルツの場合は1秒)当たりに繰り返される回数)、c は位相速度(Phase velocity、単位ラジアン)である。
  • 1次元の正弦波となるため、時刻 t における位置 x での波の振幅が導かれる。これは例えば、ワイヤーに沿った波の値と考えることが出来る。

コサイン波形(余弦)もシヌソイドと言われる。これは、正弦波が後方にシフトされたもので波形が同一だからである。

cos(x)=sin(x+1/2)

なお、正弦関数は波動方程式ヘルムホルツ方程式を満たす最も基本的な関数でもある。

 

フーリエ解析Fourier analysis

1822年、フランス人数学者のジョゼフ・フーリエは、周期的な波動をさまざまな(基本周波数の整数倍の)周波数の正弦波の重ね合わせとして表す方法を発見した。この方法はフーリエ級数またはフーリエ級数展開と呼ばれ、信号処理におけるもっとも基礎的な手法の一つである。

また、単一のパルス波や人の声による不規則な音波といった周期的でない波形も、連続的に変化する異なった周波数の波を重ね合わせて表すことができる。このような一般的で複雑な波を様々な周波数の正弦波に分解して解析する手法はフーリエ変換と呼ばれている。

音楽への応用

1950年代、正弦波音がオルガンの音に似ているということも好都合であり、電子音楽の黎明期に愛用された。この時流に沿う形で、フランスの作曲家アンリ・プッスールは「正弦波曲線が、楽曲の理想的な形式」と定義して話題となった。この理論を杓子定規に応用した作品に、篠原眞の「タンダンス(Tendance=傾向)」がある。
また、1980年代には、正弦波に対し変調を掛けることによって波形を生成するFM音源方式の楽器が発売され、少ないパラメータから算出される多彩な波形によって、一時代を築いた。一部の製品では、正弦波を単体で扱う機能や、予め変換された波形を再生することを支援する機能も実装された。

 

①その正弦波を特徴付ける主要パラメーター(parameter、媒介変数、補助変数、母数、引数、設定値などの意味を持つ英単語。 IT分野ではそのソフトウェアやシステムの挙動に影響を与える、外部から投入されるデータなどのことを指す)は「振幅Amplitude、基本円ではx<=|1|)」と「波長wavelength、基本円では2π)」/「位相Phase、単位はラジアン/秒)」/「周波数frequency、基本円では位相/波長)」の兼ね合いとなる。

統計言語Rによる実習

#グラフの色の決定
Black<-rgb(0,0,0)
Red<-rgb(1,0,0)
Magenta<-rgb(1,0,1)
Blue<-rgb(0,0,1)
Green<-rgb(0,1,0)
Cyan<-rgb(0,1,1)
Yellow<-rgb(1,1,0)
Gray<-"#888888"

#タイトル定義
Main_title<-c("Amplitude & Frequency")
x_title<-c("X")
y_title<-c("Y")

#グラフのスケール決定
gs_x<-c(-pi,pi)
gs_y<-c(-3,3)
#y=sin(x)
f1_name<-c("y=sin(x)") 
glaph_y1<-function(x) {sin(x)}

#y=2*sin(x)
f2_name<-c("y=2*sin(x)") 
glaph_y2<-function(x) {2*sin(x)}
#y=3*sin(x)
f3_name<-c("y=3*sin(x)") 
glaph_y3<-function(x) {3*sin(x)}
#y=sin(x*2)
f4_name<-c("sin(x*2)") 
glaph_y4<-function(x) {sin(x*2)}

#y=sin(x*3)
f5_name<-c("y=sin(x*3)") 
glaph_y5<-function(x) {sin(x*3)}
#y=sin(x*4)
f6_name<-c("y=sin(x*4)") 
glaph_y6<-function(x) {sin(x*4)}
#グラフ描画
plot(glaph_y1,xlim=gs_x,ylim=gs_y,type="l",col=Red, main=Main_title,xlab=x_title,ylab=y_title)
par(new=T)#上書き指定
plot(glaph_y2,xlim=gs_x,ylim=gs_y,type="l",col=Magenta, main="",xlab="",ylab="")
par(new=T)#上書き指定
plot(glaph_y3,xlim=gs_x,ylim=gs_y,type="l",col=Blue, main="",xlab="",ylab="")
par(new=T)#上書き指定
plot(glaph_y4,xlim=gs_x,ylim=gs_y,type="l",col=Green, main="",xlab="",ylab="")
par(new=T)#上書き指定
plot(glaph_y5,xlim=gs_x,ylim=gs_y,type="l",col=Cyan, main="",xlab="",ylab="")
par(new=T)#上書き指定
plot(glaph_y6,xlim=gs_x,ylim=gs_y,type="l",col=Yellow, main="",xlab="",ylab="")
#par(new=T)#上書き指定
#plot(f7,xlim=gs_x,ylim=gs_y,type="l",col=Blue, main="",xlab="",ylab="")

#基準線
abline(h=0)
abline(v=0)
#凡例描画
legend("topleft", legend=c(f1_name,f2_name,f3_name,f4_name,f5_name,f6_name),lty=c(1,1,1,1,1,1),col=c(Red,Magenta,Blue,Green,Cyan,Yellow))

f:id:ochimusha01:20190815220701p:plain

②この様に定数倍した正弦波を幾つか足し合わせただけでも、随分と複雑な波形が合成可能である。

統計言語Rによる実例1「y=sin(x)+sin(x*2)+sin(x*3)」

#グラフの色の決定
Black<-rgb(0,0,0)
Red<-rgb(1,0,0)
Magenta<-rgb(1,0,1)
Blue<-rgb(0,0,1)
Green<-rgb(0,1,0)
Cyan<-rgb(0,1,1)
Yellow<-rgb(1,1,0)
Gray<-"#888888"

#タイトル定義
Main_title<-c("Mix of Sine waves ")
x_title<-c("X")
y_title<-c("Y")

#グラフのスケール決定
gs_x<-c(-pi,pi)
gs_y<-c(-3,3)
#y=sin(x)
f1_name<-c("sin(x)") 
glaph_y1<-function(x) {sin(x)}

#y=sin(x*2)
f2_name<-c("sin(x*2)") 
glaph_y2<-function(x) {sin(x*2)}
#y=sin(x*3)
f3_name<-c("sin(x*3)") 
glaph_y3<-function(x) {sin(x*3)}
#y=sin(x)+sin(x*2)+sin(x*3)
f4_name<-c("mix of these") 
glaph_y4<-function(x) {sin(x)+sin(x*2)+sin(x*3)}

 

#グラフ描画
plot(glaph_y1,xlim=gs_x,ylim=gs_y,type="l",col=Red, main=Main_title,xlab=x_title,ylab=y_title)
par(new=T)#上書き指定
plot(glaph_y2,xlim=gs_x,ylim=gs_y,type="l",col=Magenta, main="",xlab="",ylab="")
par(new=T)#上書き指定
plot(glaph_y3,xlim=gs_x,ylim=gs_y,type="l",col=Blue, main="",xlab="",ylab="")
par(new=T)#上書き指定
plot(glaph_y4,xlim=gs_x,ylim=gs_y,type="l",col=Black, main="",xlab="",ylab="")

#基準線
abline(h=0)
abline(v=0)
#凡例描画
legend("topleft", legend=c(f1_name,f2_name,f3_name,f4_name),lty=c(1,1,1,1),col=c(Red,Magenta,Blue,Black))

f:id:ochimusha01:20190815224546p:plain

統計言語Rによる実例1「y=-sin(x)+cos(x*2)+1/2*sin(x*3)」

#グラフの色の決定
Black<-rgb(0,0,0)
Red<-rgb(1,0,0)
Magenta<-rgb(1,0,1)
Blue<-rgb(0,0,1)
Green<-rgb(0,1,0)
Cyan<-rgb(0,1,1)
Yellow<-rgb(1,1,0)
Gray<-"#888888"

#タイトル定義
Main_title<-c("Mix of Sine waves")
x_title<-c("X")
y_title<-c("Y")

#グラフのスケール決定
gs_x<-c(-pi,pi)
gs_y<-c(-3,3)
#y=-sin(x)
f1_name<-c("-sin(x)") 
glaph_y1<-function(x) {-sin(x)}

#y=cos(x*2)
f2_name<-c("cos(x*2)") 
glaph_y2<-function(x) {cos(x*2)}
#y=1/2*sin(x*3)
f3_name<-c("1/2*sin(x*3)") 
glaph_y3<-function(x) {1/2*sin(x*3)}
#y=-sin(x)+cos(x*2)+1/2*sin(x*3)
f4_name<-c("mix of these") 
glaph_y4<-function(x) {-sin(x)+cos(x*2)+1/2*sin(x*3)}

#グラフ描画
plot(glaph_y1,xlim=gs_x,ylim=gs_y,type="l",col=Red, main=Main_title,xlab=x_title,ylab=y_title)
par(new=T)#上書き指定
plot(glaph_y2,xlim=gs_x,ylim=gs_y,type="l",col=Magenta, main="",xlab="",ylab="")
par(new=T)#上書き指定
plot(glaph_y3,xlim=gs_x,ylim=gs_y,type="l",col=Blue, main="",xlab="",ylab="")
par(new=T)#上書き指定
plot(glaph_y4,xlim=gs_x,ylim=gs_y,type="l",col=Black, main="",xlab="",ylab="")

#基準線
abline(h=0)
abline(v=0)
#凡例描画
legend("bottomleft", legend=c(f1_name,f2_name,f3_name,f4_name),lty=c(1,1,1,1),col=c(Red,Magenta,Blue,Black))

f:id:ochimusha01:20190815225149p:plain

統計言語Rによる実例3「y=sin(t)+cos(t)+1/2*sin(2*t)+1/2*cos(2*t)+1/3*sin(3*t)+1/3*cos(3*t)」

#グラフの色の決定
Black<-rgb(0,0,0)
Red<-rgb(1,0,0)
Magenta<-rgb(1,0,1)
Blue<-rgb(0,0,1)
Green<-rgb(0,1,0)
Cyan<-rgb(0,1,1)
Yellow<-rgb(1,1,0)
Gray<-"#888888"

#タイトル定義
Main_title<-c("Mix of Sine waves ")
x_title<-c("X")
y_title<-c("Y")

#グラフのスケール決定
gs_x<-c(-pi,pi)
gs_y<-c(-3,3)
#y=sin(x)
f1_name<-c("sin(x)") 
glaph_y1<-function(x) {sin(x)}

#y=cos(x)
f2_name<-c("cos(x)") 
glaph_y2<-function(x) {cos(x)}
#y=1/2*sin(2*x)
f3_name<-c("1/2*sin(2*x)") 
glaph_y3<-function(x) {1/2*sin(2*x)}
#y=1/2*cos(2*x)
f4_name<-c("1/2*cos(2*x)") 
glaph_y4<-function(x) {1/2*cos(2*x)}

#y=1/3*sin(3*x)
f5_name<-c("1/3*sin(3*x)") 
glaph_y5<-function(x) {1/3*sin(3*x)}
#y=1/3*cos(3*x)
f6_name<-c("1/3*cos(3*x)") 
glaph_y6<-function(x) {1/3*cos(3*x)}

#y=sin(x)+cos(x)+1/2*sin(2*x)+1/2*cos(2*x)+1/3*sin(3*x)+1/3*cos(3*x)
f7_name<-c("mix of these") 
glaph_y7<-function(x) {sin(x)+cos(x)+1/2*sin(2*x)+1/2*cos(2*x)+1/3*sin(3*x)+1/3*cos(3*x)}

#グラフ描画
plot(glaph_y1,xlim=gs_x,ylim=gs_y,type="l",col=Red, main=Main_title,xlab=x_title,ylab=y_title)
par(new=T)#上書き指定
plot(glaph_y2,xlim=gs_x,ylim=gs_y,type="l",col=Magenta, main="",xlab="",ylab="")
par(new=T)#上書き指定
plot(glaph_y3,xlim=gs_x,ylim=gs_y,type="l",col=Blue, main="",xlab="",ylab="")
par(new=T)#上書き指定
plot(glaph_y4,xlim=gs_x,ylim=gs_y,type="l",col=Green, main="",xlab="",ylab="")
par(new=T)#上書き指定
plot(glaph_y5,xlim=gs_x,ylim=gs_y,type="l",col=Cyan, main="",xlab="",ylab="")
par(new=T)#上書き指定
plot(glaph_y6,xlim=gs_x,ylim=gs_y,type="l",col=Yellow, main="",xlab="",ylab="")
par(new=T)#上書き指定
plot(glaph_y7,xlim=gs_x,ylim=gs_y,type="l",col=Black, main="",xlab="",ylab="")

#基準線
abline(h=0)
abline(v=0)
#凡例描画
legend("topleft", legend=c(f1_name,f2_name,f3_name,f4_name,f5_name,f6_name,f7_name),lty=c(1,1,1,1,1,1,1),col=c(Red,Magenta,Blue,Green,Cyan,Yellow,Black))

f:id:ochimusha01:20190815223523p:plain

それでは実際の合成過程に目を向けて見ましょう。

ノコギリ波鋸歯状波(きょしじょうは)、英Sawtooth wave)は非正弦波的な基本的波形の一種で、波形の見た目が鋸の歯のように見えることからそのように呼ばれる。

簡単に説明すれば、その波形は時間と共に上がっていき、急降下するということを繰り返す。もちろん、逆に徐々に下がっていって急上昇することを繰り返すのこぎり波もあり、後者を「逆のこぎり波英reverse sawtooth wave、inverse sawtooth wave)」と呼ぶ。どちらの波形であってもパラメータを調整すると同じ音に聞こえる。

  • バーチャルアナログ音源や減算方式などのほとんど全てのシンセサイザーの基本となっている。
  • ブラウン管などでの画像表示のための走査線を制御する信号の波形である(ラスタースキャン)。
  • オシロスコープも時系列波形を描く際の水平方向の電子線偏向制御にこれを使っている。

音として聞いてみると、猛々しくハッキリしていて、周波数成分としては基本周波数の偶数倍音と奇数倍音の両方が含まれている。全ての整数倍音を含んでいるため、減算方式のシンセサイザーで、他の音を合成するベースとして使うのに便利である。

フーリエ級数を使用して、無限級数の形で理想的な波形を表すことができるが、この表現で興味深いのは矩形波でも起こる「ギブズ現象Gibbs phenomenon)」である。

統計言語Rによる波形合成例

Sawtooth_wave_exp<-function(x){

#グラフのスケール決定
Graph_scale_x<-c(-pi,pi)
Graph_scale_y<-c(0,pi)
#関数と凡例の決定
switch(x,
"0"= f0<-function(x) {pi/2*x/x},
"1"= f0<-function(x) {pi/2-sin(x*2)},
"2"= f0<-function(x) {pi/2-sin(x*2)-sin(x*4)/2},
"3"= f0<-function(x) {pi/2-sin(x*2)-sin(x*4)/2-sin(x/6)/3},
"4"= f0<-function(x) {pi/2-sin(x*2)-sin(x*4)/2-sin(x*6)/3-sin(x*8)/4},
"5"= f0<-function(x) {pi/2-sin(x*2)-sin(x*4)/2-sin(x*6)/3-sin(x*8)/4-sin(x*10)/5},
"6"= f0<-function(x) {pi/2-sin(x*2)-sin(x*4)/2-sin(x*6)/3-sin(x*8)/4-sin(x*10)/5-sin(x*12)/6},
"7"= f0<-function(x) {pi/2-sin(x*2)-sin(x*4)/2-sin(x*6)/3-sin(x*8)/4-sin(x*10)/5-sin(x*12)/6-sin(x*14)/7},
"8"= f0<-function(x) {pi/2-sin(x*2)-sin(x*4)/2-sin(x*6)/3-sin(x*8)/4-sin(x*10)/5-sin(x*12)/6-sin(x*14)/7-sin(x*16)/8},
"9"= f0<-function(x) {pi/2-sin(x*2)-sin(x*4)/2-sin(x*6)/3-sin(x*8)/4-sin(x*10)/5-sin(x*12)/6-sin(x*14)/7-sin(x*16)/8-sin(x*18)/9},
"10"= f0<-function(x) {pi/2-sin(x*2)-sin(x*4)/2-sin(x*6)/3-sin(x*8)/4-sin(x*10)/5-sin(x*12)/6-sin(x*14)/7-sin(x*16)/8-sin(x*18)/9-sin(x*20)/10},
"11"= f0<-function(x) {pi/2-sin(x*2)-sin(x*4)/2-sin(x*6)/3-sin(x*8)/4-sin(x*10)/5-sin(x*12)/6-sin(x*14)/7-sin(x*16)/8-sin(x*18)/9-sin(x*20)/10-sin(x*22)/11-sin(x*24)/12},
"12"= f0<-function(x) {pi/2-sin(x*2)-sin(x*4)/2-sin(x*6)/3-sin(x*8)/4-sin(x*10)/5-sin(x*12)/6-sin(x*14)/7-sin(x*16)/8-sin(x*18)/9-sin(x*20)/10-sin(x*22)/11-sin(x*24)/12-sin(x*26)/13-sin(x*28)/14},
"13"= f0<-function(x) {pi/2-sin(x*2)-sin(x*4)/2-sin(x*6)/3-sin(x*8)/4-sin(x*10)/5-sin(x*12)/6-sin(x*14)/7-sin(x*16)/8-sin(x*18)/9-sin(x*20)/10-sin(x*22)/11-sin(x*24)/12-sin(x*26)/13-sin(x*28)/14-sin(x*30)/15-sin(x*32)/16},
"14"= f0<-function(x) {pi/2-sin(x*2)-sin(x*4)/2-sin(x*6)/3-sin(x*8)/4-sin(x*10)/5-sin(x*12)/6-sin(x*14)/7-sin(x*16)/8-sin(x*18)/9-sin(x*20)/10-sin(x*22)/11-sin(x*24)/12-sin(x*26)/13-sin(x*28)/14-sin(x*30)/15-sin(x*32)/16-sin(x*34)/17-sin(x*36)/18},
"15"= f0<-function(x) {pi/2-sin(x*2)-sin(x*4)/2-sin(x*6)/3-sin(x*8)/4-sin(x*10)/5-sin(x*12)/6-sin(x*14)/7-sin(x*16)/8-sin(x*18)/9-sin(x*20)/10-sin(x*22)/11-sin(x*24)/12-sin(x*26)/13-sin(x*28)/14-sin(x*30)/15-sin(x*32)/16-sin(x*34)/17-sin(x*36)/18-sin(x*38)/19-sin(x*40)/20}
)

switch(x,
"0"=text<-c("pi/2"),
"1"=text<-c("pi/2-sin(2x)"),
"2"=text<-c("pi/2-sin(2x)-sin(4x)/2"),
"3"=text<-c("pi/2-sin(2x)-sin(4x)/2-sin(6x)/3"),
"4"=text<-c("pi/2-sin(2x)-sin(4x)/2-sin(6x)/3-sin(8x)/4"),
"5"=text<-c("pi/2-sin(2x)-sin(4x)/2-sin(6x)/3-sin(8x)/4-sin(10x)/5"),
"6"=text<-c("pi/2-sin(2x)-sin(4x)/2-sin(6x)/3-sin(8x)/4-sin(10x)/5-sin(12x)/6"),
"7"=text<-c("pi/2-sin(2x)-sin(4x)/2-sin(6x)/3-sin(8x)/4-sin(10x)/5-sin(12x)/6-sin(14x)/7"),
"8"=text<-c("pi/2-sin(2x)-sin(4x)/2-sin(6x)/3-sin(8x)/4-sin(10x)/5-sin(12)/6-sin(14x)/7-sin(16x)/8"),
"9"= text<-c("pi/2-sin(x*2)-sin(x*4)/2-sin(x*6)/3-sin(x*8)/4-sin(x*10)/5-sin(x*12)/6-sin(x*14)/7-sin(x*16)/8-sin(x*18)/9"),
"10"=text<-c("pi/2-sin(x*2)-sin(x*4)/2-sin(x*6)/3-sin(x*8)/4-sin(x*10)/5-sin(x*12)/6-sin(x*14)/7-sin(x*16)/8-sin(x*18)/9-sin(x*20)/10"),
"11"= text<-c("pi/2-sin(x*2)-sin(x*4)/2-sin(x*6)/3-sin(x*8)/4-sin(x*10)/5-sin(x*12)/6-sin(x*14)/7-sin(x*16)/8-sin(x*18)/9-sin(x*20)/10-sin(x*22)/11-sin(x*24)/12"),
"12"= text<-c("pi/2-sin(2x)-sin(4x)/2-sin(6x)/3-sin(8x)/4-sin(10x)/5-sin(12x)/6-sin(14x)/7-sin(16x)/8-sin(18x)/9-sin(20x)/10-sin(22x)/11-sin(24x)/12-sin(26x)/13-sin(28x)/14"),
"13"= text<-c("pi/2-sin(2x)-sin(4x)/2-sin(6x)/3-sin(8x)/4-sin(10x)/5-sin(12x)/6-sin(14x)/7-sin(16x)/8-sin(18x)/9-sin(20x)/10-sin(22x)/11-sin(24x)/12-sin(26x)/13-sin(28x)/14-sin(30x)/15-sin(32x)/16"),
"14"= text<-c("pi/2-sin(2x)-sin(4x)/2-sin(6x)/3-sin(8x)/4-sin(10x)/5-sin(12x)/6-sin(14x)/7-sin(16x)/8-sin(18x)/9-sin(20x)/10-sin(22x)/11-sin(24x)/12-sin(26x)/13-sin(28x)/14-sin(30x)/15-sin(32x)/16-sin(34x)/17-sin(36x)/18"),
"15"= text<-c("pi/2-sin(2x)-sin(4x)/2-sin(6x)/3-sin(8x)/4-sin(10x)/5-sin(12x)/6-sin(14x)/7-sin(16x)/8-sin(18x)/9-sin(20x)/10-sin(22x)/11-sin(24x)/12-sin(26x)/13-sin(28x)/14-sin(30x)/15-sin(32x)/16-sin(34x)/17-sin(36x)/18-sin(38x)/19-sin(40x)/20")
)

Sawtooth_wave_x<-seq(-pi,pi,length=60)
Sawtooth_wave_y<-rep(seq(0,pi,length=30),times=2)

plot(Sawtooth_wave_x,Sawtooth_wave_y,xlim=Graph_scale_x,ylim=Graph_scale_y,type="l",col=rgb(0,0,0),main="FT(Fourier transform) Sawtooth wave", xlab="Frequency", ylab="Amplitude")
par(new=T)#上書き指定
plot(f0,xlim=Graph_scale_x,ylim=Graph_scale_y,type="l",col=rgb(1,0,0),main="", xlab="", ylab="")

legend("bottomleft", legend=c("f(t)=t-mπ,mπ<=t<(m+1)π",text),lty=c(1,1),col=c(rgb(0,0,0),rgb(1,0,0)))
}


library("animation")
Time_Code=c("0","0","0","0","0","1","1","1","2","2","2","3","3","3","4","4","4","5","5","6","7","8","9","10","11","12","13","14","15")
saveGIF({
for (i in Time_Code){
 Sawtooth_wave_exp(i)
}
}, interval = 0.1, movie.name = "TEST01.gif")

f:id:ochimusha01:20190815072440g:plain

ギブズ現象Gibbs phenomenon)」は、区分的連続微分可能な周期関数のフーリエ級数において、その関数が第1種不連続 (discontinuity of the first kind 又は jump discontinuity) となる点付近では、フーリエ級数のn 次部分和が大きく振動して、部分和の最大値が関数自体の最大値より大きくなってしまうことがあるという振る舞いのことを指す。

この超過量は、高調波の周波数(つまり、部分和の項数)が増えても無くならず、ある有限極限値に近付く。日本語表記として「ギブズの現象」、「ギブス現象」、「ギブスの現象」とされることもある。

  • 米国の物理学者アルバート・マイケルソン(Albert Abraham Michelson, 1852年〜1931年)により、グラフ作成機において最初に発見された。マイケルソンは、1898年に、フーリエ級数を計算・再合成する機械的装置を開発したが、矩形波を装置に入力すると、グラフは、不連続点付近で行ったり来たりしようとするのだった。これは、発生すると、フーリエ係数の個数が無限大に近付いても持続するようだった。

    f:id:ochimusha01:20190815203914j:plain

  • この現象を始めて数学的に説明したのが、ジョシュア・ウィラード・ギブズ(Josiah Willard Gibbs, 1839年1903年)だった。大まかな表現をするなら、この現象は、不連続関数を連続関数である正弦波関数および余弦波関数からなる級数で近似することに内在する困難の現れである。それは、また、ある関数のフーリエ係数が次数の増大に応じて減衰していく仕方が、その関数の滑らかさに従うという原則に、緊密に関係している。非常に滑らかな関数では、そのフーリエ係数は非常に急速に減衰する(そして、フーリエ級数は非常に急速に収束する)。これに対し、不連続関数では、フーリエ係数の減衰は非常に緩やかである(従って、フーリエ級数の収束は非常に緩慢である)。例えば、不連続である上記の矩形波フーリエ係数 1, 1/3, 1/5, ... は、絶対収束級数ではない調和級数程度の速さでしか減衰しない。実際、上記のフーリエ級数は、変数x のほとんど全ての値で、条件収束するだけであることが分っている。このことは、ギブス現象が何故起こるのかということの一端を説明する。それは、絶対収束するフーリエ係数を有するフーリエ級数は、ワイエルシュトラスの判定法により一様収束するから、上述のような振動を起こすことはありえないからである。同じ理由で、不連続関数は、絶対収束するフーリエ係数を持つのは不可能である。何故なら、もしそうした関数が存在したとしたら、それは、連続関数列の一様極限になるので、連続関数でなければならなくなり、矛盾が生じるからである。

    アメリコネチカット州ニューヘイブン出身の数学者・物理学者・物理化学者で、エール大学(イェール大学)教授。

    • 熱力学分野で熱力学ポテンシャル、化学ポテンシャル概念を導入し、相平衡理論の確立、相律の発見など、今日の化学熱力学の基礎を築いた。統計力学の確立にも大きく貢献した。ギブズ自由エネルギーやギブズ-デュエムの式、ギブズ-ヘルムホルツの式等にその名を残している。 ベクトル解析の創始者の一人として数学にも寄与している。

    • 科学者としての経歴は、4つの時期に分けられる。①1879年まで熱力学理論を研究。②1880年から1884年まではベクトル解析分野を研究。③1882年から1889年までは光学と光理論を研究。④1889年以降は、統計力学の教科書作成に関わる。

    • 19世紀中葉、米国の大学は、科学に殆ど関心を示さず、古典に偏重していたから、ギブズの講義は、学生の興味を殆ど引かなかった。彼の業績に興味を持ったのは、他の科学者、特に、スコットランドの物理学者ジェームズ・クラーク・マクスウェルだった。ギブズが論文を発表したのが、ヨーロッパでは余り読まれていない無名雑誌であったため、評価されるようになるのは遅く、ギブズの考えがヨーロッパで広く受け入れられたのは、論文がヴィルヘルム・オストヴァルトにより書籍の形でドイツ語訳され(1888年)、アンリ・ルシャトリエによりフランス語訳されて(1899年)からだった。

    N・ウィーナーは、ギブズが物理学の根本に統計を導入したことを極めて高く評価して「アインシュタインハイゼンベルクプランクよりギップズのほうにこそ、二十世紀物理学の最初の大革命の功績は帰せられるべきだと思う」と評した。なお、彼の功績を称えて、小惑星(2937)ギブズが彼の名を取り命名されている。

実際上は、ギブズ現象による問題は、フェイエール総和法またはリース総和法 (Riesz summation) 等のフーリエ級数の総和法における平滑化を行ったり、シグマ近似を行ったりするなら、改善できる。また、フーリエ変換の代わりに、ウェーブレット変換を用いるなら、ギブズ現象は発生しなくなる。

三角波Triangular wave)も、非正弦波的な基本的な波形の一種で、波形の見た目が三角形になっているそれ(波形)のことである。無限フーリエ級数で収束。

矩形波と同様に奇数倍音のみを含むが、矩形波に比べて高い倍音成分は急速に小さくなる(倍音の次数の逆数の自乗に比例する)。そのため、矩形波よりも聴きやすく正弦波により近い。回路で生成した三角波をローパスフィルタに通すことで正弦波を得るシンセサイザーもあった。

基本周波数に奇数倍音を合成していくことで三角波の近似を得ることができる。4n−1番目の倍音に−1をかけるか (2m±1)π または (2n±1)π によって位相をずらし、基本周波数からの相対周波数の自乗の逆数で振幅を小さくすればよい。

統計言語Rによる波形合成例

Triangular_wave_exp<-function(x){

#グラフのスケール決定
Graph_scale_x<-c(-pi,pi)
Graph_scale_y<-c(0,pi/2)
#関数と凡例の決定
switch(x,
"0"= f0<-function(x) {pi/4*x/x},
"1"= f0<-function(x) {pi/4-2*pi/pi^2*cos(2*pi*x/pi)},
"2"= f0<-function(x) {pi/4-2*pi/pi^2*(cos(2*pi*x/pi)+cos(6*pi*x/pi)/3^2)},
"3"= f0<-function(x) {pi/4-2*pi/pi^2*(cos(2*pi*x/pi)+cos(6*pi*x/pi)/3^2+cos(10*pi*x/pi)/5^2)},
"4"= f0<-function(x) {pi/4-2*pi/pi^2*(cos(2*pi*x/pi)+cos(6*pi*x/pi)/3^2+cos(10*pi*x/pi)/5^2+cos(14*pi*x/pi)/7^2)},
"5"= f0<-function(x) {pi/4-2*pi/pi^2*(cos(2*pi*x/pi)+cos(6*pi*x/pi)/3^2+cos(10*pi*x/pi)/5^2+cos(14*pi*x/pi)/7^2+cos(18*pi*x/pi)/9^2)}
)

switch(x,
"0"=text<-c("π/4"),
"1"=text<-c("π/4-2π/π^2*cos(2πx/π)"),
"2"=text<-c("π/4-2π/π^2*(cos(2πx/π)+cos(6πx/π)/3^2)"),
"3"=text<-c("π/4-2π/π^2*(cos(2πx/π)+cos(6πx/π)/3^2+cos(10πx/π)/5^2)"),
"4"=text<-c("π/4-2π/π^2*(cos(2πx/π)+cos(6πx/π)/3^2+cos(10πx/π)/5^2+cos(14πx/π)/7^2)"),
"5"=text<-c("π/4-2π/π^2*(cos(2πx/π)+cos(6πx/π)/3^2+cos(10πx/π)/5^2+cos(14πx/π)/7^2+cos(18πx/π)/9^2)"),
)

Triangular_wave_x<-seq(-pi,pi,length=60)
Triangular_wave_y<-c(seq(0,pi/2,length=15),seq(pi/2,0,length=15),seq(0,pi/2,length=15),seq(pi/2,0,length=15))

plot(Triangular_wave_x,Triangular_wave_y,xlim=Graph_scale_x,ylim=Graph_scale_y,type="l",col=rgb(0,0,0),main="FT(Fourier transform) Triangular wave", xlab="Frequency", ylab="Amplitude")
par(new=T)#上書き指定
plot(f0,xlim=Graph_scale_x,ylim=Graph_scale_y,type="l",col=rgb(1,0,0),main="", xlab="", ylab="")

legend("bottomleft", legend=c("f(t)=t(0=<t<T/2), T-t(T/2<=t<T))",text),lty=c(1,1),col=c(rgb(0,0,0),rgb(1,0,0)))
}


library("animation")
Time_Code=c("0","0","0","0","0","1","1","1","2","2","2","3","3","3","4","4","4","5","5")
saveGIF({
for (i in Time_Code){
 Triangular_wave_exp(i)
}
}, interval = 0.1, movie.name = "TEST01.gif")

f:id:ochimusha01:20190815110144g:plain

矩形波くけいは、Square wave)も、非正弦波的な基本的な波形の一種で、電子工学や信号処理の分野で広く使われている。方形波とも呼ばれる。

理想上は2レベルの間を規則的かつ瞬間的に変化するが、その2レベルにはゼロが含まれることも含まれないこともある。

デジタルスイッチング回路で広く使われており、2進法(2レベル)の論理回路から生成される。厳密に定められた間隔で同期論理回路を動作させるためには矩形波の高速な遷移が適しているので、タイミングの基準や「クロック信号」に使用されている。しかしながら周波数領域グラフからわかるように、矩形波は広帯域の周波数成分を含んでいる。これらは電磁放射や電流のパルスを発生させてしまい、近くの回路に影響を及ぼしその結果、雑音や誤りを引き起こす。精密なAD変換器のような非常に敏感な回路などではこの問題を避けるために、矩形波の代わりに正弦波をタイミング基準として使用する。

フーリエ級数を使用して、無限級数の形で理想的な波形を表すことができるが、この表現で興味深いのは矩形波でも起こる「ギブズ現象Gibbs phenomenon)」である。のこぎり波でも起こる。

  • 理想的ではない矩形波で生じるリンギングは、この現象と関連している。シグマ近似(sigma-approximation)を用いることで防ぐことができ、系列がより滑らかに収束するように「Lanczos sigma factors」を使うこともできる。
  • 理想的な矩形波には、信号の高値から低値への変化が切れ味よく瞬間的であることが求められるが、そのためには無限の帯域幅が必要となるので、現実の世界での達成は不可能である。

現実の世界には有限の帯域幅しか存在しないため、ギブズ現象に似たリンギングやシグマ近似に似たリップル (電気)が起きることがよくある。

統計言語Rによる波形合成例

Square_wave_exp<-function(x){

#グラフのスケール決定
Graph_scale_x<-c(-pi,pi)
Graph_scale_y<-c(-1/2,1)
#関数と凡例の決定
switch(x,
"0"= f0<-function(x) {1/2*x/x},
"1"= f0<-function(x) {1/2+2/pi*sin(x)},
"2"= f0<-function(x) {1/2+2/pi*sin(x)+2/(3*pi)*sin(3*x)},
"3"= f0<-function(x) {1/2+2/pi*sin(x)+2/(3*pi)*sin(3*x)+2/(5*pi)*sin(5*x)},
"4"= f0<-function(x) {1/2+2/pi*sin(x)+2/(3*pi)*sin(3*x)+2/(5*pi)*sin(5*x)+2/(7*pi)*sin(7*x)},
"5"= f0<-function(x) {1/2+2/pi*sin(x)+2/(3*pi)*sin(3*x)+2/(5*pi)*sin(5*x)+2/(7*pi)*sin(7*x)+2/(9*pi)*sin(9*x)},
"6"= f0<-function(x) {1/2+2/pi*sin(x)+2/(3*pi)*sin(3*x)+2/(5*pi)*sin(5*x)+2/(7*pi)*sin(7*x)+2/(9*pi)*sin(9*x)+2/(11*pi)*sin(11*x)},
"7"= f0<-function(x) {1/2+2/pi*sin(x)+2/(3*pi)*sin(3*x)+2/(5*pi)*sin(5*x)+2/(7*pi)*sin(7*x)+2/(9*pi)*sin(9*x)+2/(11*pi)*sin(11*x)+2/(13*pi)*sin(13*x)},
"8"= f0<-function(x) {1/2+2/pi*sin(x)+2/(3*pi)*sin(3*x)+2/(5*pi)*sin(5*x)+2/(7*pi)*sin(7*x)+2/(9*pi)*sin(9*x)+2/(11*pi)*sin(11*x)+2/(13*pi)*sin(13*x)+2/(15*pi)*sin(15*x)},
"9"= f0<-function(x) {1/2+2/pi*sin(x)+2/(3*pi)*sin(3*x)+2/(5*pi)*sin(5*x)+2/(7*pi)*sin(7*x)+2/(9*pi)*sin(9*x)+2/(11*pi)*sin(11*x)+2/(13*pi)*sin(13*x)+2/(15*pi)*sin(15*x)+2/(17*pi)*sin(17*x)}
)

switch(x,
"0"=text<-c("1/2"),
"1"=text<-c("1/2+2/π*sin(x)"),
"2"=text<-c("1/2+2/π*sin(x)+2/3π*sin(3x)"),
"3"=text<-c("1/2+2/π*sin(x)+2/3π*sin(3x)+2/5π*sin(5x)"),
"4"=text<-c("1/2+2/π*sin(x)+2/3π*sin(3x)+2/5π*sin(5x)+2/7π*sin(7x)"),
"5"=text<-c("1/2+2/π*sin(x)+2/3π*sin(3x)+2/5π*sin(5x)+2/9π*sin(9x)"),
"6"=text<-c("1/2+2/π*sin(x)+2/3π*sin(3x)+2/5π*sin(5x)+2/9π*sin(9x)+2/11π*sin(11x)"),
"7"=text<-c("1/2+2/π*sin(x)+2/3π*sin(3x)+2/5π*sin(5x)+2/9π*sin(9x)+2/11π*sin(11x)+2/13π*sin(13x)"),
"8"=text<-c("1/2+2/π*sin(x)+2/3π*sin(3x)+2/5π*sin(5x)+2/9π*sin(9x)+2/11π*sin(11x)+2/13π*sin(13x)+2/15π*sin(15x)"),
"9"=text<-c("1/2+2/π*sin(x)+2/3π*sin(3x)+2/5π*sin(5x)+2/9π*sin(9x)+2/11π*sin(11x)+2/13π*sin(13x)+2/15π*sin(15x)+2/17π*sin(17x)")
)

Square_wave_x<-seq(-pi,pi,length=60)
Square_wave_y<-c(rep(0,times=30),rep(1,times=30))
plot(Square_wave_x,Square_wave_y,xlim=Graph_scale_x,ylim=Graph_scale_y,type="l",col=rgb(0,0,0),main="FT(Fourier transform) Square wave", xlab="Frequency", ylab="Amplitude")
par(new=T)#上書き指定
plot(f0,xlim=Graph_scale_x,ylim=Graph_scale_y,type="l",col=rgb(1,0,0),main="", xlab="", ylab="")

legend("bottomleft", legend=c("f(t)=0(-π=<t<0), 1(0<=t<π))",text),lty=c(1,1),col=c(rgb(0,0,0),rgb(1,0,0)))
}


library("animation")
Time_Code=c("0","0","0","0","0","1","1","1","2","2","2","3","3","3","4","4","4","5","5","6","7","8","9")
saveGIF({
for (i in Time_Code){
 Square_wave_exp(i)
}
}, interval = 0.1, movie.name = "TEST01.gif")

f:id:ochimusha01:20190815123000g:plain

放物線(希:παραβολή/parabolē、羅/英:parabola、独:Parabel)、すなわちf(t)=t^2(-a<=t<=a)合を合成する場合。

統計言語Rによる波形合成例

Parabola_wave_exp<-function(x){

#グラフのスケール決定
Graph_scale_x<-c(-pi,pi)
Graph_scale_y<-c(-pi/2,pi^2)
#関数と凡例の決定
switch(x,
"0"= f0<-function(x) {pi^2/3*x/x},
"1"= f0<-function(x) {pi^2/3+4*(-cos(x))},
"2"= f0<-function(x) {pi^2/3+4*(-cos(x)+cos(2*x)/2^2)},
"3"= f0<-function(x) {pi^2/3+4*(-cos(x)+cos(2*x)/2^2-cos(3*x)/3^2)},
"4"= f0<-function(x) {pi^2/3+4*(-cos(x)+cos(2*x)/2^2-cos(3*x)/3^2+cos(4*x)/4^2)},
"5"= f0<-function(x) {pi^2/3+4*(-cos(x)+cos(2*x)/2^2-cos(3*x)/3^2+cos(4*x)/4^2-cos(5*x)/5^2)},
"6"= f0<-function(x) {pi^2/3+4*(-cos(x)+cos(2*x)/2^2-cos(3*x)/3^2+cos(4*x)/4^2-cos(5*x)/5^2+cos(6*x)/6^2)},
"7"= f0<-function(x) {pi^2/3+4*(-cos(x)+cos(2*x)/2^2-cos(3*x)/3^2+cos(4*x)/4^2-cos(5*x)/5^2+cos(6*x)/6^2-cos(7*x)/7^2)},
"8"= f0<-function(x) {pi^2/3+4*(-cos(x)+cos(2*x)/2^2-cos(3*x)/3^2+cos(4*x)/4^2-cos(5*x)/5^2+cos(6*x)/6^2-cos(7*x)/7^2+cos(8*x)/8^2)},
"9"= f0<-function(x) {pi^2/3+4*(-cos(x)+cos(2*x)/2^2-cos(3*x)/3^2+cos(4*x)/4^2-cos(5*x)/5^2+cos(6*x)/6^2-cos(7*x)/7^2+cos(8*x)/8^2)-cos(9*x)/9^2}
)

switch(x,
"0"=text<-c("a^2/3"),
"1"=text<-c("a^2/3+4a^2/π^2(-cos(πx/a))"),
"2"=text<-c("a^2/3+4a^2/π^2(-cos(πx/a)+cos(2πx/a)/2^2"),
"3"=text<-c("a^2/3+4a^2/π^2(-cos(πx/a)+cos(2πx/a)/2^2-cos(3πx/a)/3^2"),
"4"=text<-c("a^2/3+4a^2/π^2(-cos(πx/a)+cos(2πx/a)/2^2-cos(3πx/a)/3^2+cos(4πx/a)/4^2"),
"5"=text<-c("a^2/3+4a^2/π^2(-cos(πx/a)+cos(2πx/a)/2^2-cos(3πx/a)/3^2+cos(4πx/a)/4^2-cos(5πx/a)/5^2"),
"6"=text<-c("a^2/3+4a^2/π^2(-cos(πx/a)+cos(2πx/a)/2^2-cos(3πx/a)/3^2+cos(4πx/a)/4^2-cos(5πx/a)/5^2+cos(6πx/a)/6^2"),
"7"=text<-c("a^2/3+4a^2/π^2(-cos(πx/a)+cos(2πx/a)/2^2-cos(3πx/a)/3^2+cos(4πx/a)/4^2-cos(5πx/a)/5^2+cos(6πx/a)/6^2-cos(7πx/a)/7^2"),
"8"=text<-c("a^2/3+4a^2/π^2(-cos(πx/a)+cos(2πx/a)/2^2-cos(3πx/a)/3^2+cos(4πx/a)/4^2-cos(5πx/a)/5^2+cos(6πx/a)/6^2-cos(7πx/a)/7^2+cos(8πx/a)/8^2"),
"9"=text<-c("a^2/3+4a^2/π^2(-cos(πx/a)+cos(2πx/a)/2^2-cos(3πx/a)/3^2+cos(4πx/a)/4^2-cos(5πx/a)/5^2+cos(6πx/a)/6^2-cos(7πx/a)/7^2+cos(8πx/a)/8^2-cos(9πx/a)/9^2")
)

Parabola_x<-seq(-pi,pi,length=60)
Parabola_y<-Parabola_x^2
plot(Parabola_x,Parabola_y,xlim=Graph_scale_x,ylim=Graph_scale_y,type="l",col=rgb(0,0,0),main="FT(Fourier transform) Parabola function", xlab="Frequency", ylab="Amplitude")
par(new=T)#上書き指定
plot(f0,xlim=Graph_scale_x,ylim=Graph_scale_y,type="l",col=rgb(1,0,0),main="", xlab="", ylab="")

legend("bottomleft", legend=c("f(t)=t^2(-a<=t<=a)",text),lty=c(1,1),col=c(rgb(0,0,0),rgb(1,0,0)))
}


library("animation")
Time_Code=c("0","0","0","0","0","1","1","1","2","2","2","3","3","3","4","4","4","5","5","6","7","8","9")
saveGIF({
for (i in Time_Code){
 Parabola_wave_exp(i)
}
}, interval = 0.1, movie.name = "TEST01.gif")

f:id:ochimusha01:20190815171145g:plain