超音波技術（R言語）－－超音波の音圧測定データ解析手順－－

超音波技術（R言語） 超音波の音圧測定データ解析手順 音圧測定データ 音圧解析結果 超音波システム研究所

音圧測定

１． 準備 解析用データの確認 測定データ（PSDATA ファイル）から 解析用の Microsoft Office Excel CSV ファイル がホルダーにあることを確認してください ファイルがない場合は、CSVファイルに保存してください 例 C:/us-data2/20100101-0004/20100101-0004_01.csv 上記のサンプルデータを希望される方はメールで連絡してください

２． 解析ソフトの立ち上げ ダブルクリックして立ち上げる ＞｜ 左記のようなプロンプト表示が行われます エラー表示が行われた場合には、 その他の作業ファイル・・・を終了してから もう一度立ち上げてください

３． 解析ソフトの読み込み ３－１：パッケージ－＞ ローカルにあるｚｉｐファイルからの・・・・ ３－２：デスクトップのｚｉｐファイル ＴＩＭＳＡＣ1.2.1 を選択する

３－３：読み込み画面の確認 ３－４：パッケージの読み込み パッケージ－＞パッケージ読み込み・・・ 選択する

３－５：ｔｉｍｓａｃの選択 ３－６：「ＯＫ」選択により読み込み 以上で解析準備完了です

参考 バイスペクトル バイスペクトルは以下のように 周波数 f1、f 2、f1 + f 2のスペクトルの積で表すことができる。 B( f1 , f 2 ) = X ( f1 )Y( f 2 )Z ( f1 + f 2 ) 主要周波数が f1であるとき、 f1 + f1 = f 2、f1 + f 2 = f3で表される f 2、f3という周波数成分 が存在すればバイスペクトルは値をもつ。 これは主要周波数 f1の整数倍の周波数成分を持つこと と同等であるので、バイスペクトルを評価することにより、 高調波の存在を評価できる。 詳しい説明は専門書・・・を読んで確認してください エクセルファイルのデータ列

解析 1（結果をｐｎｇファイルで保存） コマンドファイル：：解析 1（結果をｐｎｇファイルで保存）.txt data11 <- read.table("C:/us-data2/20100101-0020/20100101-0020_01.csv", skip=6, sep=",", nrows=6000) png(file="C:/us-data2/20100101-0020/20100101-0020_01.png") plot(data11$V2, main="音圧測定データ ｃｈＡ") dev.off() data11 <- read.table("C:/us-data2/20100101-0020/20100101-0020_01.csv", skip=6, sep=",", nrows=6000) png(file="C:/us-data2/20100101-0020/20100101-0020sp0001_01.png") a <- spectrum(data11$V2,method="ar") plot(a, sub="パワースペクトル") dev.off() data11 <- read.table("C:/us-data2/20100101-0020/20100101-0020_01.csv", skip=6, sep=",", nrows=6000) png(file="C:/us-data2/20100101-0020/20100101-0020bi0001_01.png") bispec(data11$V2) dev.off() data11 <- read.table("C:/us-data2/20100101-0020/20100101-0020_01.csv", skip=6, sep=",", nrows=6000) png(file="C:/us-data2/20100101-0020/20100101-0020au0001_01.png") autcor(data11$V2) dev.off() ・・・・・ ・・・・・ 実行結果（C:/us-data2/20100101-0020/）

解析 2（結果を画面にグラフで表示） コマンドファイル：：解析 2（結果を画面にグラフで表示）.txt par(mfrow=c(2,4)) data11 <- read.table("C:/us-data2/20100101-0020/20100101-0020_01.csv", skip=6, sep=",", nrows=3000) plot(data11$V2, main="音圧測定データ ｃｈＡ") data11 <- read.table("C:/us-data2/20100101-0020/20100101-0020_01.csv", skip=6, sep=",", nrows=3000) a <- spectrum(data11$V2,method="ar") data11 <- read.table("C:/us-data2/20100101-0020/20100101-0020_01.csv", skip=6, sep=",", nrows=3000) bispec(data11$V2) data11 <- read.table("C:/us-data2/20100101-0020/20100101-0020_01.csv", skip=6, sep=",", nrows=3000) autcor(data11$V2) ・・・ ・・・ 実行結果（C:/us-data2/20100101-0020/20100101-0020_01.csv C:/us-data2/20100101-0020/20100101-0020_02.csv）

解析 3（バイスペクトルを画面にグラフで表示） コマンドファイル：：解析 3（バイスペクトルを画面にグラフで表示）.txt par(mfrow=c(2,4)) data11 <- read.table("C:/us-data2/20100101-0020/20100101-0020_01.csv", skip=6, sep=",", nrows=3000) bispec(data11$V2) data11 <- read.table("C:/us-data2/20100101-0020/20100101-0020_02.csv", skip=6, sep=",", nrows=3000) bispec(data11$V2) data11 <- read.table("C:/us-data2/20100101-0020/20100101-0020_03.csv", skip=6, sep=",", nrows=3000) bispec(data11$V2) data11 <- read.table("C:/us-data2/20100101-0020/20100101-0020_04.csv", skip=6, sep=",", nrows=3000) bispec(data11$V2) data11 <- read.table("C:/us-data2/20100101-0020/20100101-0020_05.csv", skip=6, sep=",", nrows=3000) bispec(data11$V2) ・・・ ・・・ 実行結果（C:/us-data2/20100101-0020/20100101-0020_01.csv ・・・ C:/us-data2/20100101-0020/20100101-0020_32.csv）

解析コマンド dev.off() par(mfrow=c(2,2)) ： 2行 2 列のグラフ表示設定 data11 <- read.table("C:/20191220/20191220-0030/20191220-0030_12.csv", skip=6, sep=",", nrows=6000) plot(data11$V2) data11 <- read.table("C:/20191220/20191220-0030/20191220-0030_12.csv", skip=6, sep=",", nrows=6000) spectrum(data11$V2,method="ar") data11 <- read.table("C:/20191220/20191220-0030/20191220-0030_12.csv", skip=6, sep=",", nrows=6000) bispec(data11$V2) data11 <- read.table("C:/20191220/20191220-0030/20191220-0030_12.csv", skip=6, sep=",", nrows=6000) autcor(data11$V2)

dev.off() 解説 終了 plot(data11$V2) 解説 data11 の 2番目のデータ列（１ｃｈの測定データ）に対して プロット（音圧測定データのグラフ作成）を行う spectrum(data11$V2,method="ar") 解説 data11の 2 番目のデータ列（１ｃｈの測定データ）に対して ＡＲ（自己回帰）モデルによるスペクトル解析を行う bispec(data11$V2) 解説 data11の 2 番目のデータ列（１ｃｈの測定データ）に対して バイスペクトル解析を行う autcor(data11$V2) 解説 data11の 2 番目のデータ列（１ｃｈの測定データ）に対して 自己相関の解析を行う dev.off()

参考 １）以下のようにテキストデータをコピーしてＲの画面にペーストすると １ｃｈと２ｃｈのデータ比較ができます dev.off() par(mfrow=c(4,2)) data11 <- read.table("C:/20191220/20191220-0022/20191220-0022_12.csv", skip=6, sep=",", nrows=6000) plot(data11$V2) data11 <- read.table("C:/20191220/20191220-0022/20191220-0022_12.csv", skip=6, sep=",", nrows=6000) spectrum(data11$V2,method="ar") data11 <- read.table("C:/20191220/20191220-0022/20191220-0022_12.csv", skip=6, sep=",", nrows=6000) bispec(data11$V2) data11 <- read.table("C:/20191220/20191220-0022/20191220-0022_12.csv", skip=6, sep=",", nrows=6000) autcor(data11$V2) data11 <- read.table("C:/20191220/20191220-0022/20191220-0022_12.csv", skip=6, sep=",", nrows=6000) plot(data11$V3) data11 <- read.table("C:/20191220/20191220-0022/20191220-0022_12.csv", skip=6, sep=",", nrows=6000) spectrum(data11$V3,method="ar") data11 <- read.table("C:/20191220/20191220-0022/20191220-0022_12.csv", skip=6, sep=",", nrows=6000) bispec(data11$V3) data11 <- read.table("C:/20191220/20191220-0022/20191220-0022_12.csv", skip=6, sep=",", nrows=6000) autcor(data11$V3)

２）以下のようにテキストデータをコピーしてＲの画面にペーストすると すべての解析を連続的に行います dev.off() par(mfrow=c(2,2)) data11 <- read.table("C:/20191220/20191220-0022/20191220-0022_12.csv", skip=6, sep=",", nrows=6000) plot(data11$V2) data11 <- read.table("C:/20191220/20191220-0022/20191220-0022_12.csv", skip=6, sep=",", nrows=6000) spectrum(data11$V2,method="ar") data11 <- read.table("C:/20191220/20191220-0022/20191220-0022_12.csv", skip=6, sep=",", nrows=6000) bispec(data11$V2) data11 <- read.table("C:/20191220/20191220-0022/20191220-0022_12.csv", skip=6, sep=",", nrows=6000) autcor(data11$V2) ・・・・・・・・ ・・・・・・・・ data11 <- read.table("C:/20191220/20191220-0022/20191220-0022_24.csv", skip=6, sep=",", nrows=6000) plot(data11$V2) data11 <- read.table("C:/20191220/20191220-0022/20191220-0022_24.csv", skip=6, sep=",", nrows=6000) spectrum(data11$V2,method="ar") data11 <- read.table("C:/20191220/20191220-0022/20191220-0022_24.csv", skip=6, sep=",", nrows=6000) bispec(data11$V2) data11 <- read.table("C:/20191220/20191220-0022/20191220-0022_24.csv", skip=6, sep=",", nrows=6000) autcor(data11$V2)

３）以下のようにテキストデータをコピーしてＲの画面にペーストすると 各種データの関係をグラフで確認できます 下記をコピー＆ペースト setwd("D:/R-data") df <- read.csv("20170905-0054/20170905-0054_03.csv",skip=2,header=T,row.names=1) ndf <- na.omit(df)#欠損(NA)を除く ndf daBbb <- as.numeric(ndf[,1])#ndf[,1]を数字に変換 daCcc <- as.numeric(ndf[,2])#ndf[,2]を数字に変換 A1 <- ts(daBbb) B1 <- ts(daCcc) df <- read.csv("20170905-0054/20170905-0054_04.csv",skip=2,header=T,row.names=1) ndf <- na.omit(df)#欠損(NA)を除く ndf daBbb <- as.numeric(ndf[,1])#ndf[,1]を数字に変換 daCcc <- as.numeric(ndf[,2])#ndf[,2]を数字に変換 A2 <- ts(daBbb) B2 <- ts(daCcc) df <- read.csv("20170905-0054/20170905-0054_05.csv",skip=2,header=T,row.names=1) ndf <- na.omit(df)#欠損(NA)を除く ndf daBbb <- as.numeric(ndf[,1])#ndf[,1]を数字に変換 daCcc <- as.numeric(ndf[,2])#ndf[,2]を数字に変換 A3 <- ts(daBbb) B3 <- ts(daCcc) df <- read.csv("20170905-0054/20170905-0054_06.csv",skip=2,header=T,row.names=1) ndf <- na.omit(df)#欠損(NA)を除く ndf daBbb <- as.numeric(ndf[,1])#ndf[,1]を数字に変換 daCcc <- as.numeric(ndf[,2])#ndf[,2]を数字に変換 A4 <- ts(daBbb) B4 <- ts(daCcc)

par(mfrow=c(4,2)) ts.plot(A1,main="A1") ts.plot(B1,main="B1") ts.plot(A2,main="A2") ts.plot(B2,main="B2") ts.plot(A3,main="A3") ts.plot(B3,main="B3") ts.plot(A4,main="A4") ts.plot(B4,main="B4") 下記をコピー＆ペースト par(mfrow=c(2,2)) plot(A1,B1) plot(A2,B2) plot(A3,B3) plot(A4,B4)

実施例：バイスペクトルの変化を確認する dev.off() par(mfrow=c(3,2)) data11 <- read.table("C:/20191220/20191220-0022/20191220-0022_13.csv", skip=6, sep=",", nrows=6000) bispec(data11$V2) data11 <- read.table("C:/20191220/20191220-0022/20191220-0022_14.csv", skip=6, sep=",", nrows=6000) bispec(data11$V2) data11 <- read.table("C:/20191220/20191220-0022/20191220-0022_15.csv", skip=6, sep=",", nrows=6000) bispec(data11$V2) data11 <- read.table("C:/20191220/20191220-0022/20191220-0022_15.csv", skip=6, sep=",", nrows=6000) bispec(data11$V2) data11 <- read.table("C:/20191220/20191220-0022/20191220-0022_17.csv", skip=6, sep=",", nrows=6000) bispec(data11$V2) data11 <- read.table("C:/20191220/20191220-0022/20191220-0022_18.csv", skip=6, sep=",", nrows=6000) bispec(data11$V2)

実施例：自己相関の変化を確認する dev.off() par(mfrow=c(3,2)) data11 <- read.table("C:/20191220/20191220-0022/20191220-0022_13.csv", skip=6, sep=",", nrows=6000) autcor(data11$V2) data11 <- read.table("C:/20191220/20191220-0022/20191220-0022_14.csv", skip=6, sep=",", nrows=6000) autcor(data11$V2) data11 <- read.table("C:/20191220/20191220-0022/20191220-0022_15.csv", skip=6, sep=",", nrows=6000) autcor(data11$V2) data11 <- read.table("C:/20191220/20191220-0022/20191220-0022_15.csv", skip=6, sep=",", nrows=6000) autcor(data11$V2) data11 <- read.table("C:/20191220/20191220-0022/20191220-0022_17.csv", skip=6, sep=",", nrows=6000) autcor(data11$V2) data11 <- read.table("C:/20191220/20191220-0022/20191220-0022_18.csv", skip=6, sep=",", nrows=6000) autcor(data11$V2)

解析の詳細・解析結果の解釈・・・については 以下の参考書籍・・・の専門書を読んでください 参考書籍 １：統計数理 １）叩いて超音波で見る―非線形効果を利用した計測 佐藤 拓宋 (著) 出版社: コロナ社 (1995/06) ２）電気系の確率と統計 佐藤 拓宋 (著) 出版社: 森北出版 (1971/01) ３）不規則信号論と動特性推定 宮川 洋 (著), 佐藤拓宋 (著), 茅 陽一 (著) 出版社: コロナ社 (1969) ４）赤池情報量規準 AIC―モデリング・予測・知識発見 赤池 弘次 (著), 室田 一雄 (編さん), 土谷 隆 (編さん) 出版社: 共立出版 (2007/07) ５）ダイナミックシステムの統計的解析と制御 赤池 弘次 (著), 中川 東一郎 (著) 出版社: サイエンス社(1972) 参考資料 統計的な考え方を利用した超音波 http://ultrasonic-labo.com/?p=12202 超音波の洗浄・攪拌・加工に関する「論理モデル」 http://ultrasonic-labo.com/?p=3963 物の動きを読む http://ultrasonic-labo.com/?p=1074 超音波技術：多変量自己回帰モデルによるフィードバック解析 http://ultrasonic-labo.com/?p=15785