20181015

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()

ランダムノイズ(全体).png
ランダムノイズ(ズーム).png

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()
ホワイトノイズのスペクトルも、ホワイトノイズ。
全周波数にわたって、値は同程度(の振幅のノイズ)。
ホワイトノイズのスペクトルも、ホワイトノイズ.png
これがホワイトノイズの良くないところ。
人の耳は、周波数ごとに感度が違う。高周波には敏感。
どの周波数も同程度に感じるような音がピンクノイズらしい。
ピンクノイズの前に、ブラウニアンノイズを作る。

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

ブラウニアンノイズとかの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()

1/f^2.png

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()

1/f^2の修正.png

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

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()

1/f^2を、ホワイトノイズのスペクトルにかけた.png
このスペクトルは、このままだと、まずい。
まずい理由は、

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

試しに、ifftかけると

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

低周波が優勢すぎる.png
のっぺりしてしまっている。

figure_2-2.png

スピーカーの再生できる周波数に絞る(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()

可聴域のみに加工.png

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

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()

可聴域のみに加工 高さを合わせた.png

逆変換すると、

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()

逆フーリエ変換した(全体).png

逆フーリエ変換した(ズーム).png

 

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
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()

1/f.png


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

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()

1/fを、ホワイトノイズのスペクトルにかけ可聴域に加工.png

ピンクノイズは、ブラウニアンよりも減りがなだらか。
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()

逆フーリエ変換した(ピンクノイズ).png

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

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

この後、wav作ったけど…

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

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

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

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

posted by yuchan at 07:00 | Comment(0) | python