R言語でサウンド入門まとめておきます。
pythonはこれまでも何度かまとめていますが、比較のために記述します。
合わせて、非周期系時系列データである、日経225の爆上げの音を聞いてみようと思います。
参考は以下のとおりです。参考③にtuneRのインストールの仕方があります。
①How to get started with sound in R (pdf)
②Play sound in Python
③Basic Sound Processing with R
###・音にして出力する
環境
Windows10
RStudio;Version 1.2.5042
python3.6.10
mic;hdmi
###pythonで音を出す
これまでも何度かやってきましたが、前回用意したcsvファイルや数値列を以下のコードで音として、wavファイルで保存し、出力します。
####csv書込みと読込
※csv書込み前段の処理は前回の記事を見てください
csv書込みは以下のとおり、headerをつけてしかも音が持続するように、20回同じ数値列を書込みすることとしました。
import csv
with open('pl3_seasonal.csv', 'w') as f:
fieldnames = ['seasonal']
writer = csv.DictWriter(f, fieldnames=fieldnames)
writer.writeheader()
writer = csv.writer(f, delimiter='\n')
for i in range(20):
writer.writerow(pl3)
pandasでもう少し汎用的に4つのデータを格納して利用することが以下で出来ました。
df = pd.DataFrame({'seasonal': pl3,
'trend': pl2,
'resid':pl4,
'observed':pl1})
df.to_csv('seasonal_trend.csv', index=False)
以下で20回分を追記しています。
for i in range(20):
df.to_csv('seasonal_trend.csv', mode='a', header=False, index=False)
csvファイル読込は以下のとおりです。
df = pd.read_csv('seasonal_trend.csv')
pl3_seasonal = df['seasonal']
####音ファイル作成と再生
本編の音ファイル作成と保存は以下のコードで実施します。
import wave
import numpy as np
import struct
import pyaudio
fs = 800 #サンプリング周波数
#サイン波を-32768から32767の整数値に変換(signed 16bit pcmへ)
pl3_n = 32767*pl3_seasonal/max(pl3_seasonal)
swav = [int(x ) for x in pl3_n]
#バイナリ化
binwave = struct.pack("h" * len(swav), *swav)
#サイン波をwavファイルとして書き出し
w = wave.Wave_write("pl3_seasonal.wav")
params = (1, 2, fs, len(binwave), 'NONE', 'not compressed')
w.setparams(params)
w.writeframes(binwave)
w.close()
そして音の再生は以下のコードで実施できます。以下は上で保存した'pl3_seasonal.wav'を再生します。
※いろいろ変更もできますが過去記事見てください
from playsound import playsound
playsound('pl3_seasonal.wav')
####Rで音を出す
Rで音出力は初物ですが、簡単にできるんですね。
単純にディレクトリにある音声ファイル’mywave.wav’を読み込んで音を発生させます。
library(tuneR) #load tuenR package
mwav = readWave("mywave.wav")
play(mwav)
以下はサイン波で波を計算して、wavファイルに変換し、それを出力しています。
#make a simple sine wave and play
t = seq(0, 3, 1/8000)
u = (2^15-1)*sin(2*pi*440*t)
w = Wave(u, samp.rate = 8000, bit=16)
play(w)
以下で一度'w_1.wav'として保存し、さらに読み込んで再生しています。
writeWave(w, 'w_1.wav')
rw = readWave("w_1.wav")
play(rw)
上記で数値列があればそこからwavファイルを作成して保存できることが分かりました。
そこで以下では、時系列の要素分解したdecompose_seasonalデータをwavファイルに変換して音として出力してみます。そして、最後にwavファイルを保存します。
ds <- read.csv(file = 'decompose_seasonal.csv')
ds_w = Wave(32678/max(abs(ds))*ds, samp.rate = 8000, bit=16)
play(ds_w)
writeWave(ds_w, 'ds_w.wav')
※データを非周期的なデータに変えて以下に全体を説明します
###日経平均の爆上げの音を聞く
非周期的なデータだけど、一応振動しつつ変化するものの代表格に株価があります。ということで、このところの株価変動をdecomposeして音に変換して、この音を聞いてみましょう。
※これは楽しみでやっているのでマネしないでください
一応、Rのコードで実施します。
nikkei_225の最近のデータをダウロードして、準備します。
hz=18は、周期はないですがいろいろな周期間隔で分析できるように変数として定義しておきます。データ数がhzの整数倍とは限らないので、lsを導入してデータ数を調整しています。
※hzを変更していろいろなグラフや音を発生させます
decompose()は今までと同様です。
data <- read.csv("./industrial/nikkei_225_.csv",encoding = "utf-8")
hz=18
ls = round(908/hz)*hz
IIP <- ts(data,start=c(1),frequency=1)
IIP=IIP[0:ls]
IIP <- ts(IIP,start=c(1),frequency=hz)
decompose_IIP <- decompose(IIP)
trend, random, seasonal, observed(IIP)は以下のとおり
※seasonalは補助線入れときます
このseasonalを音にしようということです。
plot(decompose_IIP$trend)
plot(decompose_IIP$random)
plot(decompose_IIP$seasonal[0:72])
lines(decompose_IIP$seasonal[0:72])
| seasonal | random |
|---|---|
 |
 |
| trend | observed |
 |
 |
| 上と同じ絵だが、以下のようにm_beerを平均することにより、seasonal_beerを求め、これを3個出力してみる。 |
m_beer = t(matrix(data = decompose_IIP$seasonal, nrow = hz))
seasonal_beer = colMeans(m_beer, na.rm = T)
plot(as.ts(rep(seasonal_beer,3)))
上のグラフと異なり、線グラフで得られた。中身は同じだ。
これを、130個分としてcsvに出力した。
write.csv(as.ts(rep(seasonal_beer,130)), file = 'decompose_seasonal130.csv')
ここからは、上記の音生成と同じコードである。すなわち以下のようにcsvファイルを読み込む。
ds <- read.csv(file = 'decompose_seasonal130.csv')
ここでも、以下音発生にはtuneRライブラリを利用する。
こうして、音が発生できた。
音は最大にしたいので、pythonと同じように数値を規格化している。
library(tuneR)
ds_w = Wave(32678/max(abs(ds[2]))*ds[2], samp.rate = 150*hz, bit=16)
play(ds_w)
最後に、ファイルに保存する。
file_n <- paste('ds_',hz,'.wav')
writeWave(ds_w, file_n)
こうして、hzを変化すると、その代表的な音の波形は以下のとおり変化する。つまり、あるところで母音を再現したように見える。
実際の音は波形以上に類似して聞こえるから文字通り耳を疑う。
| hz | 近い母音 | 波形 | 母音 |
|---|---|---|---|
| 18 | あ |  |
 |
| 40 | え |  |
 |
| 128 | い |  |
 |
| 16 | お |  |
 |
| 17 | う |  |
 |
| ####得られた音データを以下に置きました | |||
| ・R_Audio/data/ | |||
| 上記配下の | |||
| R_Audio/data/1_a_ds_ 18 .wav | |||
| などをダウンロードして再生してみてください。 | |||
| ####データについて | |||
| なんでこんなことになっているかは不思議な気がしますので、ちょっとデータ構造を見てみます。 | |||
| 以下のようにhz=16として、16個ずつのグループにしてdecompose()します。 | |||
| まず、データ出力すると |
hz=16
ls = round(908/hz)*hz
IIP <- ts(data,start=c(1),frequency=1)
IIP=IIP[0:ls]
IIP <- ts(IIP,start=c(1),frequency=hz)
print(IIP[0:ls])
> source('C:/R_data/decompose_audio.R')
[1] 16820 17200 16930 16970 16820 16870 17030 16560 16800 16790 17410 17510 16870 16790 17230 17280
[17] 17460 17780 18180 17980 17700 17120 17220 17560 17520 17710 17920 17920 17200 17140 17240 17270
．．．
[897] 20070 19940 19920 19930 20010 19940 19990 20010 20050 20000 20020 20010 NA NA NA NA
そして、肝心な部分は以下のとおり
decompose_IIP <- decompose(IIP)
print(decompose_IIP$seasonal)
つまり、decompose_IIP$seasonalは以下のように16個ずつ平均取られています。
Time Series:
Start = c(1, 1)
End = c(57, 16)
Frequency = 16
[1] 16.659497 17.697443 16.413961 16.955256 -14.905540 -7.621449 13.458097 10.605824
[9] -2.302557 -11.895191 4.550122 17.479809 -19.166396 -26.945414 -21.738941 -9.244521
．．．
[897] 16.659497 17.697443 16.413961 16.955256 -14.905540 -7.621449 13.458097 10.605824
[905] -2.302557 -11.895191 4.550122 17.479809 -19.166396 -26.945414 -21.738941 -9.244521
グラフにプロットすると、以下のようになっています。
プロットを数えてみれば、16個で周期になっているのが分かります。
このデータを再生しているわけです。
そして、大切なことは株価データを周期的とみて適当な個数でグループとして平均化すると、あるデータでは母音と同じようなデータ列を与えてくれるということになりました。
完全な乱数列に対して同じ処理をしても、このような結果になりませんが、たまたま株価データはこういう構造を内包していると考えることが出来ると思います。
※本当のところは、もっと深堀しないとわかりませんが…
ああ~、この方法だとオリジナルの爆上げ音は聞けないですね．．．
※もっとも、そんなもの元々ないと思いますけどね．．．振幅が大きくなるのは音に表せると思いますが、．．．次回に予知・予兆として是非やってみたいと思います
###まとめ
・RとPythonのサウンド入門を比較してみた
・非周期的なデータである株価をdecomposeし、周期的なデータとなった季節変動を音に変換して生成してみた
・適当な周期でサンプリングすることにより、母音と似たような音発生をすることが分かった
・もう少し工夫して、decompose()と音発生を利用した予知・予兆をやってみようと思う