この記事は、58th TokyoRでのLT内容の焼き直しです。

intro

音の解析って面白そうですよね。

Aさん(ソプラノ)の「あ」とBさん（テノール）の「あ」は、どちらも「あ」として脳に認識されます。それは「い」とは異なる。当たり前に思いますが、では「あ」と「い」の境目に線を引いてくださいと言われると、悩みが大きそうです。例えば、「æ」(aとeが組み合わさった文字)は、日本語の「あ」と「え」の中間的な発音をするわけですよ。catのaの発音ですね。

これは、音の特徴量をどうやって定義するか、という問題に汎化できます。

という事で、Webに転がっている音のサンプルをネタに、音解析のデモを少々。

dataset

HeartSound & Murmur Library @ Univ. Michigan
心拍音の公開データです(.mp3形式)。
ダウンロードできます。

取り扱うのは一番上の「Apex Area - Supine, Listening with the bell of stethoscope」にします。心臓の心尖部（僧帽弁領域）の音をbell型聴診器で聞いた音だそうです。5種類の患者さんの音と、対照2名の音があります。

> head(list.files("Michigan_Heart_Sounds"), 7)
[1] "01 Apex, Normal S1 S2, Supine, Bell.mp3"  # Normal
[2] "02 Apex, Split S1, Supine, Bell.mp3"      # Normal
[3] "03 Apex, S4, LLD, Bell.mp3"               # 僧帽弁逸脱症(MVP)
[4] "04 Apex, Mid Sys Click, Supine, Bell.mp3" # 急性僧帽弁逆流症
[5] "05 Apex, S3, LLD, Bell.mp3"               # 冠動脈疾患による僧帽弁逆流症(MR)
[6] "06 Apex, Early Sys Mur, Supine, Bell.mp3" # MVPによるMR
[7] "07 Apex, Mid Sys Mur, Supine, Bell.mp3"   # 古典的MRもしくは心室中隔欠損症

では、この8データを分離できるか試してみましょう。

私は心臓に関してはド素人なので、興味のある方はダウンロード先や、下記をご参照ください。

参考：心音の聴診

library

今回は、下記を使います。いずれもCRAN版ですので、install.packages("hoge")で落とせます。

音を使うパッケージ2つと、何かと捗る{tidyverse}。

library("tuneR")
library("seewave")
library("tidyverse")

参考：俺たちのtidyverseはこれからだ！

ちなみに、

> library(tidyverse)

--------------------------------------------
filter(): dplyr, stats
lag():    dplyr, stats
select(): dplyr, MASS
spec():   readr, seewave

と出てきますね。
これらの関数を使うときは、明示的にpackage名::関数名とした方が安全です。

data import

A <- list.files("Michigan_Heart_Sounds")
k <- 1

# read file
dat <- readMP3(paste("Michigan_Heart_Sounds/", A[k], sep=""))

# data structure
str(dat)
> Formal class 'Wave' [package "tuneR"] with 6 slots
  ..@ left     : int [1:2873088] 0 0 0 0 0 0 0 0 0 0 ...
  ..@ right    : num(0) 
  ..@ stereo   : logi FALSE
  ..@ samp.rate: num 44100
  ..@ bit      : num 16
  ..@ pcm      : logi TRUE

この音は、stereoじゃなくてmonauralで、そのamplitudeが@leftに入っていて、16bitデータで、時間分解能は44100Hzだという事がわかります。65秒間ぐらいのデータですね。

切り出し

データの先頭を切り出して、プロットしてみましょう。

start <- 0 * dat@samp.rate
end <- 1.6 * dat@samp.rate

dat1 <- dat[start:end]

# amplitude
par(mar = c(3, 5, 0.5, 0.5), tcl = -0.2, mgp = c(1.5, 0.3, 0))
plot(dat1)

# spectrogram
spectro(dat1, 
        collevels = seq(-50, 0, length=20),  # dB range
        palette = reverse.gray.colors.1,     # color palette
        flim = c(0, 1))                      # kHz range

「ドッ」(I音)と「ックン」(II音)が２回繰り返されていますね。

extract signal

{seewave}のtimerを使って背景ノイズとシグナルを分離します。
このデータはかなりノイズが少ないので楽ですね。
ここではcutoff値を4.5%にしています。

# cutoff: 4.5% of signal
threshold <- 4.5

par(mar = c(3, 3, 0.5, 0.5), tcl = -0.2, mgp = c(1.5, 0.3, 0))
dat_hb <- timer(dat1, threshold = threshold, msmooth = c(50, 0))

長い無音部分は取れていますが、シグナル部分が細く分離していますね。
あと「ドッ」(I音)と「ックン」(II音)が分かれているのでこれは一塊で取り出したい。

> str(dat_hb)
List of 6
 $ s      : num [1:34] 0.04879 0.00454 0.0034 0.00113 0.00113 ...
 $ p      : num [1:35] 0.31319 0.00113 0.00227 0.00454 0.24397 ...
 $ r      : num 0.155
 $ s.start: num [1:34] 0.313 0.363 0.37 0.378 0.623 ...
 $ s.end  : num [1:34] 0.362 0.368 0.373 0.379 0.624 ...
 $ first  : chr "pause"

parameter
s	duration of signal period(s) in seconds
p	duration of pause period(s) in seconds
r	ratio between the signal and silence periods(s)
positions	a list containing four elements:
s.start	start position(s) of signal period(s)
s.end	end position(s) of signal period(s)
first	whether the first event detected is a pause or a signal

いろいろ調整すると、無音継続時間が<0.275secぐらいだと、ちょうど「ドックン」が取り出せました。泥臭いやり方です↓。

# ignor first and last pause
dat_hb$p[1] <- dat_hb$p[1]+1
dat_hb$p[length(dat_hb$p)] <- dat_hb$p[length(dat_hb$p)]+1

# ignor pause dulation < 0.275 sec
threshold_pause <- 0.275
a <- c(1:length(dat_hb$p))[dat_hb$p > threshold_pause]

# start and end of signals
list_start <- dat_hb$s.start[a]
list_end <- dat_hb$s.end[a-1]

プロットすると、

x <- 1:length(dat1@left)/dat@samp.rate
y <- x

# cutin
for(ii in 1:length(list_end)){
  y[y >= list_start[ii] & y <= list_end[ii]] <- 100
}
# cutoff
y[y < 100] <- 0

# plot
par(mar = c(3, 3, 0.5, 0.5), tcl = -0.2, mgp = c(1.5, 0.3, 0))
plot(x, y, type = "l", col = "Red", yaxt = "n")
par(new = T)
plot(dat1)

このシグナル部分を別の.wavにoutputしておきましょう。

i <- 1
dat_i <- dat[list_start[i]:list_end[i]]
writeWave(dat_i, paste("hoge", i, ".wav", sep =""))

7種類のデータからそれぞれ「ドックン」を取り出すとこんな感じ。
（全て同じ値のthresholdで取り出せます。）

statistical properties

{seewave}のspecpropを使って音の特徴量を算出しましょう。
これを全部の「ドックン」について実行して、結果を集計します。

prop_dat <- specprop(seewave::spec(dat_i, f = dat_i@samp.rate, plot =F))

str(prop_dat)
> List of 14
 $ mean    : num 336
 $ sd      : num 577
 $ median  : num 321
 $ sem     : num 5.99
 $ mode    : num 61.8
 $ Q25     : num 200
 $ Q75     : num 411
 $ IQR     : num 212
 $ cent    : num 336
 $ skewness: num 9.41
 $ kurtosis: num 108
 $ sfm     : num 0.0028
 $ sh      : num 0.595
 $ prec    : Named num 2.38
  ..- attr(*, "names")= chr "x"

parameters
mean	mean frequency (see mean)
sd	standard deviation of the mean (see sd)
sem	standard error of the mean
median	median frequency (see median)
mode	mode frequency, i.e. the dominant frequency
Q25	first quartile (see quantile)
Q75	third quartile (see quantile)
IQR	interquartile range (see IQR)
cent	centroid, see note
skewness	skewness, a measure of asymmetry, see note
kurtosis	kurtosis, a measure of peakedness, see note
sfm	spectral flatness measure (see sfm)
sh	spectral entropy (see sh)
prec	frequency precision of the spectrum

multivariate analysis

集計するとこんな感じになるはずです(以後の計算に使ったパラメータだけ表示)。
1rowが1ドックンですね。

> head(round(dat_all,3))
  datID    mean       sd   mode kurtosis   sfm    sh dulation
1     1 310.228 1773.836 77.604 1023.705 0.006 0.428    0.387
2     1 313.497 1788.627 77.833 1014.725 0.007 0.425    0.386
3     1 304.034 1750.256 77.833  997.389 0.007 0.425    0.386
4     1 306.773 1761.094 78.062  988.200 0.008 0.426    0.384
5     1 313.893 1791.810 77.604 1025.223 0.007 0.427    0.387
6     1 305.961 1761.365 77.604 1012.312 0.006 0.426    0.387

ここまできたら、もうお好きにどうぞ、ですね。

prcompでPCAをかけてもいいし、isoMDSでMDSをかけてもいいですね。
方法論とRでの実装はここでは説明しません。
資料は色々ありますし、教科書ですと多次元データ解析法がお気に入りです。

どちらでも良い分離が得られていますね。
NormalのID1, ID2が接近しているのは納得ですね。

appendix

> sessionInfo()
R version 3.3.1 (2016-06-21)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.10.5 (Yosemite)

locale:
[1] ja_JP.UTF-8/ja_JP.UTF-8/ja_JP.UTF-8/C/ja_JP.UTF-8/ja_JP.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] purrr_0.2.2     readr_1.0.0     tidyr_0.6.0     tibble_1.2     
 [5] tidyverse_1.0.0 ggplot2_2.2.1   dplyr_0.5.0     tuneR_1.3.1    
 [9] seewave_2.0.5   MASS_7.3-45