サンプル

ここには、音声関連の話で使用したR言語のスクリプトを載せます。

レイアウト

「グラフィック窓」は、リサイズすると、中の描画もきれいに追従します。余白を保っているようです。線幅は変更されません。

■余白

最初は、かなり大きな余白があります。コピーして使うのに余白を詰めておきます。

# グラフの表示のマージンを小さく設定しておく
par(mar=c(4,4,1,1),mfcol=c(1,1))

■「グラフィック窓」の分割

「グラフィック窓」を、縦横任意の数に分割して表示できるらしいが、1つの図を描くように初期化しておきます。

#　レイアウトを初期化しておく
layout(matrix(c(1,1), 1, 1, byrow = TRUE))

複素数のデータ列の2段表示用の関数

実部、虚部をそれぞれ2段で表示する。表示後、段組を解除するため関数としました。

#　複素数を2段で表示する関数
plot2 <- function(t,fft){
  layout(matrix(c(1,1,2,2), 2, 2, byrow = TRUE))
  mn<-min(Im(fft))
  mx<-max(Im(fft))
  plot(t,Re(fft),ylim=c(mn,mx),col=3)
  plot(t,Im(fft),ylim=c(mn,mx),col=2)
  layout(matrix(c(1,1), 1, 1, byrow = TRUE))
}

等差数列を作り、sin()の値を求める

FFTのテストデータを生成します。seq()で等差数列を作り、その各値のsin()をテストデータに使いました。
plot()のtypeは、線、マークの形状をしています。colは色を、xlabはX軸の名前を与えます。
cal="red"のようにも書けます。

# 下記のような16サンプルを得たとする。
data<-sin(seq(0,by=2*pi/16,length.out=16))
plot(0:15,data,type="o",col=2,xlab="")

同じデータをX軸の表示を変える

前述のデータは、1Hzのサイン波を1秒間、16Hzでサンプリングしたも」のとも、「100Hzのサイン波を0.01秒間、1600Hzでサンプリングしたもの」とも見なせますが、前者でX軸のメモリを入れてみます。

# このデータは、1Hzのサイン波を1秒間、16Hzでサンプリングしたものと言える
sec<-seq(0,by=1/16,length.out=16)
plot(sec,data,type="o")

前述のテストデータをfft()した結果を2段で表示

前半に書いた関数plot2()を呼んで、fft()の結果を、実部、虚部に分けてプロットします。

# このデータを、fftすると以下のようになる
fft<-fft(data)/length(data)
plot2(0:15,fft)

複数の折れ線を引く

同じ周期のcos()を合成して、sin()、cos()、合成結果を書きます。

# おなじ周期のcosを合成する
cos<-cos(seq(0,by=2*pi/16,length.out=16))
data2<-(data+cos)/2
plot(0:15,data,type="l",col=2)
lines(0:15,cos,type="l",col=3)
lines(0:15,data2,type="o",col=1)

exp()を描画する

(実部、虚部)の位置にマークを書きます。

# 円を8等分するような角度を生成
a<-seq(0,by=2*pi/16,length.out=16)
# その値を虚部とする複素数を生成
c_a<-complex(imaginary=a)
# そのexp()を求める
exp_ti<-exp(c_a)
#　その(実部、虚部)にマークを書く
plot(Re(exp_ti),Im(exp_ti),col="red")
#　(0,0)を通る線を引く
abline(v=0,h=0,col="gray")

スクリプト全体

# グラフの表示のマージンを小さく設定しておく
par(mar=c(4,4,1,1),mfcol=c(1,1))
# レイアウトを初期化しておく
layout(matrix(c(1,1), 1, 1, byrow = TRUE))
#　複素数を2段で表示する関数
plot2 <- function(t,fft){
  layout(matrix(c(1,1,2,2), 2, 2, byrow = TRUE))
  mn<-min(Im(fft))
  mx<-max(Im(fft))
  plot(t,Re(fft),ylim=c(mn,mx),col=3)
  plot(t,Im(fft),ylim=c(mn,mx),col=2)
  layout(matrix(c(1,1), 1, 1, byrow = TRUE))
}

# 下記のような16サンプルを得たとする。
data<-sin(seq(0,by=2*pi/16,length.out=16))
plot(0:15,data,type="o",col="red",xlab="")

# このデータは、1Hzのサイン波を1秒間、16Hzでサンプリングしたものと言える
sec<-seq(0,by=1/16,length.out=16)
plot(sec,data,type="o")

# また、100Hzのサイン波を0.01秒間、1600Hzでサンプリングしたものとも言える
sec<-seq(0,by=1/1600,length.out=16)
plot(sec,data,type="o")

# このデータを、fftすると以下のようになる
fft<-fft(data)/length(data)
plot2(0:15,fft)

# おなじ周期のcosを合成する
cos<-cos(seq(0,by=2*pi/16,length.out=16))
data2<-(data+cos)/2
plot(0:15,data,type="l",col=2)
lines(0:15,cos,type="l",col=3)
lines(0:15,data2,type="o",col=1)

# これをfft
fft2<-fft(data2)/length(data2)
plot2(0:15,fft2)

# 円を8等分するような角度を生成
a<-seq(0,by=2*pi/16,length.out=16)
# その値を虚部とする複素数を生成
c_a<-complex(imaginary=a)
# そのexp()を求める
exp_ti<-exp(c_a)
#　その(実部、虚部)にマークを書く
plot(Re(exp_ti),Im(exp_ti),col="red")
#　(0,0)を通る線を引く
abline(v=0,h=0,col="gray")

##########################################
# 100,300,800Hzのサイン波の合成
# 100Hzのサイン波を1000Hzで1秒間サンプリングしたデータ
d1<-sin(seq(0,by=2*pi*100/1000,length.out=1000))
# 300Hzのサイン波を1000Hzで1秒間サンプリングしたデータ
d2<-sin(seq(0,by=2*pi*300/1000,length.out=1000))
# 800Hzのサイン波を1000Hzで1秒間サンプリングしたデータ
d3<-sin(seq(0,by=2*pi*800/1000,length.out=1000))
# 合成
d<-(d1+d2+d3)/3
# プロット
plot(1:10,d[1:10],type="o",ylim=c(-1,1))
lines(1:10,d1[1:10],col="red")
lines(1:10,d2[1:10],col="green")
lines(1:10,d3[1:10],col="blue")
##########################################
# 補助線用に800Hzのサイン波を10KHzで10msサンプリングした状態
d4<-sin(seq(0,by=2*pi*800/10000,length.out=100))
x<-seq(1,by=0.1,length.out=100);
plot(1:10,d3[1:10],type="o",col="blue")
lines(x,d4,type="l",col="green")
##########################################
# 合成波形のfft()
fft<-fft(d)/length(d)
mat <- matrix(c(1,1,2,2), 2, 2, byrow = TRUE)
layout(mat)
plot(0:999,Re(fft),col="green4")
plot(0:999,Im(fft),col="red")
##########################################
# パワースペクトルのプロット
Hz<-c(0:498)
PowerSpectrum=matrix(500)
for(i in 0:499){
y<-abs(fft[i])
PowerSpectrum[i]<-y*y}
plot(Hz,PowerSpectrum,type="l",col="red")

軸の表示

# サンプルデータ
sin<-sin(seq(0,by=2*pi/8,length.out=8))
x<-c("0","(1/4)π","(2/4)π","(3/4)π","(4/4)π","(5/4)π","(6/4)π","(7/4)π")
plot(0:7,sin,type="o",axes=FALSE,xlab="",col="firebrick")
axis(1,at=0:7,labels=x)
axis(2,at=-1:1)
abline(h=0,col="gray")

音の話で使ったサンプル

# グラフの表示のマージンを小さく設定しておく
par(mar=c(4,4,1,1),mfcol=c(1,1))

# 平均律を9オクターブ分作る
N<- -57:50 #基準から-57音、50音。計108音
equal<- 440*(2^(N/12)) #440Hzの前後を計算

# 純正律を9オクターブ分作る
r<- c(1,16/15,9/8,6/5,5/4,4/3,17/12,3/2,8/5,5/3,9/5,15/8)
c<- 264/16 # 440Hzを含むところから4オクターブ下のド
just<-matrix(nrow=1)
for(i in 1:9) {
  for(j in 1:12) {
      just[(i-1)*12+j] <- c * r[(j-1)%%12+1]
  }
  c <- c * 2 # 次のド
}

# 差を百分率でプロット
percentage <- (just-equal)/equal*100
plot(N,percentage,type="l")

ラウドネスレベル曲線のようなもの

# 計算する周波数の並び
Hz=c(20,25,31.5,40,50,63,80,100,125,160,200,250,315,400,500,630,800,1000,1250,1600,2000,2500,3150,4000,5000,6300,8000,10000,12500)
# 周波数依存係数
ll=c(-31.6,-27.2,-23.0,-19.1,-15.9,-13.0,-10.3,-8.1,-6.2,-4.5,-3.1,-2.0,-1.1,-0.4,0.0,0.3,0.5,0.0,-2.7,-4.1,-1.0,1.7,2.5,1.2,-2.1,-7.1,-11.2,-10.7,-3.1)
fa=c(0.532,0.506,0.480,0.455,0.432,0.409,0.387,0.367,0.349,0.330,0.315,0.301,0.288,0.276,0.267,0.259,0.253,0.250,0.246,0.244,0.243,0.243,0.243,0.242,0.242,0.245,0.254,0.271,0.301)
f1=c(78.5,68.7,59.5,51.1,44.0,37.5,31.5,26.5,22.1,17.9,14.4,11.4,8.6,6.2,4.4,3.0,2.2,2.4,3.5,1.7,-1.3,-4.2,-6.0,-5.4,-1.5,6.0,12.6,13.9,12.3)
loudness_curve <- function(phon){
  D <- 0.00447 * (10^(0.025*phon) - 1.15)
  d1 <- (0.4*10^(((f1+ll)/10)-9 ))
  d1 <- d1^fa
  D <- D + d1
  f <- (10/fa) * log10(D) - ll + 94
}
dB<-loudness_curve(80)
plot(Hz,dB,type="l",log="x",ylim=c(0,120))
lines(Hz,loudness_curve(60))
lines(Hz,loudness_curve(40))
lines(Hz,loudness_curve(20))

十二平均律の各音程の周波数の差

N <- c(-57:50) # 108音分の音程番号
f <- 440 * 2^(N/12)
diff <- f[2:108]-f[1:107]
plot(N[1:107],diff,type="l")

音程を判定するのに必要な時間

N <- c(-57:50) # 108音分の音程番号
f <- 440 * 2^(N/12)
diff <- f[2:108]-f[1:107]
d2<-diff/2
sec<-1/d2
plot(N[1:108],sec,type="l")

re <- matrix(0,16)
im <- matrix(0,16)
s <- sin( seq(0,by=2*pi/16,length.out=32) )
for(i in 1:16){
  k <- i+15
  f <- fft( s[i:k] )/16
  re[i] <- Re(f[2])
  im[i] <- Im(f[2])
}
plot(re,type="l",col="green",ylab="")
par(new=T)
plot(im,type="l",col="red",ylab="FFT Result[2] Real / Imaginary")
legend(13,y=0.4,c("実部","虚部"), col = c("green","red"), lty=1, merge = TRUE, bg='gray90')
title(main="位相の表現")

re <- matrix(0,16)
im <- matrix(0,16)
s <- sin( seq(0,by=2*pi/16,length.out=32) )
for(i in 1:16){
  k <- i+15
  f <- fft( s[i:k] )/16
  re[i] <- Re(f[2])
  im[i] <- Im(f[2])
}
plot(re,im)
segments(0,-0.52,0,0.52,col="gray")
segments(-0.52,0,0.52,0,col="gray")
title(main="位相を変えたFFTの結果の実部、虚部")

#################################################################
#データを作る。
d <- sin(seq(0,by=3*pi/7,length.out=16)) #サイン波の16サンプル
s <- matrix(0,48)                        #s[1:48]にデータを作る
for(k in 1:16) s[k] <- d[k]              #最初の16サンプル
for(k in 1:16) s[16+k] <- 0              #中間の16サンプルをゼロに
for(k in 1:16) s[32+k] <- d[k]           #最後の16サンプル
windows(width=4,height=3)
plot(s,type="l")                         #確認ためプロット
title(main="サンプリングデータのつもり")
#################################################################
#FFT（サンプル数で割る）
f0 <- fft(s)/length(s)
windows(width=4,height=3)
plot(Re(f0),col="darkgreen")
title(main="FFTの結果-実部")
windows(width=4,height=3)
plot(Im(f0),col="red")
title(main="FFTの結果-虚部")
#################################################################
#FFTの逆変換で、復元を確認
f1 <- fft(f0,inv=T)
windows(width=4,height=3)
plot(Re(f1),type="l")
title(main="FFTの逆変換で、復元を確認")
#################################################################
#FFTの結果から、合成（全部）
i <- complex(im=1)
y <- matrix(0,48)
t <- 0
for(k in 1:48){
  for(j in 1:48){
    y[k] <- y[k] + Re(f0[j])*exp(i*(j-1)*t) + Im(f0[j])*exp(i*(j-1)*t)*i
  }
  t <- t + 2*pi/48
}
windows(width=4,height=3)
plot(Re(y),type="l")
title(main="FFTの結果を全部使って合成")
#################################################################
#FFTの結果から、合成(部分)
i <- complex(im=1)
y <- matrix(0,48)
t <- 0
for(k in 1:48){
  y[k] <-        Re(f0[10])*exp(i* 9*t) + Im(f0[10])*exp(i* 9*t)*i
  y[k] <- y[k] + Re(f0[11])*exp(i*10*t) + Im(f0[11])*exp(i*10*t)*i
  y[k] <- y[k] + Re(f0[12])*exp(i*11*t) + Im(f0[12])*exp(i*11*t)*i
  y[k] <- y[k] + Re(f0[13])*exp(i*12*t) + Im(f0[13])*exp(i*12*t)*i
  y[k] <- y[k] + Re(f0[37])*exp(i*36*t) + Im(f0[37])*exp(i*36*t)*i
  y[k] <- y[k] + Re(f0[38])*exp(i*37*t) + Im(f0[38])*exp(i*37*t)*i
  y[k] <- y[k] + Re(f0[39])*exp(i*38*t) + Im(f0[39])*exp(i*38*t)*i
  y[k] <- y[k] + Re(f0[40])*exp(i*39*t) + Im(f0[40])*exp(i*39*t)*i
  t <- t + 2*pi/48
}
windows(width=4,height=3)
plot(Re(y),type="l",col=1)
title(main="FFTの結果の一部を合成(8,10箇所)")
#もう2箇所追加
t <- 0
for(k in 1:48){
  y[k] <- y[k] + Re(f0[ 9])*exp(i* 8*t) + Im(f0[ 9])*exp(i* 8*t)*i
  y[k] <- y[k] + Re(f0[41])*exp(i*40*t) + Im(f0[41])*exp(i*40*t)*i
  t <- t + 2*pi/48
}
par(new=T)
plot(Re(y),type="l",col=2,axes=F,ylab="",xlab="")
#################################################################
#FFTの結果から、合成(半分で波形が決まる)
i <- complex(im=1)
y <- matrix(0,48)
t <- 0
for(k in 1:48){
  y[k] <-        Re(f0[10])*exp(i* 9*t) + Im(f0[10])*exp(i* 9*t)*i
  y[k] <- y[k] + Re(f0[11])*exp(i*10*t) + Im(f0[11])*exp(i*10*t)*i
  y[k] <- y[k] + Re(f0[12])*exp(i*11*t) + Im(f0[12])*exp(i*11*t)*i
  y[k] <- y[k] + Re(f0[13])*exp(i*12*t) + Im(f0[13])*exp(i*12*t)*i
  t <- t + 2*pi/48
}
windows(width=4,height=3)
plot(Re(y),type="l")
title(main="FFTの結果の半分で波形が決まる")
#################################################################
#FFTの結果から、合成(三角関数に)
i <- complex(im=1)
y <- matrix(0,48)
t <- 0
for(k in 1:48){
  y[k]= (2*Re(f0[10])*cos(9*t) + 2*Im(f0[10])*sin(9*t)
      + 2*Re(f0[11])*cos(10*t) + 2*Im(f0[11])*sin(10*t)
      + 2*Re(f0[12])*cos(11*t) + 2*Im(f0[12])*sin(11*t)
      + 2*Re(f0[13])*cos(12*t) + 2*Im(f0[13])*sin(12*t))
  t <- t + 2*pi/48
}
windows(width=4,height=3)
plot(y,type="l")
title(main="三角関数で合成")

x<-seq(from=-1,to=1,by=0.05)
abs1<-abs(x)
log1<-log(abs1)
plot(x,x,type="l",axes=T,ylim=c(-3,1),xlab="スペクトルの値",ylab="変換前後の値")
par(new=T)
plot(x,log1,type="l",axes=F,ylim=c(-3,1),xlab="",ylab="",col=2)
abline(h=-3:1,v=seq(-1,1,by=0.5),col="gray")
legend(0.4,y=-1.5,c("FFT結果の値","算出値"), col = c(1,2), lty=1, merge = TRUE, bg='gray100')
title(main="log(abs(x))の効果")

par(mar = c(4, 4, 5, 4))　#　右の軸のためにマージンを広く
x1 <- seq(from=0,to=32767,by=1000)
x2 <- seq(from=0,to=1,by=1/(length(v1)-1))
plot(x1, log(x1),type="l",col=1,
axes=F,xlab="非正規化振幅",ylab="非正規化振幅の対数")
axis(1) # X軸(上)
axis(2) # Y軸(左)
par(new=T) # 重ねがきをする
plot(x2,log(x2),type="p",col=2,axes=F,xlab="",ylab="")
axis(3,col=2,col.axis=2) # X軸(上) 色を赤に。col.axisは文字の色
axis(4,col=2,col.axis=2) # Y軸(右) 色を赤に。colは、線の色
mtext("正規化振幅",3,line=2,col=2) # X軸(上)のタイトル
mtext("正規化振幅の対数",4,line=2,col=2) # Y軸(右)のタイトル
title(main="入力の正規化とその対数",line=3)

# 0から7の連番をデータとする
c2 <- complex(re=seq(0,by=1,length.out=8),im=0)
# 指数関数部分を計算しておく
ep <- exp(complex(im=-seq(0,by=2*pi/8,length.out=8)))
f1 <- c2*ep # 1周期目は、指数関数を順番に掛ける
ep2 <- c(ep[1],ep[8],ep[7],ep[6],ep[5],ep[4],ep[3],ep[2])
f7 <- c2*ep2
windows(3.5,4) # FFTの結果の状態の表示 #############################
plot(c2,pch=c("0","1","2","3","4","5","6","7"),xlim=c(-8,8), ylim=c(-8,8),xlab="実部",ylab="虚部")
par(new=T)
plot(f1,pch=c("0","1","2","3","4","5","6","7"),font=2,xlim=c(-8,8),ylim=c(-8,8),axes=F,col=2,xlab="",ylab="")
par(new=T)
plot(f7,pch=c("0","1","2","3","4","5","6","7"),font=2,xlim=c(-8,8),ylim=c(-8,8),col=3,axes=F,xlab="",ylab="")
abline(h=0,v=0,col="gray")
title(main="連番に[1][7]の指数関数乗算")
windows(3.5,4) # 共役を考慮して、再度、指数関数を乗算すると #########
plot(f1,pch=c("0","1","2","3","4","5","6","7"),xlim=c(-8,8),ylim=c(-8,8),col=2,xlab="実部",ylab="虚部")
par(new=T)
plot(f1*ep2,pch=c("0","1","2","3","4","5","6","7"),font=2,xlim=c(-8,8),ylim=c(-8,8),col=1,axes=F,xlab="",ylab="")
abline(h=0,v=0,col="gray")
title(main="再度、指数関数を乗算[1]")
windows(3.5,4) # 共役を考慮して、再度、指数関数を乗算すると #########
plot(f7,pch=c("0","1","2","3","4","5","6","7"),xlim=c(-8,8),ylim=c(-8,8),col=3,xlab="実部",ylab="虚部")
par(new=T)
plot(f7*ep,pch=c("0","1","2","3","4","5","6","7"),font=2,xlim=c(-8,8),ylim=c(-8,8),col=1,axes=F,xlab="",ylab="")
abline(h=0,v=0,col="gray")
title(main="再度、指数関数を乗算[7]")

リストの要素、文字列の分割

C#でカナ文字入力ダイアログを書いたのですが、縦長すぎるので「アイオエオ」を「アカサタナ」にしたいと思いました。
R言語で並べ替えて使いました。

> d <- c(" アイウエオ",
+        " カキクケコ",
+        " サシスセソ",
+        " タチツテト",
+        " ナニヌネノ",
+        " マミムメモ",
+        " ヤユヨ＿＿",
+        " ラリルレロ",
+        " ワン＿＿＿",
+        " ガギグゲゴ",
+        " ザジズゼゾ",
+        " ダヂヅデド",
+        " パピプペポ")
> e <- strsplit(d," ")
> for(i in 2:6){
+   s <- ""
+   for(j in 1:13){
+     s <- paste(s,e[[j]][i])
+   }
+   print(s)
+ }
[1] " アカサタナマヤラワガザダパ"
[1] " イキシチニミユリンギジヂピ"
[1] " ウクスツヌムヨル＿グズヅプ"
[1] " エケセテネメ＿レ＿ゲゼデペ"
[1] " オコソトノモ＿ロ＿ゴゾドポ"
>

３Dの分布図

bm_read() は、作成したエクステンション。
scatterplot3dは、パッケージで「パッケージ」「パッケージのインストール」「パッケージの読み込み」をして使う。

#3D ヒストグラム
source("Bitmap_R.R")
img <- bm_read("12.jpg")

s <- array(0,dim=c(256,256,256))
for(i in 1:dim(img$rgb)[2]){
  x<-img$rgb[1,i]+1
  y<-img$rgb[2,i]+1
  z<-img$rgb[3,i]+1
  s[x,y,z]<-s[x,y,z]+1
}
limit <- 200
i <- 0
cx <- NULL
cy <- NULL
cz <- NULL
count<-NULL
for(x in 1:256) for(y in 1:256) for(z in 1:256){
  if(s[x,y,z]>limit) {
    cx[i] <- x
    cy[i] <- y
    cz[i] <- z
    count[i] <- s[x,y,z];
    i <- i + 1
  }
}
pt <- rainbow(max(count))
cols <- pt[count]
scatterplot3d(cx, cy, cz, color = cols,
      xlab="R:赤", ylab="G：緑", zlab="B：青",
      box = FALSE, angle = 120,
      xlim = c(min(cx), max(cx)), ylim = c(min(cy), max(cy)), zlim = c(min(cz), max(cz)),
      main=paste("RGBの度数 780x767 ",as.integer(limit),"超 ",
                 as.integer(sum(count)),"/",as.integer(img$width*img$height)))

Lab*色度図

# CIE L*a*b* 色度図
w <- 50
h <- 50
for(l in seq(25,95,length.out=6))
{
  bm <- array(0,c(3, w * h))
  i <- 0
  for(a in seq(-100,100,length.out=h))
    for(b in seq(-100,100,length.out=w)) {
      c <- hex(LAB(l,a,b))
      if(!is.na(c))
        bm[,i] <- col2rgb(c)
      i = i + 1
    }
  windows(4,4)
  bm_image(bm, w, h)
}

Lab* 3D分布図、分離画像

#3D LAB化
source("Bitmap_R.R")
img <- bm_read("12.jpg")
cols <- rgb(img$rgb[1,],img$rgb[2,],img$rgb[3,],maxColor=255)
rgbc <- hex2RGB(cols)
labc <- as(rgbc,"LAB")
#度数
s <- array(0,dim=c(100,200,200))
for(i in 1:length(labc@coords[,1])){
  l <- round(labc@coords[i,1])
  a <- round(labc@coords[i,2]+100)
  b <- round(labc@coords[i,3]+100)
  s[l,a,b]<-s[l,a,b]+1
}
limit <- 60
i <- 0
cL <- NULL
ca <- NULL
cb <- NULL
count<-NULL
for(l in 1:100) for(a in 1:200) for(b in 1:200){
  if(s[l,a,b]>limit) {
    cL[i] <- l
    ca[i] <- a - 100
    cb[i] <- b - 100
    count[i] <- s[l,a,b];
    i <- i + 1
  }
}
pt <- rainbow(max(count))
cols <- pt[count]
scatterplot3d(ca, cb, cL, color = cols,
      xlab="a*", ylab="b*", zlab="L*",
      box = FALSE, angle = 120,
      xlim = c(min(ca), max(ca)), ylim = c(min(cb), max(cb)), zlim = c(min(cL), max(cL)),
      main=paste("L*a*b*の度数 780x767 ",as.integer(limit),"超 ",
                 as.integer(sum(count)),"/",as.integer(img$width*img$height)))

windows(4,4)
plot(count,ylab="度数",main="L*a*b*での度数")
windows(4,4)
plot(count,ylim=c(0,3000),ylab="度数",main="L*a*b*での度数(3000以下）")

# a*がゼロ未満の部分を描画。その他は白
windows(5,5)
b2 <- matrix(0,nrow=3,ncol=length(img$rgb[1,]))
b2[1,] <- ifelse(labc@coords[,2]<0,img$rgb[1,],255)
b2[2,] <- ifelse(labc@coords[,2]<0,img$rgb[2,],255)
b2[3,] <- ifelse(labc@coords[,2]<0,img$rgb[3,],255)
bm_image(b2,img$width,img$height)
# a*がゼロ以上の部分を描画。その他は白
windows(5,5)
b3 <- matrix(0,nrow=3,ncol=length(img$rgb[1,]))
b3[1,] <- ifelse(labc@coords[,2]>=0,img$rgb[1,],255)
b3[2,] <- ifelse(labc@coords[,2]>=0,img$rgb[2,],255)
b3[3,] <- ifelse(labc@coords[,2]>=0,img$rgb[3,],255)
bm_image(b3,img$width,img$height)
# L*をグレースケールと見なして描画
windows(5,5)
b4 <- matrix(0,nrow=3,ncol=length(img$rgb[1,]))
b4[1,] <- labc@coords[,1]*2.55
b4[2,] <- labc@coords[,1]*2.55
b4[3,] <- labc@coords[,1]*2.55
bm_image(b4,img$width,img$height)

RGLパッケージで3D作図

図の上でマウスを動かすと視点が変えられる。
左ボタンを押してマウスを動かすと移動、右ボタンだとズーミングする。

z<-c(0,0,0)
x<-c(5,0,0)
y<-c(0,5,0)
rz<-array(0,41)
r<-seq(0,2*pi,by=((2*pi)/40))
rx<-10*cos(r)
ry<-10*sin(r)
open3d()
triangles3d(x,y,z,col='magenta')
lines3d(c(-10,10),c(0,0),c(0,0),col='blue')
lines3d(c(0,0),c(-10,10),c(0,0),col='red')
lines3d(c(0,0),c(0,0),c(-10,10),col='green')
lines3d(c(0,0),c(-10,0),c(0,10),col='darkgray')
lines3d(c(0,0),c(10,0),c(0,10),col='darkgray')
lines3d(rx,ry,rz,col='gray')
text3d(0,0,10,'eye',col='Red')
axes3d()
view3d(theta = 30)

RGB色相

#########################################
# 色相のサンプル
#########################################
#上り坂
up_y <- seq(0,by=(255/60),length.out=60)
#下り坂
down_y <- seq(255,by=-(255/60),length.out=60)
#水平
py0 <- array(0,120)
py1 <- array(255,120)
#R
r <- c(py1[1:60],down_y,py0,up_y,py1[1:60])
#G
g <- c(up_y,py1,down_y,py0)
#B
b <- c(py0,up_y,py1,down_y)
#LINES
windows(5,4)
plot(c(0:359), r, type="l", lwd=2, col="red", xlab="色相(度)",ylab="RGBの各値")
lines(c(0:359), g, type="l", lwd=2, col="green")
lines(c(0:359), b, type="l", lwd=2, col="blue")
#Colors
for(x in 0:359)
{
lines(array(x,32),c(112:143), col=rgb(r[x+1],g[x+1],b[x+1],max=255), type="l", lwd=2)
}
#
title(main="RGBの各値と色相")

彩度

color <- rainbow(11)
windows(4,4)
plot(0:255,0:255,type="n",xlab="max",ylab="min",main="レンジと彩度（分母：最大値）")
for(x in 1:255)
{
for(y in 0:x)
{
v <- (x-y)/x
points(x,y,col=color[v*10])
}
}
windows(4,4)
plot(0:255,0:255,type="n",xlab="max",ylab="min",main="レンジと彩度（分母：255-|255-(M+m)|）"
)
for(x in 1:255)
{
for(y in 0:x)
{
v <- (x-y)/(255 - abs(255 - (x+y)))
points(x,y,col=color[v*10])
}
}
n<-11
windows(4,3)
image(0:n,0:0,matrix(1:n,ncol=1),col=rainbow(n),
main="彩度の色表示",
axes=F,xlab="",ylab="")
axis(1, c(0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5), col.axis = "blue")

Amazon.co.jp ウィジェット

mikeo_410