ガレスタさんのDIY日記

電子回路、Web、組み込み、自作エフェクターを語るblog

Python:ディレイを実装してみる:1

今回はディレイを実装していきます。
いわゆる遅延系エフェクトでやまびこ効果とか言われるやつです。


ディレイの実装には参考文献に示すようにフィードバックの仕方によっていろいろな種類がありますが一般的なフィードバックありミックスありの通常のタイプをやって行きます。

f:id:gsmcustomeffects:20180820231509p:plain

まず時系列信号処理をするにあたってデータを保存するバッファを確保する必要があるのでリングバッファを使ってそれを実現します。
リングバッファの説明は参考文献の二つ目を参照してくれると理解しやすいと思う。

実装

というわけでさっそく実装していく

まずはエフェクトのパラメータであるfeedbackとmixとtimeを定義する。

delay_buffer_idx    = 0                                 #リングバッファのインデックス
feedback            = 0.4                               #フィードバック(0~1)
delay_ms            = 465                               #ディレイタイム[ms]
delay_buf_size      = int(frame_fs * delay_ms / 1000)   #バッファ確保[sample]
level                 = 0.1                               #level(0~1)

frame_fs は処理系によって異なりますが基本48000(48kHz)ぐらいだと思います。
48kHzサンプリングなので48000サンプルで1秒のディレイタイムが確保できる。
今回はPC上での実装なのでメモリはほぼ無限にあるがDSPだったり組み込みでの実装の場合限られた中で実装する必要もあるので頭に入れておくといいです。

次に入出力関係の整備を行う
ここはいつも通り入力バッファと出力バッファを確保しそれと今回使うディレイバッファを確保する。

input_buf    = in_data#wavデータだったり正弦波だったりを格納しておく
output_buf   = np.zeros(frame_size,dtype=np.float)#
delay_buffer = np.zeros(delay_buf_size, dtype=np.float)#パラメータのところで計算したサイズ分だけ確保する。

信号処理部分はコメントに示す通りの流れになる。

"""signal  processing"""
for i in range(frame_size):
    l_x = input_buf[i]*1.0+ delay_buffer[delay_buffer_idx] * level #入力(ゲインは1)+ディレイ
    delay_buffer[delay_buffer_idx] = input_buf[i] + delay_buffer[delay_buffer_idx]*feedback#入力とフェードバック分をバッファに書き込み
    delay_buffer_idx = int((delay_buffer_idx + 1) % delay_buf_size)#リングバッファをすすめる
    output_buf[i] = l_x#出力

これで処理自体は完了となる。
処理としてはメモリに入れて所定サンプル数経過後に原音と加算して出力してるだけなのでかなりコードが短い

最後に処理前後の波形を示しておきます。
f:id:gsmcustomeffects:20180821001507p:plain

音が切れてから残響効果により音が残る波形となっていることがわかる

サンプル音源

元音源

ディレイ

各種設定

"""delay param"""
delay_buffer_idx    = 0                                 #リングバッファのインデックス
feedback            = 0.4                               #フィードバック
delay_ms            = 465                               #ディレイタイム[ms]
delay_buf_size      = int(frame_fs * delay_ms / 1000)   #バッファ確保[sample]
level                 = 0.8                               #ミックス(エフェクトレベル)

コード全体

かなり走り書きなので参考程度にお願いします。

import numpy as np
import matplotlib.pyplot as plt
import create_func as cf
import wave
import audio_func as af
import struct
import scipy

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

"""delay param"""
delay_buffer_idx    = 0                                 #リングバッファのインデックス
feedback            = 0.4                               #フィードバック
delay_ms            = 400                              #ディレイタイム[ms]
delay_buf_size      = int(frame_fs * delay_ms / 1000)   #バッファ確保[sample]
level                 = 0.5                               #ミックス(エフェクトレベル)
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_buf    = in_data
output_buf   = np.zeros(frame_size,dtype=np.float)
delay_buffer = np.zeros(delay_buf_size, dtype=np.float)


"""signal  processing"""
for i in range(frame_size):
    l_x = input_buf[i] + delay_buffer[delay_buffer_idx] * level #入力+ディレイ
    delay_buffer[delay_buffer_idx] = input_buf[i] + delay_buffer[delay_buffer_idx]*feedback#入力とフェードバック分をバッファに書き込み
    delay_buffer_idx = int((delay_buffer_idx + 1) % delay_buf_size)#リングバッファをすすめる
    output_buf[i] = l_x#出力

plt.subplot(211)
plt.plot(input_buf[0:frame_size])
plt.subplot(212)
plt.plot(output_buf[0:frame_size])
plt.show()

output = [int(x * 32768.0) for x in output_buf]
output = struct.pack("h" * len(output), *output)
af.play(output,wf.getsampwidth(),1,wf.getframerate()) #再生
#af.save(output, frame_fs,1,"delay2.wav")

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基板でディジタル信号処理初体験を参照

音声信号処理まとめ

ここ一年でオーディオ処理について学んできたのでいろいろ書いていこうと思います。
最近は個人で扱えるマイコンも高性能化してきたり音声ライブラリも豊富にあり逆にどれからやったらいいかわからない人も多いはずです。

そういうわけで自分がやっている例でも紹介できたらと思い記事のまとめを作りました。

2018/8/13
最近は東京と新潟行ったり来たりなのでハードウエアネタは少な目でPythonでエフェクトアルゴリズムについて個人で研究した内容を書いてます。
そのうち書籍まとめも出来たらいいと思ってます。


マイコン上での実装

エフェクターは実機で動かしてこそなのでココでは高性能ARMマイコン上にオーディオフレームワークを実装し処理する例をまとめてみた。
主にSTM32F746(単精度FPU),STM32F767(倍精度FPU),NXP i.MX RT1020関連ネタが多め

SAIペリフェラルでトークスルーを動かす - ガレスタさんのDIY日記(マイコンで信号処理するときの基本)
CMSIS DSP:FIRフィルタをやってみる - ガレスタさんのDIY日記
CMSIS DSP:biquad-low-passをやってみる - ガレスタさんのDIY日記
STM32F7:トレモロ(リングモジュレータ)の実装 - ガレスタさんのDIY日記

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

まとめ

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

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