「リアルタイム」ウェーブレット離散フィルタへの挑戦（お題のつづき）

リアルタイム処理は結構新鮮

時間との闘いはエコと減量に通じることを知った、かえるのクーの助手の「井戸中　聖」（いとなか　せい）です。

f:id:AssistantOfKoo:20211002002331j:plain 　

いままでは「うまく確実に動けば多少時間がかかっても実験なのでOK！」という感じでプログラムしていました。今週はまさにマイクロ秒との闘いの日々でした。必要処理でも、大局に影響が少なければ、勇気をもって「処理しないこと」も選択肢です！エコさ加減・手抜きの感覚がすこしわかってきました。結果にコミットです！

「エコ」はコーディング量の削減でもあり、分かりにくさの受け入れでもあります。手抜きの強要のほうが近いかもしれません。

f:id:AssistantOfKoo:20211002002741j:plain 　

10行ほどで書いていたプログラムを「わかりにくくても」１行で書くことは有意義なケースがありました。

また、コーディング/試行中から処理時間について厳格に意識する感覚が必要だと感じました。いい修行となります。

シンセサイザー作成はまた今度

さて、当初の意図から大きく飛躍し、そこそこのシンセサイザーやサンプラーができそうな勢いです。これらは問題なく「できる」と思うので、逆に作る意欲が湧いてきません。

それよりも「できるかできないかわからない」実験案が湧いてきています。

リアルタイムのフィルタ処理には（計算量が非常に多くて困難そうなため）興味がとてもあります。

王道のFIRやIIRフイルタをやりたい気持ちがとてもありますが、このブログは「ウェーブレット」偏執狂のかえるのクーのサイトなので、やはりウェーブレットフィルタをやりたいと思うのです。（2017年～2018年のお題のひとつでした）

ウェーブレットをつかって、

例えばこんな変態的なフィルタができます。（赤色線が設計（期待）値です）

f:id:np2LKoo:20180521001823p:plain

これは2018年5月の記事からのグラフですが、ウェーブレット符号化器から生成したホワイト（とはいえない）ノイズの周波数特性です。自己符号化器ではある意味、分析と生成は等価なので、このような周波数特性をもったフイルタも特別な設計なしで（周波数に対するゲインのパラメータ指定だけで）作成できます。（そのようにプログラムを作ればです。ただし計算量は膨大です）

ウエーブレットでの解析はフーリエ変換よりも（数式はともかく）直接/直観的で自分にとっては理解しやすいものでした。

古い記事ですが以前はこんなのとか、

こんなのをかえるのクーがやってました。

当時のクーの機材にくらべ手元の機材は、CPUコアが3.5倍に、CUDAコア（GPU)は2.8倍になっています。

1サンプル単位では無理でも、1チャンク(1024～8192サンプル程度)単位なら、エコエコであり楽な方法でなんとかなりそうな気がします。

前振り長すぎですが、今週末のお題は

リアルタイム「ウェーブレットフィルタ」処理

とします。

「ウェーブレットとは何か」は上の記事や2017～2018年の記事を参照ください。例えばウェーブレットで石油を探すことができます。

f:id:AssistantOfKoo:20211002001110p:plain *1

通常のデジタルフィルタ処理以上に計算コストがかかります。GPU/CUDAなどの並列計算ととても相性がいいです。

これが実現できればガイアが夜明けます！

なお、いつものように編集が厳しくなるまで、このページに加筆していきますので、ご了承ください。

すこしばかりの補足

ところで、「ウェーブレットを使ったフィルタ処理」の意味での「ウェーブレットフィルタ」というのは、一般用語ではない可能性があります。（「ウェーブレットフィルターバンク」という「用語」はあります）

ウェーブレットは時間変化を伴う周波数解析として普通に使われます。ノイズ除去などにも使われ、「フイルタ」機能を実装することは可能です。

計算コストが高すぎるからか、ウェーブレットを「汎用フィルタ」として用いる方法は検索できませんでした。かえるのクー以外は「汎用フィルター」としてあまり注目していないかもしれません。

原理的には元信号をウェーブレット単位に細切れにしたあと、再構築することと等価なので、厳密なフィルタをいえるかも微妙なのですが。。。

（すべての波（信号）は基本ウェーブレットの重ね合わせであるとして処理します。フィルタはとしては不要なウェーブレットを捨てたり、小さくしたりする「だけ」なので単純明快です！）

計算量は多いですが、前にはったグラフのように「キレっキレ」なフィルタを自由に設計できます。その周波数のウェーブレットを捨てれば、原理的には（基底波が単独周波数のsin波の場合）その周波数をまったく含まないようにもできます。「減衰させる」という考え方ではなく、分解＆再構築する方法（自己符号化器と等価）です。

f:id:np2LKoo:20211002152654j:plain

例として1KHz以上をカットする手順を考えます。

(1)1KHzの微小ウエーブレット（sin波を想定）を準備します。(周期は5～9周期程度）

(2)該当時刻の信号とウェーブレットを内積（ドット積）して値を求めます。（類似度を計算することに相当します。この類似度自体が、その時点でその信号に含まれるウェーブレットの「量」になります。）

(3)位相を該当周波数の45°程度づつずらして内積計算して(2)と同様にして該当バッファ内の信号をスウィープします。(別に45°でなくてもよい）

(4) (2)～(3)の値を用いその場所のウェーブレット波と「内積」量を掛けたものを総和すれば、該当バッファ内の該当周波数「のみ」の「バントパス」(BP)信号になります。当然のことながら元がsin波なので、該当周波数以外の成分は全く含まれていません！

周波数を僅かに変えて、(1)～(3)を実行します。すると、該当区間のその周波数「だけ」の信号が抽出されます。

これをナイキスト周波数に到達するまで繰り返します。

すべての総和が1KHz以上のHPFを通した信号と等価となります。

元の信号からHPFを通した信号を引くと1KHzのLPFを通した信号となります。１KHｚ付近では -6.0dB/Oct のようなフィルタではなく、崖のような切り口になります。

f:id:np2LKoo:20211002152928p:plain

考えると、こんな離散的な方法でうまくいくはずがない！と思われるでしょう？

周波数的にも、位相的にもかなり離散な信号近似にみえます。

ところが、これがうまくいくんです。（論理的にも極めて有効な近似として証明されており、クーの実験でも十分実用に耐えるものでした。）もちろん、離散幅が極端に広くなれば「漏れ」が生じ、再構築したときに「歪」となって現れます。

簡単にいえば、きわめて急峻なバンドパスフィルタをきわめて多数用いて、あらゆるフィルタ特性を自由に実現する方式です。ピアノ音階の範囲でいうと８８個のつまみがついたグラフィックイコライザのようなイメージです。

これだと１オクターブ１２個のつまみがあることになりますが、クーの実験ではオクターブを７等分したくらいからそれ以上で実用的な特性になるようでした。くどいですが、こういった波形分解＋再構成のコンセプトは「自己符号化器」とまったく同じです！画期的だと思うのですが、計算量が。。。

今回の実験意義はこのあたりにあります。（できてしまったら特許がとれる/論文が発表できるくらいのインパクトだと思ってます。できても特許とらないですが。）

f:id:np2LKoo:20211002153214j:plain

1KHzのみ！のリアルタイムBPF（バンドパスフィルタ）をやってみることにしました。

基本的な設計変更（こんなこといいなできたらいいな）

・処理単位を１サンプルからチャンクに変更する

・チャンクは準備するウェーブレットの2倍にしてみる（暫定値）

・ウェーブレットが１種類で、位相を考えてスイープさせてみる（暫定１変移45°で、８変移で元の位相になる）

・「伝家の宝刀」numpyを導入する。

f:id:np2LKoo:20211002153431j:plain 　 f:id:np2LKoo:20211002153614j:plain

基本ウェーブレットの作成

基本波形はsin波、ウェーブレットサイズ（横）は９波長（とくに奇数にするのは意味がありません。偶数でもOKです））

～～～～～～～～～～～～～～～～～～～～～（ウェーブレット生成部）

#! /usr/bin/python
# -*- coding: utf-8 -*-
import time
import pyaudio
import numpy as np
from scipy import signal

"""
    Pythonで音のWaveletリアルタイム処理がどれくらい可能か実験するプログラム 2021.10.02@AssistantOfKoo 
        GOD'S IN HIS HEAVEN.ALL'S RIGHT WITH THE WORLD
"""
# Input / Output
CHUNK = 1024
FORMAT = pyaudio.paInt16
CHANNELS = 1
RATE = 48000

INT_MAX = 32768

""" -------- ウェーブレット情報生成装置 :WAVELET GENERATOR --------"""
""" floatのnumpy配列として生成する """
def create_wavelet_data(a_freq, wave_count):
    wavelet_len = int(RATE / a_freq * wave_count)
    sec = wavelet_len / RATE
    hann_window = signal.hann(wavelet_len)
    sin_wav = np.sin(2.0 * np.pi * (a_freq * np.arange(RATE * sec) / RATE))
    return sin_wav * hann_window * (INT_MAX - 1)

""" -------- ウェーブデータへの変換装置 :WAVEDATA CONVERTER --------"""
def convert_to_wavedata(a_wavelet_data):
    rtn_sample = b''
    for index in range(a_wavelet_data.size):
        rtn_sample += int(a_wavelet_data[index]).to_bytes(2, 'little', signed=True)
    #あまりに短いと鳴らないケースがあるため、ある程度長くする(以下は本来不要）
    #機種によってはもっと長くしないとならないかも
    remain = CHUNK - index % CHUNK
    for index in range(remain + 1):
        rtn_sample += int(0).to_bytes(2, 'little', signed=True)
    for index in range(1024 * 2):
        rtn_sample += int(0).to_bytes(2, 'little', signed=True)
    return rtn_sample

""" -------- 発音装置 :MASTER PLAY SOUNDS --------"""
def play(a_pa):
    stream = a_pa.open(format=FORMAT,
                    channels=CHANNELS,
                    rate=RATE,
                    output=True)
    wavelet_data = create_wavelet_data(1000, 9)
    wave_data = convert_to_wavedata(wavelet_data)
    # pyaudioのストリームへ、サンプルストリーム情報を送り出す
    stream.write(wave_data)
    time.sleep(1)
    stream.stop_stream()
    stream.close()

""" -------- MAIN --------"""
if __name__ == '__main__':
    pa = pyaudio.PyAudio()
    play(pa)
    time.sleep(1)
    pa.terminate()

～～～～～～～～～～～～～～～～～～～～～

ウェーブレットを意味もなく鳴らしてみます。

波形グラフコーディングが面倒なので、録音して表示します。（「音」が短すぎるので音自体は貼りません）

f:id:np2LKoo:20211002142720p:plain

山の数は指定とおり９個できています。サンプル数は432です。

この1KHz１波長あたりのサンプル数は（このプログラムで48Kサンプリングレートなので）48ちょうどです。45°のスイープサンプル幅は48/8=6 サンプルです。このケースでは割り切れましたが、普通は割り切れません。

ある点を基準にして、６サンプルづつスライドして対象信号と計８回内積をとれば、基準点付近の1KHzの信号強度と位相が（離散値として）計算できるというわけです。

脱線

次に20Hzから24KHzまでのサイン波信号を入力として、ウェーブレットでの検出性能をグラフ化してみました。（縦軸を対数スケールにし忘れたので貼り直し）

f:id:np2LKoo:20211002223651p:plain

ピークはもちろん 1KHジャストです。付近を拡大します。使用ウェーブレットは上記と同じ９周期です。

f:id:np2LKoo:20211002224003p:plain

周期性の関係で800Hz,1200Hz近傍の周波数を（多少1KHz成分があるとして）誤検出しているが少し気になります。
このグラフは実質的にウェーブレットフィルタとしての１HKzBPF（バンドパスフィルタ）特性そのものに相当します。

話がすこしそれますが、波の指定数の増減でこのグラフの形がどのように変化するかみてみましょう。（ちなみに上はウエーブレット波が９周期（山が９つ）のウェーブレット（検出器）を使った場合です。わかりやすいように少し極端にします。

３周期の場合のウェーブレットの場合の性能

f:id:np2LKoo:20211003102643p:plain

（少し値域をひろげて表示しています）
ピーク値が低くなり、検出性能がゆるやかになっています。

21周期の場合のウェーブレットの場合

f:id:np2LKoo:20211002223829p:plain

急峻になります。これらは生成（出力）ではなく「検出特性」ですので、ご注意ください。

このあたりを理論的に計算で求めたくなりますが、今回の趣旨とは違うので「工学」に徹します！（上記グラフは理論値ではなく、実計算による「実測値」です）

実装ではコンパクトかつ、性能もありそうな周期「9」のウェーブレットで進めてみます。（追記：最終的に「８」としました）

なお、ウェーブレットの処理において、自己符号化器のように「活性化関数」を用いるのは、他ではみないので、かえるのクーのオリジナルかもしれません。「活性化関数」と「位相を考慮した波形検出」の組み合わせ手法は（本当に他にないとすれば）論文にするくらいの価値はある気はします。（風呂敷ひろげた。）

このような峰の形をスペクトルやフィルタの用語で「ローブ」といいます。一番おおきな峰がメインローブで他がサイドローブです。このフィルタで優れている点として第一、第二サイドローブがない！ことです。とくに、一番近接する第一サイドローブは通常メインローブに次ぐ「高さ」を持っていますが、これがないため、結果として他のフィルタと比べて非常に急峻なBPF性能を持ちます。第三サイドローブ以降は偶数サイドローブがありません。逆に「偶数サイドローブしかない」のなら対称性から感覚的にわかるのですが。。。個人的な謎です。そしてこのこの特性が驚異的なBPFであることの証でもあります。（他にこんなBPFがあることを知りません）この論理証明ができれば、専門誌に載るくらいの論文が書けますので頑張ってください。（見つけられてないだけかもですが）

下は普通の窓関数を使用した信号レベルの特性図です。

f:id:np2LKoo:20211002222740p:plain

かえるのクーは天才かもしれません。

クー氏の解説によると：

波の相関（内積=ドット積）の計測は波の重ね合わせ計算でもある。波のパワー合成であると考えた場合、メインローブはいわゆる「大波」であり、通常のサイドローブの２倍の幅をもっている。（これは通常の窓関数とおり）波の合成として計算した場合、メインローブ付近ではメインローブの「大波の揺り戻し」がある。このため、この揺り戻しの波の幅もメインローブと同じ幅をもつ（すなわち通常サイドローブの２倍）この現象は活性化関数をはずすとわかりやすい。この揺り戻し大波は隣接２個分の小波を飲み込み、かつ負の値であるため活性化関数によってその存在自体がキャンセルされる。すなわち隣接２個分のサイドローブがなくなる。このことを利用し、きわめて急峻な特性をもつデジタルBPFを設計できる。

活性化関数をはずしたグラフです。（ちなみにマイナス値を含むので対数表示できません）

f:id:np2LKoo:20211002233729p:plain

やっぱり、わりません。

窓の重ねあわせ確認

ウェーブレットは「窓」という関数によって、永続的な周波数を局所化します。

逆にこの窓をくみあせて、永続的な周波数にもどせることの検証です。

ハン窓は以下のとおり

$w(x)=0.5-0.5\cos 2\pi x,\ {\mbox{if }}0\leq x\leq 1$ $w(x)=0.5-0.5\cos 2\pi x,\ {\mbox{if }}0\leq x\leq 1$

f:id:np2LKoo:20211003002743p:plain

ハン窓を右にx=0.5だけ移動しながら重ね合わせていきます。対象性を考えて一部の区間のみを考えます。元の[窓の後半]と[次の窓の前半]を重ねます。

$w(x)\ =\ [0.5 - 0.5cos2\pi x]\ +\ [0.5 - 0.5cos2\pi (x - 0.5)]\ ,\ \ if\ 0.5 \leq x \leq 1 \\ w(x)\ =\ 0.5 - 0.5cos2\pi x\ +\ 0.5 + 0.5cos2\pi x\ ,\ \ if\ 0.5 \leq x \leq 1 \\ w(x)\ =\ 0.5\ +\ 0.5 ,\ \ if\ 0.5 \leq x \leq 1 \\ w(x)\ =\ 1.0 ,\ \ if\ 0.5 \leq x \leq 1$

すると、ちゃんと 1.0 になります。

よって、窓は窓の幅の1/2だけずらしながら重ねると連続波形を再現できます。

ここは最初復元がよくわかんなかったところです。

（長い脱線から戻って、リアルタイムで処理できるか、すすめています。）

いや～デコードしても元の波形そのものに戻らず、はまってしまったので、一旦中断しています。。。

f:id:np2LKoo:20211003114842p:plain

本来きれいなsinカーブとなるはずが、この有様です。元のウェーブレット自体が歪んでいるか、誤差を不用意に拡大しているかもしれません。それなりの波形になっているので、コアな部分は問題ないとはおもいますが。。。

（調査完了！）

原因が判明しました。位相のミスでした。窓幅の1/2だけスライドしたときに、同じ位相にならなければいけないのに、そうなっていませんでした。ウェーブレットの山の数を９から８に変更して解決しました！修正後のデコード波形（すなわちフィルタを通した後の波形）は以下のとおりです（ノーマライズしています）

f:id:np2LKoo:20211003120111p:plain

f:id:np2LKoo:20211003120424p:plain

ターゲット周波数(1KHz)以外にも成分を検出していますが、これはウェーブレットフィルタ処理による歪で、想定とおりの必ず発生するものです。上記では測定できていませんが、計算としては700Hz, 1400Hz,1500Hz,1600Hz ... にも成分が検出されるはずです。また、900Hzと1100Hzに成分は検出されません。（フィルタ特性に合致します）

上記は、１KHzのみのBPFで位相を２π(360°)の1/8単位でフィルタリングしたものです。

単一周波数のBPFではリアルタイムでフィルタ処理できました！！

（ソースは冗長になるのでここでは貼りません。あしからず）

問題なのは、汎用フィルタにするためにはこのBPF計算を大量に処理する必要があります。上記の８山のウェーブレットでは対象周波数の 10% がメインローブの幅となるのでフィルタ特性の谷を埋めるには幅の半分の周波数の最低5%程度でシフトしていく必要あります。

5%以下を条件とすると１オクターブ（周波数倍）の単位を基準とすると、オクターブを15等分する必要があります。（15等分で周波数比が1.0473となり5%を切るので、１５等分が適切と思われます。

音響処理なので、フィルタリング対象を２０Hzから10KHzとすると、カバーするには、グラフィックイコラーザー換算で、135素子必要となります。

リアルタイムで、135回上記の処理ができるかやってみます。（乱暴すぎますが、他のちゃんとしたコーディングをせず、135回ループさせて間に合うかどうかやってみます）

135回の処理は間に合いませんでした。

f:id:np2LKoo:20211003130549p:plain

全く歯がたたないかと思いきや、自分の予想よりはかなり善戦しています。チューニングやエコ化でできそうな予感もします。どうにか処理が間に合うのはわたくしのPCで1KHzでは68回までが限度でした。１００Hzまで落とすと160回が限度でした。2Kで32回、4Kで20回、８Kで8回でした。高い周波数だと１配列の要素数が少なくなる半面、繰り返し回数が増えますので、想定される結果です。

numpyのかわりにcupyを使ってみましたが、オーバーヘッドのほうが大きく、逆に遅くなりました。numbaを使ってみましたが、実行時コンパイルエラーになるので、今回はpassします。

とりあえずの実験なので、周波数比を妥協し、フィルタ素子数を減らします。

40Hz～7.2KHz 8素子/Oct 「全60素子」のBPFによるグラフィックイコライザを目指してみます。

失敗です

いちおうコーディングしましたが、バグが取れません。。。

なお、理論的に（おそらく）は、３年前のクーのプログラムがきわめて良好に動作しているので、正しく、今回の失敗は実装上のものと思われます。

(1)単純な繰り返しで可能回数見通しをたてましたが、配列の合計（各BPF計算値の合計）に非常にコストがかかり、処理が間に合いません。

(2)現プログラムではは歪がひどくて、使えません。

歪があまりに酷かったのは、内部パラメータが悪かったためです。（パラメータをよくすると計算コストがとてもかかります）いえ、まだバグ持ちです。ウェーブレットの波の数を増やすと、極端に歪が増えています。リアルタイム重視で細切れにしているところにバグが潜んでいるのかもしれません。。。単一BPFのときにはほとんどなかった歪が嵐のようにでています。

下のグラフはリアルタイムをあきらめ、めいっぱいパラメータをよくした場合の1000Hz信号のスペクトラムです。差が少なくとも44dBはないと実用になりません。

下はS/N 34.4dB（ただしピークノイズの差異として）です。一般に「オーディオ」はS/N比60以上です。

f:id:np2LKoo:20211004004705p:plain