この投稿以降「音の話」が気になって…

#Rで離散フーリエ変換を試す。

準備としてとりあえず、以下のお題をそのまま解いて見ます。

Rによる信号処理

振幅が1で周波数10[Hz]の正弦波と，振幅が1/4で周波数30[Hz]の余弦波の合成波信号データについて，離散フーリエ変換（FFT）によって得られた振幅スペクトルを取得しなさい．ただし，サンプリング周波数FSは160[Hz]，標本間隔DTは6.25[ms]（1/FS），標本の個数Nは16個とする．

統計言語Rによるプログラミング例（初期値設定）

FS <- 160    # サンプリング周波数 FS=160[Hz]
DT <- 1/FS  # 標本間隔 DT=1/160=0.00625[s]
N <- 16        # データ数 N=16
i <- 0:(N-1)   # 添字(index)の範囲 i

DT16=0.0062516=0.1[s]を基本波の周期とする。
従ってその逆数である10[Hz]が基本波の周波数となる。

スペクトラム表示関数

Spectrum_disp<-function(x,N){
i <- 0:(N-1) 
spec <- fft(x)/N
plot(i, Re(spec), type='h', xaxt='n', yaxt='n', col='blue', lty=3, lwd=2, ylim=c(-0.8,0.8), ylab='', xlab='')
par(new=T) #上書き
plot(i, Im(spec), type='h', xaxt='n', yaxt='n', col='red', lty=1, lwd=2, ylim=c(-0.8,0.8), ylab=expression(paste('Spectrum ', X[k])), xlab='Index k')
axis(1, pos=0, at=0:N, adj=0)
axis(2, at=seq(-1, 1, by=0.25))
legend(0.5, 0.7, legend=c(expression(Re(X[k])), expression(Im(X[k]))), lty=c(3,1), col=c('blue', 'red'), lwd=c(2,2))
}

振幅1，10[Hz]の正弦波

siga <- sin(2*pi*10*DT*i)
plot(DT*i, siga, xlab='Time', type='o')  
#sigaのFFT結果
Spectrum_disp(siga,N)

k=1の場合が基本周波数10[Hz]
奇関数Sinの成分なので虚部にスペクトルが現れる。
ちょうどk=8を折り返し周波数として半分ずつの数値がk=1とk=15の位置に対照的に現れる。

振幅0.25，30[Hz]の余弦波(sigb)

sigb <- 0.25*cos(2*pi*30*DT*i)
plot(DT*i, sigb, xlab='Time', type='o') 
#sigaのFFT結果
Spectrum_disp(sigb,N)

k=3の場合が周波数30[Hz]
偶関数Cosの成分なので実部にスペクトルが現れる。
ちょうどk=8を折り返し周波数として半分ずつの数値がk=3とk=13の位置に対照的に現れる。

合成波(sig)

sig <- siga + sigb
plot(DT*i, sig, xlab='Time', type='o') 
#sigのFFT結果
Spectrum_disp(sig,N)

スペクトルの内容はまさにsiga+sigbになっている。

スペクトルの3D表示

library(rgl)
Samples_z<-i
Real_x<-Re(fft(sig))
Imagin_y<-Im(fft(sig))
plot3d(Real_x,Imagin_y,Samples_z)  
movie3d(spin3d(axis=c(0,0,1),rpm=5),duration=10,fps=25,movie="~/Desktop/test")
#Movieファイル出力ディレクトリーは環境依存

ここで一応、用語解説。
サンプリング周波数
 帯域幅 - Wikipedia

#平均律と純正律
まずは基本中の基本から抑えて行きます。
和音と波形 | You Look Too Cool

とりあえず音の高さを1、周期を1とします。長３度和音ドミソは純正律（Just intonation）だときれいな協和音になりますが、この場合の周波数は次のとおりです。

ド 1.0000
ミ 5/4=1.2500
*ソ 3/2=1.5000

一方、平均律のドミソは次のとおりです。

ド 1.0000
ミ 2^(4/12)=1.2599
ソ 2^(7/12)=1.4983

若干、ずれが生じます。したがって音にも影響します。ちなみにこれは明るい感じのするメジャー和音ですが、ドレソ#は暗い感じがします。

純正律（Just Intonation）

ド 1.0000
レ 9/8=1.2500
ソ# 10/9=1.111111

平均律（Equal Temperament）

ド 1.0000
レ 2^(2/12) =1.122462
ソ# 2^(8/12) =1.587401

純正律 - Wikipedia

統計言語Rによる純正律再生

library(tuneR)
#純正律再生関数
Just_Intonation<-function(note){
f1<-function(x){(440*2^(-3/12))*x} #チューンダウン
notes<-f1(note)
for (i in notes) {
beep<- sine(i);
play(beep)
}
}
#ドレミファソラシド、ドシラソファミレド
c0<-c(1,9/8,5/4,4/3,3/2,5/3,15/8,2)
c1<-c(c0,rev(c0))
Just_Intonation(c1)
#ドミソ、ドミソ
c0<-c(1,5/4,3/2)
c1<-c(c0,c0)
Just_Intonation(c1)
#ドレソ#、ドレソ#
c0<-c(1,9/8,10/9+16/15-1)
c1<-c(c0,c0)
Just_Intonation(c1)
#ドミ♭ソ、ドミ♭ソ
c0<-c(1,5/4+16/15-1,3/2)
c1<-c(c0,c0)
Just_Intonation(c1)

統計言語Rによる平均律再生

library(tuneR)
#平均律再生関数
Equal_Temperament<-function(note){
f1<-function(x){(440*2^(-3/12))*2^x} #チューンダウン
notes<-f1(note)
for (i in notes) {
beep<- sine(i);
play(beep)
}
}

#ドレミファソラシド、ドシラソファミレド
c0<-c(0,2/12,4/12,5/12,7/12,9/12,11/12,1)
c1<-c(c0,rev(c0))
Equal_Temperament(c1)
#ドミソ、ドミソ
c0<-c(1,4/12,7/12)
c1<-c(c0,c0)
Equal_Temperament(c1)
#ドレソ#、ドレソ#
c0<-c(1,2/12,8/12)
c1<-c(c0,c0)
Equal_Temperament(c1)
#ドミ♭ソ、ドミ♭ソ
c0<-c(1,3/12,7/12)
c1<-c(c0,c0)
Equal_Temperament(c1)

東京都高等学校軽音楽連盟

平均律は周波数を12回、同じ倍率で大きくすると２倍になるような各１回を半音としています。具体的な倍率は２の12乗根で、1.059463倍です。(計算機でこの数字を12回掛けてみてください)つまり、各半音の幅がきっちり同じなわけです。

純正律は、ある調専用(例えばハ長調限定)の調律法で、長３度和音(ドミソド)の周波数が
完全に整数倍になるように合わせます。平均律で合わせた場合はこの周波数が合いません。

ギターで弾くとそれほどの違いは感じられませんし、単音で聞いてこの違いがわかる人もほとんどいないと思われますが、正弦波を発振する装置で両方の和音を聞き比べるとはっきり違いがわかります。チューニングの項で書いた「うなり」が、純正律だと全く発生しないわけです。ただし、ハ長調の純正律に合わせた楽器で別の調の曲を弾くと、平均律以上に狂いが目立ってきますので、いろいろな曲をひとつの楽器で弾く場合は平均律に合わせるのが常識となっています。

なおギターのチューニングは、各フレットが平均律で作られていて、ハーモニックスチューニングをすると弦間が純正律という中途半端な形になります。

それでは離散フーリエ変換によって解析してみましょう。

統計言語Rによるプログラミング例（FFT解析結果表示関数）

Spectrum_disp<-function(x){
N<-length(x)
i <- 0:(N-1) 
spec <- fft(x)/N
plot(i, Re(spec), type='h', xaxt='n', yaxt='n', col='blue', lty=3, lwd=2, ylim=c(-0.8,0.8), main='Spectrum by FFT',ylab='Spectrum', xlab='Samples')
par(new=T) #上書き
plot(i, Im(spec), type='h', xaxt='n', yaxt='n', col='red', lty=1, lwd=2, ylim=c(-0.8,0.8), main='',ylab='', xlab='')
axis(1, pos=0, at=0:N, adj=0)
axis(2, at=seq(-1, 1, by=0.25))
legend(0.5, 0.7, legend=c("Real Number","Imaginary Number"), lty=c(3,1), col=c('blue', 'red'), lwd=c(2,2))
}

統計言語Rによるプログラミング例（初期値設定）

N<-12# データ数 N=12
i<-0:(N-1)# 添字(index)の範囲 i
j<-seq(0,2*pi,length=12)# 周期関数のX

基本波形「ド(Do)」の解析

Do<-sin(j)
plot(i, Do, xlab='Time', xlim=c(0,11),ylim=c(-1,1),main="Sampling of Do",type='o')

純正律(Just intonation)「ドミソ」の合成

Do_j<-sin(j)
Mi_j<-sin(j*5/4)
plot(i, Mi_j, xlab='Samples', xlim=c(0,11),ylim=c(-1,1),main="Sampling of Mi(Just intonation)",type='o')
So_j<-sin(j*3/2)
plot(i, So_j, xlab='Samples', xlim=c(0,11),ylim=c(-1,1),main="Sampling of So(Just intonation)",type='o')
DoMiSo_j<-Do_j+Mi_j+So_j
plot(i, DoMiSo_j, xlab='Samples', xlim=c(0,11),ylim=c(-3,3),main="Sampling of DoMiSo(Just intonation)",type='o')
Spectrum_disp(DoMiSo_j)

#スペクトラムの3D表示
library(rgl)
Samples_z<-i
Real_x<-Re(fft(DoMiSo_j))
Imagin_y<-Im(fft(DoMiSo_j))
plot3d(Real_x,Imagin_y,Samples_z)  
movie3d(spin3d(axis=c(0,0,1),rpm=5),duration=10,fps=25,movie="~/Desktop/test")
#Movieファイル出力ディレクトリーは環境依存

平均律(Equal Temperament)「ドミソ」の合成

Do_e<-sin(j)
Mi_e<-sin(j*2^(4/12))
plot(i, Mi_e, xlab='Samples', xlim=c(0,11),ylim=c(-1,1),main="Sampling of Mi(Equal Temperament)",type='o')
Spectrum_disp(Mi_e)
So_e<-sin(j*2^(7/12))
plot(i, So_e, xlab='Samples', xlim=c(0,11),ylim=c(-1,1),main="Sampling of So(Equal Temperament)",type='o')
Spectrum_disp(So_e)
DoMiSo_e<-Do_e+Mi_e+So_e
plot(i, DoMiSo_e, xlab='Samples', xlim=c(0,11),ylim=c(-3,3),main="Sampling of DoMiSo(Equal Temperament)",type='o')
Spectrum_disp(DoMiSo_e)

#3D表示
library(rgl)
Samples_z<-i
Real_x<-Re(fft(DoMiSo_e))
Imagin_y<-Im(fft(DoMiSo_e))
plot3d(Real_x,Imagin_y,Samples_z)  
movie3d(spin3d(axis=c(0,0,1),rpm=5),duration=10,fps=25,movie="~/Desktop/test")
#Movieファイル出力ディレクトリーは環境依存

とりあえず自分史上、このジャンル最初の投稿を果たしたという事で、以下続報…
離散フーリエ変換(DFT)の仕組みを完全に理解する
 離散時間フーリエ変換の困るところ
これからしばらくはこういう記事に目を通し続ける日々が続きます。