「Ｒ」フリーな統計処理言語かつ環境を利用した、超音波の音圧測定データ解析

「R」を利用した、 超音波の音圧解析操作資料 超音波テスター（音圧測定・解析システム） 超音波システム研究所 1

音圧測定データ（***.csv）の解析操作 「R」をダブルクリックして立ち上げる https://youtu.be/2RcXz_xtNu4 2

解析用テキストファイル（解析用テキストｃｈＡ．ｔｘｔ）を開く https://youtu.be/uo2PCRF2xIo https://youtu.be/GR7NLxUfU4A 3

解析プログラム１（測定データ32個の平均・分散・最大・最小の変化） par(mfrow=c(3,2)) data08 <- read.table("C:/sample/sample-0001_01.csv", skip=6, sep=",", nrows=9000) data11 <- na.omit(data08)#欠損(NA)を除く m01 <- mean(data11$V2)#平均 v01 <- var(data11$V2)#分散 s01 <- min(data11$V2)#最小値 t01 <- max(data11$V2)#最大値 data08 <- read.table("C:/sample/sample-0001_02.csv", skip=6, sep=",", nrows=9000) data11 <- na.omit(data08)#欠損(NA)を除く m02 <- mean(data11$V2)#平均 v02 <- var(data11$V2)#分散 s02 <- min(data11$V2)#最小値 t02 <- max(data11$V2)#最大値 ・・・・・・ data08 <- read.table("C:/sample/sample-0001_32.csv", skip=6, sep=",", nrows=9000) data11 <- na.omit(data08)#欠損(NA)を除く m32 <- mean(data11$V2)#平均 v32 <- var(data11$V2)#分散 s32 <- min(data11$V2)#最小値 t32 <- max(data11$V2)#最大値 gdata1 <- c(m01,m02,m03,m04,m05,m06,m07,m08,m09,m10,m11,m12,m13,m14,m15,m16,m17,m1 8,m19,m20,m21,m22,m23,m24,m25,m26,m27,m28,m29,m30,m31,m32) gdata2 <- c(v01,v02,v03,v04,v05,v06,v07,v08,v09,v10,v11,v12,v13,v14,v15,v16,v17,v1 8,v19,v20,v21,v22,v23,v24,v25,v26,v27,v28,v29,v30,v31,v32) gdata3 <- c(s01,s02,s03,s04,s05,s06,s07,s08,s09,s10,s11,s12,s13,s14,s15,s16,s17,s1 8,s19,s20,s21,s22,s23,s24,s25,s26,s27,s28,s29,s30,s31,s32) gdata4 <- c(t01,t02,t03,t04,t05,t06,t07,t08,t09,t10,t11,t12,t13,t14,t15,t16,t17,t1 8,t19,t20,t21,t22,t23,t24,t25,t26,t27,t28,t29,t30,t31,t32) plot(gdata1,type="l", xlab="時間", ylab="平均") plot(gdata2,type="l", xlab="時間", ylab="分散") plot(gdata3,type="l", xlab="時間", ylab="最小値") plot(gdata4,type="l", xlab="時間", ylab="最大値") 4

5

解析プログラム２（測定データのバイスペクトル・自己相関） dev.off() par(mfrow=c(2,2)) data11 <- read.table("C:/sample/sample-0001_01.csv", skip=6, s ep=",", nrows=14000) plot(data11$V2, main="音圧測定データ ｃｈＡ")#画面表示 a <- spectrum(data11$V2,method="ar")#パワースペクトル bispec(data11$V2)#バイスペクトル autcor(data11$V2)#自己相関 png(file="C:/sample/sample-0001_01s.png")#ファイル保存 par(mfrow=c(2,2)) plot(data11$V2, main="音圧測定データ ｃｈＡ") a <- spectrum(data11$V2,method="ar") bispec(data11$V2) autcor(data11$V2) 6

7

サンプリング時間の表示方法（解析結果 周波数０．５の値） 入力 data11 <- read.table("C:/2011/201101.csv", skip=0, sep="," , nrows=10 ) data11 応答（パソコンの画面表示） V1 V2 V3 1 Time Channel A Channel B 2 (ms) (V) (V) 3 0.00000000 0.33310950 -0.07290872 4 0.00256000 0.07199316 0.03616443 グラフ ０．５＝１９５ｋＨｚ 5 0.00512000 -0.16211430 -0.00927763 6 0.00768000 0.06299020 -0.09109775 7 0.01024000 0.05398724 0.13611260 8 0.01280000 -0.16211430 0.15430160 9 0.01536000 0.10800500 -0.22745450 10 0.01792000 0.27008880 -0.23654900 > 1 秒／0.00256000ms ＝ 390．625ｋHz 390／2＝ 195ｋHz data11 <- read.table("C:/2011/2011103.csv", skip=0, sep="," , nrows=10 ) data11 V1 V2 V3 1 Time Channel A Channel B 2 (ms) (V) (V) 3 0.00000000 -0.03604236 -0.11838130 4 0.00064000 -0.04504532 -0.06381420 グラフ ０．５＝７８1ｋＨｚ 5 0.00128000 -0.05404828 -0.02746666 6 0.00192000 -0.07205420 0.02706992 7 0.00256000 -0.05404828 0.09067049 8 0.00320000 -0.03604236 0.10885950 9 0.00384000 -0.01803644 0.09976501 10 0.00448000 -0.03604236 0.07248146 1 秒／0.00064000ms ＝ 1562．5ｋHz 1562／2＝ 781ｋHz 注：解析には下記ツールを利用します 注：OML(Open Market License) https://www.ism.ac.jp/ismlib/jpn/ismlib/license.html 注：TIMSAC(TIMe Series Analysis and Control program) https://jasp.ism.ac.jp/ism/timsac/ 注：「R」フリーな統計処理言語かつ環境 https://cran.ism.ac.jp/ 8

＜サンプリング時間の設定は自動調整されます＞ サンプリング時間 解析グラフ 0．5 の周波数（最大解析周波数） 1 秒／1s ＝ 1Hz 1／2＝ 0.5Hz ・・・ 1 秒／0.02ms ＝ 50ｋHz 50／2＝ 25ｋHz 1 秒／0.01ms ＝ 100ｋHz 100／2＝ 50ｋHz 1 秒／0.0050048ms ＝ 200ｋHz 200／2＝ 100ｋHz 1 秒／0.0020032ms ＝ 500ｋHz 500／2＝ 250ｋHz 1 秒／0.0010048ms ＝ 995ｋHz 995／2＝ 497ｋHz 1 秒／0.0005056ms ＝ 1977ｋHz 1977／2＝ 988ｋHz 1 秒／0.0002048ms ＝ 4882ｋHz 4882／2＝ 2441ｋHz（2．4MHｚ） 1 秒／0.0001024ms ＝ 9765ｋHz 9765／2＝ 4882ｋHz（4．8MHｚ） 1 秒／0.0000512ms ＝ 19531ｋHz 19531／2＝ 9765ｋHz（9．7MHｚ） 1 秒／0.0000256ms ＝ 39062ｋHz 39062／2＝ 19531ｋHz（20MHｚ） 1 秒／0.0000128ms ＝ 78125ｋHz 78125／2＝ 39062ｋHz（39MHｚ） ・・・ 1 秒／0.016μs ＝ 62.5MHz 62.5／2＝ 31.25MHｚ 1 秒／0.008μs ＝ 125MHz 125／2＝ 62.5MHｚ 1 秒／0.004μs ＝ 250MHz 250／2＝ 125MHｚ 1 秒／0.002μs ＝ 500MHz 500／2＝ 250MHｚ ・・・ 9

音圧レベルの表示 入力 data11 <- read.table("C:/20111/2011102.csv", skip=6, sep=",") mean(data11$V2) mean(data11$V3) var(data11$V2) var(data11$V3) range(data11$V2) range(data11$V3) 応答（パソコンの画面表示） > data11 <- read.table("C:/20111022w/20111022-0412.csv", skip=6, sep=",") > mean(data11$V2) CH1 の平均値 [1] -0.001047526 > mean(data11$V3) CH2 の平均値 [1] 3.430622e-05 > var(data11$V2) CH1 の分散値 [1] 0.009286384 > var(data11$V3) CH２の分散値 [1] 0.001448241 > range(data11$V2) CH1 の最小・最大値 [1] -0.4412366 0.4141362 > range(data11$V3) CH２の最小・最大値 [1] -0.1547288 0.1361126 > 注意 統計処理を行うために、測定値が自動的に、規格化（正規化）されています バイスペクトルについて理解が深まるまでは 最大・最小値、分散値、平均値 を利用することを推奨します 絶対値としての音圧は、測定データのグラフから読み取ってください その値に対する平均や分散を上記の処理で推定して利用します 10

具体例 data11 <- read.table("C:/20191220/191220-0018/20191220-0018_15.csv", skip=6, sep=",") mean(data11$V2) mean(data11$V3) var(data11$V2) var(data11$V3) range(data11$V2) range(data11$V3) グラフ青 音圧レベル ３６０ｍＶ グラフ赤 音圧レベル １６ｍＶ 11

音圧・周波数の変化を目視確認する 12

音圧・周波数の変化を音圧データの解析で確認する 13

様々な条件で、繰り返し確認すると、 音圧データのグラフを見ることで解析結果は推測できます 以上 14