Python:ビブラートとコーラス:2

前回の記事ではPythonでビブラートエフェクトとコーラスエフェクトを実装した。
今回はFM変調とBessel関数に関して面白い記述*1があったのでほとんどパクりだけど(情報は多いほうがいいに決まってるw)それについてみていく

FM変調

FM変調は以下の式で示されるのだった。

  • y(n) = Asin(2 \pi f_{c} \frac{n}{f_{s}} + Isin(2\pi f_{m}\frac{n}{f_{s}}))
  • A:原音信号振幅
  • fc:原音信号周波数
  • fs:サンプリング
  • fm:変調周波数
  • I:変調深度(エフェクターによくあるdepthつまみ)

このままでもいいのだが書くのがめんどいので2 \pi fの部分を\omegaに置き換えさせてもらう
置き換えるとこうなると思う。

\begin{align*}
y(t) = a_c \sin(\omega_c t + a_m \sin(\omega_m t))
\end{align*}

ここで
\begin{align*}
t &= n / fs\\
\omega_c&= carrier freq\\
\omega_m&=modulation freq
\end{align*}

キャリアは原音信号のこと

この式を加法定理を使って展開しておく。

\begin{align*}
\sin (\alpha+\beta)=\sin\alpha\cos\beta+\cos\alpha\sin\beta
\end{align*}


\begin{align*}
y(t) = a_c (\sin(\omega_c t)\cos(a_m \sin(\omega_m t)) + \cos(\omega_c t)\sin(a_m \sin(\omega_m t)) )
\end{align*}

第一種ベッセル関数

第一種ベッセル関数はベッセル微分方程式の解のひとつである。
第二種とか、特殊とかあるので詳しくは参考文献を見てくれたほうがいいと思う。
形としてはこんなグラフです。

f:id:gsmcustomeffects:20180819070610p:plain

ベッセル関数とFM基本式

参考文献よりベッセル関数から以下の関係式が得られます。

\begin{align*}
\cos(x \sin(y)) &= \sum_{n=-\infty}^{n=\infty} J_n (x) \cos (ny) \\
\sin(x \sin(y)) &= \sum_{n=-\infty}^{n=\infty} J_n (x) \sin (ny)
\end{align*}


それを用いて上の加法定理の式をそれで置き換えると
\begin{align*}
y(t)=&
a_c
\sum_{n=-\infty}^{n=\infty} J_n (a_m) (\sin(\omega_c t)\cos(n \omega_m t) + \cos(\omega_c t)\sin(n \omega_m t)) \\
=&
a_c
\sum_{n=-\infty}^{n=\infty} J_n (a_m) \sin ( (\omega_c + n \omega_m)t )
\end{align*}

となる。

もう一度まとめてみると

\begin{align*}
a_c \sin(\omega_c t + a_m \sin(\omega_m t)) = a_c
\sum_{n=-\infty}^{n=\infty} J_n (a_m) \sin ( (\omega_c + n \omega_m)t )
\end{align*}

ここからFM変調はベッセル関数と基音+倍音の積で表すことが可能だとわかる。

変調指数の大きさと倍音の関係

次に変調指数についてみていく


第一種ベッセル関数の\(J_{n}(\alpha)\)は一般的に変調指数と呼ばれる。
上で導いたようにFMの変調深度もそれにあたる。

というわけでこいつを大きくしていけばエフェクトとしては深くかかっていく理論。
その時ベッセル関数での挙動はどうなっているのか出してみた。
x軸を第一ベッセル関数の次数にしてy軸をJn(固定)にして表示すればよい

Jn(0.5)の場合

f:id:gsmcustomeffects:20180819071931p:plain

Jn(5)

f:id:gsmcustomeffects:20180819072006p:plain

Jn(10)

f:id:gsmcustomeffects:20180819072016p:plain

変調指数が低いと原音がきちんと聞こえるエフェクトとなるがJn(0.5)の時を見ると元スペクトルが大きく、側波帯が低くなっています。
第一種ベッセル関数は徐々に振幅が下がっていく特性を持つので側波帯のスペクトルは速やかに減衰する。すなわちエフェクトの効果通り原音の色が強い音色となる。


この記事で使ったPythonコード

scipy.specialモジュールを使います

from scipy.special import jv, yv
import numpy as np
import matplotlib.pyplot as plt

z = np.arange(-20,20,0.01)#x軸の範囲

for i in  range(0,5):#次数の決定(別に何個でもいい)
    bessel = (jv(i,z))
    plt.plot(z,bessel,label="J"+str(i))

plt.legend()
plt.show()

nを横軸にしたグラフのコード

from scipy.special import jv, yv
import numpy as np
import matplotlib.pyplot as plt
z = np.arange(0,50,0.01)
x = []
bessel = []
jn = 0.5
jn_index = int(jn*100)
for i in  range(0,30):
    bessel.append(jv(i,z[jn_index]))
    x.append(i)

plt.xlabel("bessel order")
plt.ylabel("Jn(" +str(jn) + ")")
plt.scatter(x, bessel)
plt.show()

参考文献

*1:Analyzing an FM signal [John d .cook氏]

Python:ビブラートとコーラス:1

今回はFM変調を用いたビブラートとコーラスを実装していきたいと思います。

原理

遅延回路(アナログで言うところのBBD)のクロックをLFO(変調波形)で変えることにより周波数方向に揺れた音を作るエフェクトです。

基本式は以下のようになります。

  •  y(n) = Asin(2 \pi f_{c} \frac{n}{f_{s}} + Isin(2\pi f_{m}\frac{n}{f_{s}}))
  • A:原音信号振幅
  • fc:原音信号周波数
  • fs:サンプリング
  • fm:変調周波数
  • I:変調深度(エフェクターによくあるdepthつまみ)

sinの中にsinがあって.....うーんって感じがしますが、そのまま書いてグラフ出すと原理は理解しやすい
f:id:gsmcustomeffects:20180815160538p:plain

画像は上から原信号、変調波系、被変調波となっています。
時間によって周波数が変化している。
一方でトレモロの時はAM変調で単なる乗算だったので
\begin{align}
y = sin a\cdot sin b
\end{align}

な基本式が信号処理っぽく詳しく書くと
\begin{align}
y(n) = Asin(2 \pi f_{m}\frac{n}{f_{s}}) \cdot sin(2 \pi f_{c}\frac{n}{f_{s}})
\end{align}

こうなってx[n]を入力サンプルとすると

\begin{align}
x[n]\cdot Asin(2 \pi f_{m}\frac{n}{f_{s}})
\end{align}

に書き直せば一般的な問題として考えることが可能だった。

というわけで?

こいつらをどうやって実装していくかって話。
デジタル信号処理では遅延処理はメモリ(バッファ)を使う事で実現できる。
ということは入力x[n]をメモリM[n]に入れてnを変調してやればどうにかなりそう?

\begin{align}
y[n] = x[M[n]]
\end{align}

このように入れ子の構造にしてnを変調する

図で書いてみるとこんな感じになった。
f:id:gsmcustomeffects:20180815215038p:plain

次に変調ありとなしの場合のデータの流れを見ていく

f:id:gsmcustomeffects:20180815220144p:plain

変調なしの場合は順番にデータが取り出されるのでインデックス番号は一定になるが変調ありの場合はインデックスを変調したためインデックス番号が小数のところが現れる(というか整数になるとこが稀)。
まあC言語をやっている人はわかると思うが配列のインデックス番号は整数でないといけない

float buffer[10];

a = buffer[1]    // OK
b = buffer[1.1]//NG

というわけでキャストとかを使うと変調後のnは0,01,1,2,3.....のようになり複数回同じデータをyに出力することになる。
入力周波数が低い場合はいいが周波数が高いとサンプル点が少なくなるのでなめらかな波形でなくなる。
別にこの状態でもコーラスエフェクトとして使えるが一般にはこの現象を防ぐために補間を使って間のデータを補う手法が用いられる。*1

というわけで実装していきましょうか
まずは入力サンプルを作っていくので正弦波合成からですね。いきなりwaveファイルを通して実験してもいいのですがグラフだした時に確認しにくいのでこの方式にします。
正弦波の作り方はここに書いてますのでそれを参考にしてください

input = cf.create_sin(1.0,440,frame,sampling=48000)

まあこんな感じで良いでしょう
次はインデックスの変調を行う部分を作ります。

def index_mod(frame):
    n = np.zeros(frame, dtype=np.float)
    i_num = 0
    for i in range(0,frame):
        n[i] = i + depth_ms * (np.sin(i_num*(2.0 *np.pi * freq/fs)))
        i_num += 1
    return n

ここは別にどう書いてもいいですけど変調後のインデックスを返す関数を一個作っておきました。

あとはこういうふうにつかえばOK

n = index_mod(frame)#インデックスの変調
output = np.zeros(frame,dtype=np.float)#空の配列作成
for i in range(0,frame):
    output[i] = mix*input[int(n[i])] + (1-mix)*input[i] #原音と変調音を足している

サンプルサウンド

原音

エフェクト音(コーラス)

エフェクト音(えぐいビブラート)

実際のwavファイルのをつかった例

結構試行錯誤的な部分があるのでくそコードなのは勘弁していただきたい。
補間については別の記事でやりたいので今回は無補間で処理してます
無処理な割にはそこそこきれいな音になっているはず

import numpy as np
import matplotlib as pl
pl.use('TKagg')
import matplotlib.pyplot as plt
import wave
import audio_func as af
import struct
import scipy

def index_mod(frame):
    n = np.zeros(frame, dtype=np.float)
    i_num = 0
    for i in range(0,frame):
        n[i] = i + depth_ms * (np.sin(i_num*(2.0 *np.pi * freq/fs)))
        i_num += 1
    return n

wf = wave.open("GS05.wav", "r")#wav 扱うならお決まりのやつ
num_data = scipy.fromstring(wf.readframes(wf.getnframes()),dtype = "int16") / 32768.0#正規化
fs = wf.getframerate()
frame = wf.getnframes()

freq        = 9     #変調周波数
depth       = 6     #変調深度
depth_ms    = int(fs * depth /1000)
mix = 0.99          #ミックス(コーラスで言うエフェクト量高いほどヴィブラートになる)

if(wf.getnchannels() == 2):
    left = num_data[::2]#1 スライス
    right= num_data[1::2]
    in_data = left*0.5 + right*0.5#2ステレオモノラル化
    #1:スライスの説明
    #a[1,2,3,4,5]ていうリストがあったとして
    #a[::2]  -> 1,3,5
    #a[1::2] -> 2,4

    #2:ステレオ->モノラル化
    #もともとギターの録音で左右に同じ音ふってるのでLR分解して半分にして足せば同じようなもん


input = in_data

n = index_mod(frame)#インデックスの変調
output = np.zeros(frame,dtype=np.float)#空の配列作成


for i in range(0,frame-9000):#変調をあげすぎるとインデックスオーバーするから適宜調整する
    output[i] = mix*input[int(n[i])] + (1-mix)*input[i]
plt.subplot(211)
plt.plot(input[0:frame])
plt.subplot(212)
plt.plot(output[0:frame])
plt.show()

#正規化したデータを元に戻す
output = [int(x * 32768.0) for x in output]
output = struct.pack("h" * len(output), *output)
#af.play(output,wf.getsampwidth(),1,wf.getframerate()) #再生
af.save(output, fs,1,"custom_vib3.wav")

*1:音遊び!Blackfin DSP基板でディジタル信号処理初体験を参照

Python:正弦波合成

デジタル信号処理での正弦波の作り方

create_func.py

import numpy as np

def create_sin(amplitude,freq,frame,*,sampling = 48000):
    data = np.zeros(frame,dtype = np.float)
    for n in np.arange(1,frame):
       sine = amplitude * np.sin(2 * np.pi * freq * n / sampling)
       data[n] = sine
    return data

def create_cos(amplitude,freq,frame,*,sampling = 48000):
    data = np.zeros(frame,dtype = np.float)
    for n in np.arange(1,frame):
       cosine = amplitude * np.cos(2 * np.pi * freq * n / sampling)
       data[n] = sine
    return data

使い方

import audio_func as af
sinewave = cf.create_sin(1.0,1000,frame,sampling=48000)




参考文献

Python:リングモジュレーターの実装②

前回はリングモジュレータを実装したのでもう少し細かく見てみようということで数式をあげながら理論を追ってみる

前回の記事は此方
gsmcustomeffects.hatenablog.com


原理

被変調波(原音)を sin(2 \pi f_{c}\frac{n}{f_{s}})と置くとリングモジュレーションは以下の式で示すことができる。

  •  y(n) = Asin(2 \pi f_{m}\frac{n}{f_{s}}) \cdot sin(2 \pi f_{c}\frac{n}{f_{s}})

これを積和の公式(高校の加法定理とこで習う)を使って書き直すと

  •  sin\alpha sin\beta = -\frac{1}{2}\{cos(\alpha + \beta) - cos(\alpha - \beta)    \}
  •  y(n) = \frac{1}{2}[cos(2\pi(f_{c}-f_{m})\frac{n}{f_{s}})-cos(2\pi(f_{c}+f_{m})\frac{n}{f_{s}})   ]

fcを1kHz,fmを2kHzとした時のことを考えてみる。

  •  y(n) = \frac{1}{2}[cos(2\pi(1000-2000)\frac{n}{f_{s}})-cos(2\pi(1000+2000)\frac{n}{f_{s}})   ]
  •  y(n) = \frac{1}{2}[cos(2\pi(-1000)\frac{n}{f_{s}})-cos(2\pi(3000)\frac{n}{f_{s}})   ]

cosは偶関数なのでf(x) = f(-x)

  •  y(n) = \frac{1}{2}[cos(2\pi(1000)\frac{n}{f_{s}})-cos(2\pi(3000)\frac{n}{f_{s}})   ]

というわけで0.5倍された1kHzのcos と3kHzのcosが出てきそうである。

Pythonでグラフだしてみると同じになった

f:id:gsmcustomeffects:20180811004424p:plain

実験室!

理論がだいたいわかったとこで変調かけた音がどういうふうに聞こえるか検証していく
今回もゆきょん(@Yukyoooon)君に音源とってもらいました。
twitter.com

変調とのことでアルペジオ音源です.

ちなみにグラフは

  • 上段 : 原音
  • 中段 : 変調波形
  • 下段 : 加工音

原音

まずは原音から

100Hzのsin波形

高い周波数をかけすぎるとほんとに音程がわからなくなるので100Hz

[f:id:gsmcustomeffects:20180812204119p:plain]

10Hzのsin波形

[f:id:gsmcustomeffects:20180812204007p:plain]

5Hzのsin波形

[f:id:gsmcustomeffects:20180812212939p:plain]

1Hzのsin波形

[f:id:gsmcustomeffects:20180812212533p:plain]

これに関しては変調波形が遅いこともあってトレモロっぽい

5+10Hzの波形

  • wave1 : 振幅0.5 5Hz
  • wave2 : 振幅0.5 10Hz

[f:id:gsmcustomeffects:20180812214705p:plain]

1+10+100Hzの波形

  • wave1 : 振幅0.3 1Hz
  • wave2 : 振幅0.3 10Hz
  • wave3 : 振幅0.3 100Hz

[f:id:gsmcustomeffects:20180812215151p:plain]

1+10+100Hzの波形その2

さっきのだと全部の成分の主張が激しすぎるのでバランス調整した時の音

  • wave1 : 振幅0.7 1Hz
  • wave2 : 振幅0.1 10Hz
  • wave3 : 振幅0.1 100Hz

[f:id:gsmcustomeffects:20180812215544p:plain]

sin a * sin b * sin c + sin dの波形

  • wave1 : 振幅 1.0 1Hz
  • wave2 : 振幅 1.0 2Hz
  • wave3 : 振幅 1.0 3Hz
  • wave4 : 振幅 0.2 500Hz

wave1 * wave2 * wave3 +wave4って感じで積和の公式から変な周波数生成してそれに対して振幅弱めの高周波を加算するというアイデア

\begin{equation}
sin x sin 2x sin 3x = \frac{1}{4}(sin2x +sin4x - sin6x)
\end{equation}

前の三項はこういう感じに機能する。

波形の原型は3項がつかさどり味つけに4項目が効くという感じ

f:id:gsmcustomeffects:20180812223510p:plain

んで最後に音源


めんどくさいのでやらなかったアイデア

  • 加法定理から高い音が加えられるということはなんとなくわかったのでフィルターを組み合わせると面白いのかなと思った
  • 波形同士の引き算で倍音を消せるかも?
  • 変調波形をLFOでFMするとか?

あとこの方のdome filterも面白いのかなと

hhh-yama.hatenadiary.jp

まとめ

リングモジュレータの続きという形でもう少しまともにやってみた。
加法定理忘れてて焦った。

この辺までやっていろいろ面白いことが分かったのであとは実機でつまみいじりながら調整したいと思いました。
次回はトレモロでもやります。

Python:ScipyでFFTを行う

scipyのfftメソッドは二つありますがscipy.fftを使うことが推奨されてるようです。

コード

import matplotlib.pyplot as plt
import numpy as np
from scipy.fft import fft, fftfreq

f0 = 1000
fs = 96000
N = 5000
addnum = 5.0

def create_sin_wave(amplitude,f0,fs,sample):
    wave_table = []
    for n in np.arange(sample):
        sine = amplitude * np.sin(2.0 * np.pi * f0 * n / fs)
        wave_table.append(sine)
    return wave_table

# 信号波形の生成
wav = create_sin_wave(1.0,f0,fs,N)

# FFTの実行
X = fft(wav[0:N])
freqList = fftfreq(N, d=1.0/ fs)
amplitude = np.abs(X)/(N/2)

# 波形を描画
plt.subplot(211) 
plt.subplots_adjust(hspace=0.5)
plt.plot(range(0,N), wav[0:N],label = "wave1")
plt.axis([0, 500, -1.0, 1.0])
plt.xlabel("sample")
plt.ylabel("amplitude")

# 振幅スペクトルを描画
plt.subplot(212)
plt.plot(freqList, amplitude, marker='.', linestyle='-',label = "fft plot")
plt.axis([0, 20000, 0, 2])
plt.xlabel("frequency [Hz]")
plt.ylabel("amplitude")

plt.show()

FFT結果

f:id:gsmcustomeffects:20220305222226p:plain

Python:リングモジュレーターの実装①



今回はリングモジュレータを実装していく。

リングモジュレータ

入力と関係ない音程を出すエフェクター。マルチ・エフェクターに組み込まれていることもある。かなり前衛的な演奏ができる。(weblio引用

音としてはこんな感じ

元の音

エフェクト音

理論としては入力の音に対してある周期波形を乗算してめちゃくちゃな音程にするエフェクトです。
変調波形を正弦波として表すとこのような感じ

  • out[ n] = in[ n] \times Asin2\pi f_{0} t
  • t = n / fs

A:正弦波の振幅
f0:正弦波の周波数[Hz]
n:その時のサンプル
fs:サンプリング周波数[Hz]

エフェクトの波形

f:id:gsmcustomeffects:20180801023157p:plain

  • 上段:入力音
  • 中段:変調波形(200Hzの正弦波)
  • 下段:出力音

かけている正弦波の影響を受けて波形が変形している。

ソースコード

処理の流れはこのようになる

  1. 元データがステレオなのでwavファイルをLRにわける
  2. モノラル化する(LRを半分にして加算)
  3. 正弦波を作る
  4. エフェクトをかける
  5. wavに格納しなおす
import wave
import audio_func as af
import scipy
import struct
import numpy as np
from pylab import *

wf = wave.open("GS04.wav", "r")#wav 扱うならお決まりのやつ
num_data = scipy.fromstring(wf.readframes(wf.getnframes()),dtype = "int16") / 32768.0#正規化

if(wf.getnchannels() == 2):
    left = num_data[::2]#1 スライス
    right= num_data[1::2]
    in_data = left*0.5 + right*0.5#2ステレオモノラル化
    #1:スライスの説明
    #a[1,2,3,4,5]ていうリストがあったとして
    #a[::2]  -> 1,3,5
    #a[1::2] -> 2,4

    #2:ステレオ->モノラル化
    #もともとギターの録音で左右に同じ音ふってるのでLR分解して半分にして足せば同じようなもん

fs = wf.getframerate()
AM_WAVE = []#振幅変調波形のメモリ確保
f0 = 200      #正弦波の周波数[Hz]

#フレームの長さ分だけ正弦波を作る
for n in np.arange(wf.getnframes()):
    sine = 0.2 * np.sin(2 * np.pi * f0 * n / fs)
    AM_WAVE.append(sine*(1/0.2))

#リングモジュレータなので振幅変調する
out_data = in_data * AM_WAVE
subplots_adjust(hspace=0.5)
subplot(311)
plot(in_data[0:6000],label="input")
legend()
subplot(312)
plot(AM_WAVE[0:6000],"red")
subplot(313)
plot(out_data[0:6000],label="output")
legend()
show()
#正規化したデータを元に戻す
out_data = [int(x * 32768.0) for x in out_data]
out_data = struct.pack("h" * len(out_data), *out_data)
#af.play(out_data,wf.getsampwidth(),1,wf.getframerate()) #再生
af.save(out_data, fs,1,"ring.wav")

audio_funcは一例です。(かなり雑な実装なため

import wave
import pyaudio
from pylab import *

def printWaveInfo(wf):
    """WAVEファイルの情報を取得"""
    print("チャンネル数 : "+ str(wf.getnchannels()))
    print("サンプル幅 : "+ str(wf.getsampwidth()))
    print("サンプルレート : "+ str(wf.getframerate()))
    print("フレーム数 : "+ str(wf.getnframes()))
    print("総パラメータ(一括表示用) : "+ str(wf.getparams()))
    print("再生時間 : "+ str(float(wf.getnframes()) / wf.getframerate()))

def play (data,sampleWidth,Channel,sampleRate):
    # ストリームを開く
    p = pyaudio.PyAudio()
    stream = p.open(format=p.get_format_from_width(sampleWidth),
                    channels=Channel,
                    rate=int(sampleRate),
                    output= True)
    # チャンク単位でストリームに出力し音声を再生
    chunk = 1024
    sp = 0  # 再生位置ポインタ
    buffer = data[sp:sp+chunk]
    while buffer != '':
        stream.write(buffer)
        sp = sp + chunk
        buffer = data[sp:sp+chunk]
    stream.close()
    p.terminate()


def save(data, fs,Channel,filename):
    """波形データをWAVEファイルへ出力"""
    wf = wave.open(filename, "w")
    wf.setnchannels(Channel)
    wf.setsampwidth(2)
    wf.setframerate(fs)
    wf.writeframes(data)
    wf.close()