pythonで、ホワイトノイズやピンクノイズのwavファイルを作る。

ホワイトノイズを作る

ホワイトノイズは、ランダムノイズ。
ランダムノイズは、フーリエ変換してもほぼランダムノイズ。（対数グラフで描くと、そうとも言い切れないが。）

「〜wavファイルを作る。」と書いたが以前書いたPythonで.wavファイルを出力（標準のwaveとscipy.io.wavfileの速度比較）の写しで、wavファイルを作れるので、そこら辺は省略。なので、ホワイトノイズに関しては、ランダムノイズのプロットだけ乗せとく。ブラウンもピンクも同様。

以前は、wavファイルにsinカーブを入れたが、今回は、その部分のみランダムノイズに変えたり、周波数をいじって加工する。
ランダムノイズは、numpyの機能を使って書いた。

import numpy as np
import matplotlib.pyplot as plt


stop = 10
Srate = 44100
x1=np.linspace(0, stop, Srate*stop)
A=0.5
y1 = np.random.rand(Srate * stop)-A

plt.plot(x1, y1)
plt.show()

#plt.plot(x1[:100], y1[:100])
#plt.show()

import numpy as np
import matplotlib.pyplot as plt


stop = 10
Srate = 44100
A=0.5
y1 = np.random.rand(Srate * stop)-A
f  = np.linspace(1/float(stop), 44100, Srate * stop)
y11 = np.fft.fft(y1)/float(Srate * stop /2)
x1=np.linspace(0, stop, Srate*stop)

fig,((ax1,ax12),(ax2,ax22),(ax3,ax32))=plt.subplots(3,2)

ax1.plot(x1, y11)
ax12.plot(f, np.abs(y11))
ax22.plot(f, np.real(y11))
ax32.plot(f, np.imag(y11))
ax1.set_title("Input(white noise)")
ax12.set_title("abs")
ax22.set_title("real")
ax32.set_title("imag")
#plt.xscale("log")
#plt.yscale("log")
plt.tight_layout()

plt.show()

ホワイトノイズのスペクトルも、ホワイトノイズ。

全周波数にわたって、値は同程度（の振幅のノイズ）。

これがホワイトノイズの良くないところ。

人の耳は、周波数ごとに感度が違う。高周波には敏感。

どの周波数も同程度に感じるような音がピンクノイズらしい。

ピンクノイズの前に、ブラウニアンノイズを作る。

ブラウニアンノイズも作る。

ブラウニアンノイズとかの1/f^2とか1/fノイズを作るのは難しい。

wikipediaのピンクノイズのページを見ると、「パワーが周波数に反比例する雑音のこと。」とある。

ピンクノイズを作るための方法として、voss法というのがあるらしい。
しかし、あんまり難しいのは避けたい。「新しい勉強をするのも時間的にちょっと…。」と思ったので諦めた。

そこで、思いついたやり方でやる。

ホワイトノイズの信号を使って

ホワイトノイズをfft→周波数空間で加工→ifftで元に戻す。

もしこのやり方に名前がないのであれば、「fourie (変換おじさん)法」とでもしておく。

なぜかホワイトノイズのままになってしまう

「間違いなく作ったのに、出てくる音はホワイトノイズのそれなんだよなぁ…」と思って、スペクトル調べたら、ホワイトノイズだった。

そのあと、原因がわからなくて、四苦八苦した。
原因は、スペクトルの折り返しを消してしまっていたからだった。

スペクトルは、折り返しの所で、傾きが逆になる。

1/f^0=1（無処理：ホワイトノイズ）と1/f^2（ブラウンノイズ）の間がピンクノイズらしい。

ナイキスト周波数の所で折り返すので、fmax/2の所で、傾きを逆にしなければならない。

とりあえず、ブラウンノイズから作りたい。
ブラウンノイズっぽく1／f^2で減っている線を描く。

import matplotlib.pyplot as plt
import numpy as np


stop = 10
Srate = 44100
f  = np.linspace(1/10.0, 44100, Srate * stop)
f1 = 1.0/f**2
f1 = np.r_[f1[:len(f1)/2],f1[1:len(f1)/2+1][::-1]]


plt.plot(f,f1)
plt.yscale("log")
plt.xlabel("frq. in log scale")
plt.xscale("log")
plt.ylabel("mag. in log scale")
plt.show()

logスケールで直線になった。
ただ、1Hz以下の部分が大きすぎる。
あと、後述するが、20Hz以下はいらんので、20Hzの所が1になるようにしたい。

import matplotlib.pyplot as plt
import numpy as np


stop = 10
Srate = 44100

f  = np.linspace(1/10.0, 44100, Srate * stop)
f1 = 1.0/f**2
v  = f1[20*stop]
f1 = f1/v
f1 = np.r_[f1[:len(f1)/2],f1[1:len(f1)/2+1][::-1]]


plt.plot(f,f1)
plt.yscale("log")
plt.xlabel("frq. in log scale")
plt.xscale("log")
plt.ylabel("mag. in log scale")
plt.show()

これを、ホワイトノイズにかけ合わせたい。

import numpy as np
import matplotlib.pyplot as plt


stop = 10
Srate = 44100
A=0.5
y1 = np.random.rand(Srate * stop)-A
f  = np.linspace(1/float(stop), 44100, Srate * stop)
f1 = 1.0/f**2
v  = f1[20*stop]
f1 = f1/v
f1 = np.r_[f1[:len(f1)/2],f1[1:len(f1)/2+1][::-1]]
y11 = np.abs(np.fft.fft(y1)*f1)/float(Srate * stop /2)


plt.plot(f, y11)
plt.xscale("log")
plt.yscale("log")
plt.show()

このスペクトルは、このままだと、まずい。
まずい理由は、

• 可聴域の周波数で、強度が低すぎる。
• 低周波が強すぎる。（どうせスピーカーでは再生できない周波数が強く出すぎて、ほかの音が小さい。これはこれで、ピンクノイズではあるのだろうけど。）

試しに、ifftかけると

y11 = np.fft.ifft(np.fft.fft(y1)*f1)
plt.plot(np.real(y11))
plt.show()

のっぺりしてしまっている。

スピーカーの再生できる周波数に絞る（20〜20kHz位だと思っている。）

スペクトルの、赤で囲った領域はいらないので、0で置換。
それ以外の周波数は、20 Hzの値が、1以下に収まるようにする。

import numpy as np
import matplotlib.pyplot as plt


stop = 10
Srate = 44100
A=0.5
y1 = np.random.rand(Srate * stop)-A
f  = np.linspace(1/float(stop), 44100, Srate * stop)
f1 = 1.0/f**2
v  = f1[20*stop]
f1 = f1/v
f1[:stop*20]=0
f1[stop*20000:]=0
f1 = np.r_[f1[:len(f1)/2],f1[1:len(f1)/2+1][::-1]]
f11 = f1
y11 = np.abs(np.fft.fft(y1)*f11)

plt.plot(f, y11/float(Srate*stop/2))
plt.xscale("log")
plt.yscale("log")
plt.show()

こうなった。これを、上げる。

import numpy as np
import matplotlib.pyplot as plt


stop = 10
Srate = 44100
A=0.5
y1 = np.random.rand(Srate * stop)-A
f  = np.linspace(1/float(stop), 44100, Srate * stop)
f1 = 1.0/f**2
v  = f1[20*stop]
f1 = f1/v
f1[:stop*20]=0
f1[stop*20000:]=0
f1 = np.r_[f1[:len(f1)/2],f1[1:len(f1)/2+1][::-1]]
f11 = f1
y11 = np.fft.fft(y1)*f11
y11 = y11/abs(y11[20*stop])*float(Srate*stop/2)
y11 = np.abs(y11)/float(Srate*stop/2)

plt.plot(f, y11)
plt.xscale("log")
plt.yscale("log")
plt.show()

逆変換すると、

import numpy as np
import matplotlib.pyplot as plt


stop = 10
Srate = 44100
A=0.5
x1=np.linspace(0, stop, Srate*stop)
y1 = np.random.rand(Srate * stop)-A
f  = np.linspace(1/float(stop), 44100, Srate * stop)
f1 = 1.0/f**2
v  = f1[20*stop]
f1 = f1/v
f1[:stop*20]=0
f1[stop*20000:]=0
f1 = np.r_[f1[:len(f1)/2],f1[1:len(f1)/2+1][::-1]]
f11 = f1
y11 = np.fft.fft(y1)*f11
y11 = y11/abs(y11[20*stop])*float(Srate*stop/2)
y11 = np.real(np.fft.ifft(y11))/stop/2.0

plt.plot(x1, y11)
plt.show()

 

１より大きくなっちゃった。（泣）

ピンクノイズは、

周波数的には

import matplotlib.pyplot as plt
import numpy as np


stop = 10
Srate = 44100
f  = np.linspace(1/10.0, 44100, Srate * stop)
f1 = 1.0/f
v  = f1[20*stop]
f1 = f1/v
f1 = np.r_[f1[:len(f1)/2],f1[1:len(f1)/2+1][::-1]]

plt.plot(f,f1)
plt.yscale("log")
plt.xlabel("frq. in log scale")
plt.xscale("log")
plt.ylabel("mag. in log scale")
plt.show()

できてきたノイズのスペクトルは、

import numpy as np
import matplotlib.pyplot as plt


stop = 10
Srate = 44100
A=0.5
y1 = np.random.rand(Srate * stop)-A
f  = np.linspace(1/float(stop), 44100, Srate * stop)
f1 = 1.0/f
v  = f1[20*stop]
f1 = f1/v
f1[:stop*20]=0
f1[stop*20000:]=0
f1 = np.r_[f1[:len(f1)/2],f1[1:len(f1)/2+1][::-1]]
f11 = f1
y11 = np.fft.fft(y1)*f11
y11 = y11/abs(y11[20*stop])*float(Srate*stop/2)
y11 = np.abs(y11)/float(Srate*stop/2)

plt.plot(f, y11)
plt.xscale("log")
plt.yscale("log")
plt.show()

ピンクノイズは、ブラウニアンよりも減りがなだらか。

import numpy as np
import matplotlib.pyplot as plt


stop = 10
Srate = 44100
A=0.5
x1=np.linspace(0, stop, Srate*stop)
y1 = np.random.rand(Srate * stop)-A
f  = np.linspace(1/float(stop), 44100, Srate * stop)
f1 = 1.0/f
v  = f1[20*stop]
f1 = f1/v
f1[:stop*20]=0
f1[stop*20000:]=0
f1 = np.r_[f1[:len(f1)/2],f1[1:len(f1)/2+1][::-1]]
f11 = f1
y11 = np.fft.fft(y1)*f11
y11 = y11/abs(y11[20*stop])*float(Srate*stop/2)
y11 = np.real(np.fft.ifft(y11))/stop/2.0

plt.plot(x1, y11)
plt.show()

これも１より大きくなっちゃった。

こうしてできた波形を、wavファイルに入れるのであった。。。

この後、wav作ったけど…

この後、wavファイルにした。
したけど、、、どうも、これだけじゃダメなようだ。信号の最大値が高すぎて、破裂音がする。

もう少し、全体のレベルを調節しないと。

あと、「結局聞いているのって、高周波音」というのがある。高周波がないと、全然聞こえないし。

しばらく書かないつもりだったが、せっかく思いついたので書いておく。

本屋に行った。

酒を飲みながら書く。

当然、読み直したりなんかしない。

３〜４か月は、手を抜く。

はやく、正月になんねえかなぁ。

pythonとか、matplotlibとか、コードの類は、ぱっと書けたら載せたいが、そんな余裕は、ない。

しかし、更新は続ける。

神戸の三宮に行った。

junku堂の場所を覚えた。これで、本屋には困らない。

駅前の商業施設の７F。

junku堂で、10,000円以上買ったら、

布バッグに入れてくれた。

英国屋の500円券くれた。

使わないで大切にとっておくことにする。

期限は12月だが、知らない間に過ぎていると思う。

さいたまだったら、書楽に行っていた。

あそこは駐車券しかくれない。

１Fに、カフェがあるにも関わらず！

500円券よこせ！（おれは徒歩だぞ！？）

東京駅通ってたから、丸善にかよっててもよかった気がする。

しかし、あそこは狭い。いやだ。

その点、junku堂も、本当は若干嫌。

神戸も都会すぎる。しかし、我慢する。

junku堂になかったら、amazon。神戸にいる間はこうすることにした。

続きを読む

軸対称要素の荷重のかけ方

軸対称要素の荷重のかけ方

タイトルの先頭に、「(FEMAPの)」と書いていたが消した。

画像に書いた。もはや、文字を書く気すらおきん。
たったこれだけのために、僕のお盆休みが消えた。
本当は、「ちょちょっ」とやって「ぱっ」とTwitterなんかに投稿して、
僕の満たされない承認欲求を満たそうと画策していた。
そんなこと、もったいなくってできない。

この画像は、40880文字分のPythonコードで作られている。
行数にすると、1155行。
（不必要な空白行、スペースも多いので、実質800〜900行か？）

なんで途中でやめられなかったのか、悔やまれる。ぼくのオヴォン休み。と、最終日22:18に思うのであった。

Femapで円筒座標系を使って荷重をかける時の注意点。

有料ソフトなので、画像は載せなかった。
絵がないと、分かりにくい内容だけど。

なぜFemapを使っているのか。

計算力学技術者試験の受験資格を手に入れたい。
計算力学技術者の初級（２級の受験資格になる）は、講習で手に入る（結構高い）。

その講習で、Femapを使った。
Femapは300ノード制限バージョンを使用した。

Femapをいじっていて、使ったことのない円筒座標系を使ってみようと思った。
円筒座標系の荷重でつまづいたので、書き残す。

先に結論を書いておくと、
円筒座標系の中心は、原点。任意の場所を円筒の中心にできない。

円筒座標系で荷重をかけたかった。

こうすりゃいいんじゃない？と思って、

• 座標系を「全体直交座標系」→「全体円筒座標系」に変更。
• 方向を「ベクトル」にしdR, dT, dZ(中心からの距離Radius、角度Theta, 高さ座標Z)を(1., 0., 0.)に変更。
• 荷重を入力。

これだけでできると思っていた。
（通常は、これだけでできる。）

[あらぬ方向に荷重方向の矢印が伸びている絵をここに乗せたかった。]

うまくいっていません。
明らかに、円筒の中心から伸びてきた矢印ではない。
この矢印の根元は、どこをさしているのかというと、(0, 0, 0)

問題なのは、モデルの原点。

最初に書かなかったけど、すでにモデルはできていた。
そのとき、円筒座標系のことを何も考えずに作ってしまっていた。
円筒の中心は、原点から離れた場所に作っていた。

これが原因だった。

円筒座標系とベクトルを使って荷重をかけると、原点から線が伸びる。
だから、荷重方向がおかしなことになった。

これは、絵がないとわからないでしょう。
なのに絵は載せられない（怖くて）。

解決策

1. 円筒の中心点を、原点に移動する。
2. 座標系を「全体直交座標系」→「全体円筒座標系」に変更。
3. 方向を「ベクトル」にしdR, dT, dZ(中心からの距離Radius、角度Theta, 高さ座標Z)を(1., 0., 0.)に変更。
4. 荷重を入力。

2以降は、上と同じ。

局所円筒座標系はないのか

円筒座標系の中心は、原点。任意の場所を円筒の中心にできない。
「全体円筒座標系」の中心は、絶対原点。
なら、「局所円筒座標系」とか作って、任意の場所を中心に荷重をかけられるようにしてくれたらいいのに。
あっても多分使わないけど。

• ２個以上の筒状のタンクか何かが、連結されていて
• 筒の内部で、圧力ではない荷重が放射状にかかっている
  こういうのを模擬するときに使えるかもしれないと思った。

まぁ、そんなん無くても、三角関数（円周上の座標(cosΘ, sinΘ)）を使って値を出して、ノード一点一点にその値を当てはめればできますが。

脳ドッグに行ってきた。→MRIの画像データをもらった。（DICOMファイル）→k-spaceは入ってなかった。

先に断っておく。

• 医療情報っぽいことを書いたが、自分は医療従事者でもなんでもない。
• 嘘を書いたつもりはない。しかし、間違っていても、責任は持てない。

脳ドッグに行ってきた。→MRIの画像データをもらった。（DICOMファイル）→残念ながら、k-spaceは入ってなかった。
しょうがないので、画像データを二次元フーリエ変換してみた。

MRIとは

MRI（Magnetic Resonance Imaging: 磁気共鳴画像）のこと。
装置は、google検索でトップだった、このサイトの画像のまさにこれ。これの中に入ってきた。

MRIの結果とフーリエ変換

MRIによって、脳とか、心臓の断面画像が手に入る。

この画像の元データがk-spaceとかいうデータで、
これを二次元フーリエ変換すると、断面画像ができるということで、これが欲しかった。
（逆二次元フーリエではないらしい。空間周波数のデータを、フーリエ変換するって何かいつもと違うのがいやだけど。）
脳ドッグでは、MRIを使うそうな。
だから、無保険（16万自腹）でちょっと高めの脳ドッグを受けてみた。

しかし、遊び目的ではない。

たぶんTIA（Transit Ischemic Attack：一過性脳虚血発作）ってやつにかかったかもしれないと思う出来事があった。
だから、脳の精密検査をしたかった。こっちが本当はメイン。

先日、ジムに行っていたときのこと。
週１で、トレッドミル（ランニングマシン）で10km走るのが習慣になっていた。10.5 km/h
換気が悪いなぁと思っていた。

そしたら汗をいつも以上にかいてしまったらしく、脱水症っぽい感じになった。
さらにひどいことに、走り終わった後、手がしびれ、ろれつが回らない状態になった。
こんなん初めて。

ポカリをたくさん飲んだら治った。
先月は、これについて書こうとしていたが、pythonの話とは毛色が違い、話がまとまらなかったため、投稿を断念した。上記のようなことを書こうとしていた。

脳が悪いのか！？と思って怖くなった。

#7119（神戸で救急車が必要か判断できない人のための電話番号）に電話したら、看護師さんにつないでくれて、事情を、上記以上に話して、何個か簡易テスト的なのを指示されて行ったら。「緊急性はないですね。けど、早いうちに医療機関に受診されたほうがいいですよ」とのこと。

もともと、k-space画像が欲しかったので、いわれなくても受けるつもりではあった。
しかし、もともと脳ドッグを予約していた日時（2か月くらい先）を言うと、
看護師さんがちょっと訂正するような口調で「っ、もうちょっと、早いほうが、いいですね」と言われてしまった。

で、直近の土曜日に、脳ドッグを受けた。

その結果が、3週間後（今）に届いた。
結果は、大丈夫だった。
どっか、血管が詰まっているとかはなかった。

ただ、そのほか健康面（コレステロールとか血圧とか）で、所見があったので、結果として行ってよかった。

MRIした感想

• 狭い
• 装置の音がすごかった。
• ヘッドフォン装着→音楽は選べた。

• 体は拘束された状態だった。

• 「息を吸ってー、、、息を止めて。」と何回か指示が来た。

• 「血管を拡げるげるお薬」なるスプレーを舌にかけられた。→寝すぎた時の頭痛と同じ状態に。
  スプレー１吹きで、体にあんな変化があるなんて、「医療すごい」ってなった。

DICOMファイル

医療情報の画像データの規格らしい。
画像のみならず、音声データ等も保存できるそうな。
その他、情報がたくさん詰まっている。

「これに、k-spaceが！？」と思い、
メモ帳で開くと、ところどころ文字化け。バイナリファイルなんでしょう。きっと。

DICOMファイルは、アップロードしない。
なぜならば、

• 自分の個人情報
• 受診した病院名　所在地
• オペレーターの氏名

などが書いてあったから。
自分のは情報はいいとして、自分以外の情報を公開するのは気が引ける。
たとえやっても法律上問題ないとしても気が引ける。

なので、結果画像だけ読み込んでみた。
そしたら、すでに画像化されていた。
k-spaceのデータは含まれていないようです。

まぁそりゃそうだよなぁ。
医療従事者が二次元フーリエ変換なんで遊ぶ必要ないもん。
患者さんを助ける情報が必要なだけだもん。

というか、オペレーターの方が、「k-spaceはあげられないんですよぉー」（にこっ）
っておっしゃっていた気がする。いまさら思い出す。
MRIの技師の方は、二次元フーリエ変換知ってた。やはり使うんだろうか。

dicomファイルの中に、機械の製造者の情報（？）が、Philipsとなっていた。
Philipsといえばシェーバーの会社じゃないですか。
髭剃りだけじゃなかったんだな。
https://ja.wikipedia.org/w/index.php?title=%E3%83%95%E3%82%A3%E3%83%AA%E3%83%83%E3%83%97%E3%82%B9&oldid=69059425

pydicom

pythonでdicomファイルを開くには、このpydicomがいる。

1. dicomファイルを読む(画像のnumpy arrayを含むデータの読み込み)
2. matplotlibで、numpy arrayをimshow
   という、いつもの流れだった。

PythonでDICOM画像をなんとかするここが詳しい。

cmdで、全部小文字で、

pip install pydicom

だった気がする。間違ってたらごめん。
当然ですがpipが入ってないと、効かない。

環境

python2.7
pydicom, numpy, matplotlib, tkFileDialog

コード

• このプログラムは、DICOMファイルを開いて、画像化するプログラムです。
• 二次元フーリエ変換した画像も併せて載せていますが、これは（多分）k-spaceではありません。
• DICOM以外のファイルを開いた時の挙動は意図していません。

import pydicom as dc
import numpy as np
import matplotlib.pyplot as plt


import tkFileDialog as tk
name = tk.askopenfilename(filetypes=[("dicom",("all","*.*"))])
d = dc.read_file(name)


mx = np.max(np.real(np.fft.fft2(d.pixel_array)))/100000.

f,((ax1,ax2),(ax3,ax4))=plt.subplots(2,2)

ax1.imshow(d.pixel_array, cmap="gray")
ax2.imshow(np.abs(np.fft.fft2(d.pixel_array)), vmin=0, vmax=mx*1000, cmap="gray")
ax3.imshow(np.real(np.fft.fft2(d.pixel_array)), vmin=-mx, vmax=mx, cmap="gray")
ax4.imshow(np.imag(np.fft.fft2(d.pixel_array)), vmin=-mx, vmax=mx, cmap="gray")

ax1.set_title("pydicom pixel_data")
ax2.set_title("2DFFT absolute value")
ax3.set_title("2DFFT real part")
ax4.set_title("2DFFT imaginary part")

plt.show()

プログラムの実行方法

上記コードを、テキスト帳にコピーし、「a.py」というファイル名で保存。
↓
そのファイルが表示されているブラウザ上で、普段アドレスが表示されているところに、

cmd

と入力し、Enterキーを押す。
（もしくは、ブラウザ上の何もないところを、shift+右クリックして、「PowerShellWindowをここに開く」をする。）
↓
そうすると、黒い画面（青い画面）が出てくる。
↓
その画面上で、

python a.py

と入力すると、実行できる。ファイルダイアログが出るはず。
（pythonのフォルダが、環境変数で適切に設定してあればの話。）dicomファイルを選ぶと、画像が表示される。

• 画像データはnumpy arrayで返ってくる
• 画像は、mm感覚で撮影されているらしい
  なので、3次元でプロットすれば、立体的に見れるはず。
  大変そうだからまだやらない。