Document - 松本充弘の個人的ページ

機械理工学・マイクロエンジニアリング・航空宇宙工学専攻
熱物理工学 2016（松本担当分）講義資料
松本充弘
第４回 2016/06/20
[email protected]
*レポートの提出状況は別紙の受講者リストで確認してください．間違いがあれば松本まで．
前回のレポート課題について
課題楽曲のスペクトル解析をおこなう．
○レポートから（三上慎司さん）
ベートーベンの「交響曲第九番」（通称）「第九」の前奏と中間
奏のスペクトル解析を行った。データ数は
63 であり，八分
音符などの半端な長さの音符は，切り捨てもしくは切り上げを
行った。得られたデータは下図のようになった。「第九」スペクト
ルはおおよそ 1 / f 2 に近いように見える。
ちなみに，一つのデータあたりの長さを八分音符の長さにす
れば量はほぼ 2 倍になるが，より正確に表現できたのではない
かと思い，プログラムを実行した。得られた結果が下図である。
○レポートから（山本隆正さん）
一部を抜粋
世界に一つだけの花のサビの一部について，パワースペクト
ルを求めてグラフを描いた．ちょうど 1/f と 1/f^2 の間に点を取っ
ているように見える．
○レポートから（西岡寿朗さん）
一部を抜粋
題材として音符が簡単そうな「君が代」を用いました．メロディ
の数値化は，ドレミファ…を 1,2,3,4…に置き換えて行った．具
体的には，22112233 55332222 33556656 99776655
3355666699889999 33556655 33352222 66889999 88996655
66532222 であり，データ数 ndata は 88 個であった．
パワースペクトルは周波数相当である i の中央 i =44 まわり
に偶対称であることが確認できる．つまり片側 i=1, 2, ... , 44 (=
ndata/2) までのスペクトルの情報で，このグラフのすべてのス
ペクトルの情報がわかるということである．これはサンプリング
定理の一例であると考えられる．サンプリング定理とは，「波形
の最大周波数の 2 倍以上の周波数で標本化すれば完全再
構成」というもので，別の表現をすると，「サンプリング周波数が
fs であるとき、離散化して表現可能である周波数の上限は fs/2
に限られる」である．この fs/2 はナイキスト周波数と呼ばれる．
（中略）次にパワースペクトルの解析を行う．先程の議論から，
熱物理工学講義資料
1
中央 i = 44 より左側のみを考える．（中略）両対数表示された点に対して線形回帰を行ったところ，
Power = 𝑒 10.179 𝑖 −2.0093 という結果を得た．Fig. 2 に淡青で示す．その相関係数は 0.833 であった．こ
れよりパワースペクトルが i の二乗に逆比例していると考えられる．
（松本の補足）スペクトルが「左右対称」になることの略証
等間隔でサンプリングされた N 個のデータ Fm があるとする．Fourier 係数の定義は
N 1
Gn   Fm exp i n m,
m 0
n 
2 n
N
そこで
N 1
G N n   Fm exp i N n m
m 0
N 1
 2 ( N  n ) 
  Fm exp i
m
N


m 0
N 1
 2n 
  Fm exp 2miexp   i
m
N 

m 0
N 1
  Fm exp  in m
m 0
 Gn
よって Power (n)  Gn  Gn  GN n  GN n  Power ( N  n) が成り立つ．
（証明終）
○レポートから（原田浩行さん）
一部を抜粋
童謡「グリーングリーン」のスペクトルはおおよそ 1/f に近いよ
うに見える。
○レポートから（徐海元さん）
「Mary had a little lamb」のスペクトルはおおよそ 1/f^2
に近いように見える．
「森のくまさん」のスペクトルはおおよそ 1/f に近いように見える．
○レポートから（前田篤弥さん）
子守唄は心地よさを伴うものであり，1/f ゆらぎと関
連があると考え，江戸子守歌についてスペクトル解析
を行った．楽譜は Wikipedia に掲載されたものを参考
熱物理工学講義資料
2
にした．江戸子守歌のパワースペクトルは比較的ブラウンノイズに近いゆらぎであることが確認
できた．また，得られたデータは 31 個であったが，15 個目のデータと 16 個目のデータを境界と
して対称なデータであり，調査するデータは半分で済
むことがわかった．
○レポートから（久保田奈緒さん）
一部を抜粋
キラキラ星と、カエルの歌をプロットしています。カ
エルの歌がグラフ上で下(？)に振れているように見え
たので、10000/(x^3)の線も付け足しました。
○レポートから（上原水優さん）
一部を抜粋
「サザエさん」と「カエルの歌」の 2 曲について, フーリエ解
析を行った. （中略）両対数グラフにおいて, 𝐹 ∝ 𝑓 −1 と
𝐹 ∝ 𝑓 −2 の 2 つについて, 最小二乗法により回帰直線を
描いた. 「サザエさん」は𝐹 ∝ 𝑓 −1 については 40.56,
𝐹 ∝ 𝑓 −2 については 62.35 となり, 「カエルの歌」は
𝐹 ∝ 𝑓 −1 については 19.32, 𝐹 ∝ 𝑓 −2 については 29.81 と
なった. このことから, 「サザエさん」も「カエルの歌」も, 𝐹 ∝ 𝑓 −1 の傾向が強いと考えられる.
松本のコメント：対数プロット，なぜか，かなり水平に近く見えますねぇ．最小二乗法でフィッティ
ングしてみるのはいいことですが，上に述べた「サンプリング定理」から，最大周波数の半分ま
でのデータのみ使うべきですね．あと，各データ点にどのような重みをつけるべきかというのも，
本当は重要そうですね．
○レポートから（広田翔さん）
一部を抜粋
学校のチャイムに対してドレミファソラシを順に１２３４５６７と
いう数値データに関連付けてスペクトル分解を行った。1/f^2
には近くない。
ちなみに私は以前、自分のツイッターについてスペクトル分
解を行った。すなわち、各日ごとにツイート数を抽出してデ
ータ化し、ツイート数の多少の周期が見られないか調査した。
その結果、周期は見られなかった。つまらない。
0
1
3
1
松本のコメント：対数プロットをみると，周波数が 10 →10 でパワーが２桁ほど（10 →10 ）落ちて
2
いるように見えるので，まあ 1/f に近いんでしょうね．ところで，「ツイート数のスペクトル」，
たとえば週末に増加するなどの理由で，７日周期などは存在してもいい気がしますが，どうだっ
たでしょうか．
○レポートから（陳聰さん）
一部を抜粋
とても私を落ち着ける「青石巷」と言うピアノの曲にした．メロディの前 14
熱物理工学講義資料
3
節を 7 音階主義で数値データにして，パワーベクトルをプロットした．おおよそ 1/f に近いように見える．
○レポートから（角木亮介さん）
一部を抜粋
パッヘルベルのカノンについてパワ
ースペクトル解析を行った．（中略）
図 1 を見ると結果は左右対称になっ
ており，これはサンプリング定理に対
応している．また，同じ結果を両対数
グラフに直した図 2 を見ると，1/f^2 とおおよそ近いように見える．カノンは 1/f 揺らぎを持っていると紹介さ
れていたので，意外な結果だった．
○レポートから（伊藤圭人さん）
一部を抜粋
京都大学学歌の前半部分を選んだ。メロディの数値化は 12 音階で
行った。16 分音符を単位とした。得られたスペクトルは f^-2 に近い。
○レポートから（白水仁さん）
一部を抜粋
「ミッキーマウス・マーチ」のスペクトルはおおよそ 1/f ノイズに近いと言
えます。
○レポートから（寺本達哉さん）
一部を抜粋
某コンビニ（Ｆ社）の入店時の音楽でやってみました。やはり、入店時
に流す効果はあるようで、1/f に近い結果が得られました。
○レポートから（樋口智哉さん）
一部を抜粋
日本では「蛍の光」として知られているスコットランド民謡の「 Auld
Lang Syne」を取り上げる。（中略） 1/𝑓に近いように見える。
○レポートから（近藤諒さん）
一部を抜粋
10000000
1000000
matlab を用いて中島みゆきのファイト!を周波数領域で表示した
100000
クトルと，採譜したデータのパワースペクトルは，違う概念で
あることに注意してください．詳しくは解説時に．
Power
10000
松本のコメント：音楽の生データ（例えば midi とか）のパワースペ
10000/i
1000
蛍の光
10
1
1
0.1
○レポートから（瀧沢賢さん）
10000/(i*i)
100
10
100
1000
i(~frequency)
一部を抜粋
サンプルプログラムで，fscanf のフォーマットは %f ではな
く %lf ではないでしょうか．%f ではうまく読めませんでした．
松本のコメント：あわてて私のプログラムを見直したのですが，どこ
も正しく %lf になっているようですが... 倍精度実数型の変数に読み込むので, 当然 %lf です．
熱物理工学講義資料
4
第４回レポート課題
出題日：6 月 20 日（月）
提出期限：6 月 25 日（土）
提出方法：電子メールで，[email protected]
へ．
課題今回はもちろんフラクタル解析の練習です．これまで，
 2015 年：雲の形
 2014 年：地形図
 2013 年：細菌集団の形状
 2012 年：毛細血管
 2011 年：フィヨルドの海岸線
とさまざまなものを取り上げてきました．今年は，ある雑誌（Nature, vol 532, April 14, p.166）に載
っていたがん組織の形状を解析してみようと思います．どんな条件で，なぜ，このような形をとるの
かについて私は全く知りません．単に，フラクタル解析のテストデータとして取り上げるだけです．
この記事の冒頭にあった “がん組織の顕微鏡像” を，まず bitmap ファイルとして保存しました．これを，
付録に示すプログラムで，適当なしきい値 threshold で２値化し，白黒のデータとしました．「見た
目」と同じような形になるしきい値をうまく選ぶのは難しく，元の画像とはやや印象が違うかもしれ
ませんが，気にしないでください．
→
論文の 1 ページ目と，抽出した画像を 2 値化したもの
1.
私のホームページ
http://www.mitsuhiromatsumoto.mech.kyoto-u.ac.jp/
から，２値化画像データ tumour_truncated_mono.bmp をダウンロードしてください．
2.
解析用プログラム fractal_analysis.cpp もダウンロードしてください．これは，講義ノート 4.5.3 節で紹介
した，「図形をタイルで覆ってその枚数を数える」という発想で図形のフラクタル次元を解析するプログラ
ムです．このプログラムを完成させて，この画像データのフラクタル次元を見積もってください．
熱物理工学講義資料
5
（ヒント） for ループから「脱出」するのには，break あるいは goto が使えます．使い方は，web ででも
調べてみてください．
fractal_analysis.cpp
#include <stdio.h>
#include <math.h>
#define DFILE "tumour_truncated_mono" // 入力用画像データファイル名（拡張子は .bmp を自動的に付加する）
#define DMAX 65536
// Maximum header size
#define IMAX
4800
// Maximum image size
int Width, Height;
int Image[IMAX][IMAX];
int read_image(char* INPUT_FILE ) //=====================================================
{
FILE *fin;
unsigned char cdata[DMAX];
unsigned int Offset, HeaderSize;
unsigned short int NBit, dummy;
unsigned int Cmprs;
unsigned int idata;
int x,y;
fin=fopen(INPUT_FILE,"r");
if (fin==NULL) {
printf("Error: no input file¥n");
return -1;
}
// Bitmap File Header
fread(cdata,10,1,fin);
printf("File Type = %c%c¥n",cdata[0],cdata[1]);
if (cdata[0]!='B' || cdata[1]!='M') {
printf("Error: file is not BMP¥n");
fclose(fin);
return -1;
}
fread(&Offset,4,1,fin);
// Bitmap Info Header
fread(&HeaderSize,4,1,fin);
printf("Header Size = %d byte¥n",HeaderSize);
printf("Offset = %d byte¥n",Offset);
fread(&Width,4,1,fin);
fread(&Height,4,1,fin);
fread(&dummy,2,1,fin);
fread(&NBit,2,1,fin);
fread(&Cmprs,4,1,fin);
printf("Image Size = %d X %d (%d bit)¥n",Width,Height,NBit);
printf("Image Compression = %d¥n",Cmprs);
if (NBit!=24) {
printf("Sorry, this program works only for 24-bit image¥n");
fclose(fin);
return -1;
}
if (Cmprs!=0) {
printf("Sorry, this program does not support this compression type¥n");
fclose(fin);
return -1;
}
if ((Width>IMAX) || (Height>IMAX)) {
printf("Sorry, image is too large¥n");
fclose(fin);
return -1;
}
fclose(fin);
fin=fopen(INPUT_FILE,"rb");
fread(cdata,Offset,1,fin);
// Read Image
for (y=0;y<Height;y++) {
熱物理工学講義資料
6
for (x=0;x<Width;x++) {
if (fread(&idata, 3, 1, fin)<1) {
printf("Error: data end while reading¥n");
return -1;
}
Image[y][x]=idata;
// 黒 =0, 白 !=0
} }
fclose(fin);
return 0;
}
int fractal(char* OUTPUT_FILE) //=============================
{
int x, y, xi, yi;
int resolution;
// Tile Size
int sum;
// Number of Tiles
int res_max;
FILE *fout;
fout=fopen(OUTPUT_FILE,"w");
if (Height<Width) res_max=Height/2;
else
res_max=Width/2;
for (resolution=2; resolution<=res_max; resolution*=2) { // タイルサイズを２倍ずつ大きくする
sum=0; // 「条件」に合致するタイルの枚数
// -------------------------------------------------------------------for (x = 0; x < Width/resolution; x++) {
for (y = 0; y < Height/resolution; y++) { // タイルを１枚ずつチェックしていく
・・・
この間のコードを作成する
・・・
} }
// -------------------------------------------------------------------fprintf(fout,"%5d %10d¥n", resolution, sum);
printf("%5d %10d¥n", resolution, sum);
// タイルの枚数を出力
}
fclose(fout);
return 0;
}
int main( ) //=================================================
{
int ret;
char cfile[100];
FILE *fout;
sprintf(cfile,"%s.bmp",DFILE);
ret=read_image(cfile);
if (ret!=0) return -1;
printf("Analyzing %s¥n",cfile);
sprintf(cfile,"fractal_analysis.dat");
fractal(cfile);
// 出力用データファイル
return 0;
}
本日のレポート課題は以上です．次回は６月２７日（月）．
熱物理工学講義資料
7
（参考１）
ビットマップ形式（正式には，windows bitmap image format，ファイル拡張子は通常 .bmp）は，
画像ファイル形式の１つで，jpeg や tiff など他の画像形式に比べて容量が大きくなる傾向があるが，
無圧縮のときには，比較的単純なファイル構造をしており，初心者にとって扱いが容易である．

14 byte のファイルヘッダ，先頭は必ず “BM” である

（多くの場合）40 byte の情報ヘッダ

カラーパレット

ビットマップデータ（すなわち画像本体）
からなると規定されている．通常，１つのカラーが 3 byte (RGB 各 1 byte) あるいは 4 byte （予備 1
byte を含む）で表現される．詳細な規定については，ウェブ情報など（例えば，Wikipedia の “Windows
bitmap”）を参照されたい．
（参考２）
画像ファイルはテキストエディタでその中身を見たり編集したりすることを考えていないので，一般
にバイナリファイルである．したがって，画像の閲覧には，特殊なソフトウェア（ビューワー）が必
要なのである．C プログラム中でバイナリファイルを扱うには，バイナリモードでファイルを開く必
要があり，また読み込みも少し特殊である．（テキストデータの読み込みに使う fscanf( ) ではう
まくいかない．
）
今回のプログラムでは，以下の関数を使っている：

fopen(”ファイル名”,”rb”) ファイルをバイナリモードで開く

fread(変数ポインタ，サイズ，個数，ファイルポインタ) データを読み込む
（参考３）
今回の教材に使った画像は，論文の pdf に含まれていた 713×1120 ピクセルのカラー画像を，1 行が
4 バイトの倍数になるようにサイズを少し調整したのち，bmp 形式で保存し，さらに，以下のプログ
ラムにより２値化したものである．ここで使った bmp 形式（カラー24bit）では，各画素が r (red), g
(green), b (blue) それぞれ８ビット（2^8=256 階調）のデータを持っており，ここでは，
(r+g+b)/3 < THRESHOLD
のときにその画素を黒（r=g=b=0）とするように設定した．THRESHOLD は試行錯誤で与えた．
二値化に使ったプログラム Program for monochromatization
#include <stdio.h>
#include <math.h>
#define DFILE "tumour_truncated"
#define THRESHOLD 134
#define DMAX
#define IMAX
// File name should be given here
// threshold value (<256) to binarize the image
65536
9600
int Width, Height;
int Image[IMAX][IMAX];
unsigned int FSize, Offset, HeaderSize;
unsigned short NBit, ISize, HRes, VRes, NColor, NIndex;
unsigned int Cmprs;
int read_image(char* INPUT_FILE ) //===============================================
熱物理工学講義資料
8
{
FILE *fin;
unsigned char cdata[DMAX];
unsigned short dummy;
unsigned int idata=0;
int x,y;
fin=fopen(INPUT_FILE,"r");
if (fin==NULL) {
printf("Error: no input file¥n");
return -1;
}
// Bitmap File Header
fread(cdata,2,1,fin);
printf("File Type = %c%c¥n",cdata[0],cdata[1]);
if (cdata[0]!='B' || cdata[1]!='M') {
printf("Error: file is not BMP¥n");
fclose(fin);
return -1;
}
fread(&FSize,4,1,fin);
fread(&dummy,2,1,fin);
fread(&dummy,2,1,fin);
fread(&Offset,4,1,fin);
// Bitmap Info Header
fread(&HeaderSize,4,1,fin);
printf("Header Size = %d byte¥n",HeaderSize);
printf("Offset = %d byte¥n",Offset);
fread(&Width,4,1,fin);
fread(&Height,4,1,fin);
fread(&dummy,2,1,fin);
fread(&NBit,2,1,fin);
fread(&Cmprs,4,1,fin);
printf("Image Size = %d X %d (%d bit)¥n",Width,Height,NBit);
printf("Image Compression = %d¥n",Cmprs);
if (NBit!=24) {
printf("Sorry, this program works only for 24-bit image¥n");
fclose(fin);
return -1;
}
if (Cmprs!=0) {
printf("Sorry, this program does not support this compression type¥n");
fclose(fin);
return -1;
}
if ((Width>IMAX) || (Height>IMAX)) {
printf("Sorry, image is too large¥n");
fclose(fin);
return -1;
}
fread(&ISize, 4,1,fin);
fread(&HRes, 4,1,fin);
fread(&VRes, 4,1,fin);
fread(&NColor, 4,1,fin);
fread(&NIndex, 4,1,fin);
fclose(fin);
fin=fopen(INPUT_FILE,"rb");
fread(cdata,Offset,1,fin);
// Read Image
for (y=0;y<Height;y++) {
//
printf("%4d ",y);
for (x=0;x<Width;x++) {
if (fread(&idata, 3, 1, fin)<1) {
printf("Error: data end while reading¥n");
return -1;}
Image[y][x]=idata;
} }
fclose(fin);
return 0;
}
int write_image(char* OUTPUT_FILE ) //==============================================
熱物理工学講義資料
9
{
FILE *fout;
unsigned char cdata[DMAX];
unsigned short dummy;
unsigned int idata;
int x,y;
unsigned int ZERO=0, ndata;
unsigned int white, black, rdata, gdata, bdata;
fout=fopen(OUTPUT_FILE,"wb");
if (fout==NULL) {
printf("Error: cannot create output file¥n");
return -1;
}
// Bitmap File Header
fwrite("BM", 2,1,fout);
fwrite(&FSize, 4,1,fout);
fwrite(&ZERO, 2,1,fout);
fwrite(&ZERO, 2,1,fout);
fwrite(&Offset, 4,1,fout);
// Bitmap Info Header
fwrite(&HeaderSize, 4,1,fout);
fwrite(&Width, 4,1,fout);
fwrite(&Height, 4,1,fout);
ndata=1;
fwrite(&ndata, 2,1,fout);
fwrite(&NBit, 2,1,fout);
fwrite(&Cmprs, 4,1,fout);
fwrite(&ISize, 4,1,fout);
fwrite(&HRes, 4,1,fout);
fwrite(&VRes, 4,1,fout);
fwrite(&NColor, 4,1,fout);
fwrite(&NIndex, 4,1,fout);
// Write Image
white = 256*256*256-1;
black = 0;
for (y=0;y<Height;y++) {
//
printf("%4d ",y);
for (x=0;x<Width;x++) {
idata=Image[y][x];
rdata=idata%256;
idata/=256;
gdata=idata%256;
bdata=idata/256;
idata=(rdata+gdata+bdata)/3;
if (idata<THRESHOLD) fwrite(&black, 3, 1, fout);
else
fwrite(&white, 3, 1, fout);
} }
fclose(fout);
return 0;
}
int main( ) //=================================================
{
int ret;
char cfile[100];
FILE *fout;
sprintf(cfile,"%s.bmp",DFILE);
ret=read_image(cfile);
if (ret!=0) return -1;
printf("Transforming %s¥n",cfile);
//
return 9; // For debug only
sprintf(cfile,"%s_mono.bmp",DFILE);
write_image(cfile);
return 0;
}
熱物理工学講義資料
10

Download Report