※この記事はMMA Advent Calendar 2024 - Adventarの20日目の記事です。
はじめに
こんにちは。Udonです。
このアドベントカレンダーも20日目となりました。ご協力いただきありがとうございます。
前日の記事はfjdjさんのドライブについてダラダラ語る②でした。文字数がすごいです。勝手な印象ですが隙あらば車に乗っているような気がするので、そういった経験がしっかり文章化できていてすごいと思いました。またドライブ誘ってください(奥多摩とか行きたい)。
というわけで、自分の番になったわけなので、今回は「音声処理のプログラム」について書いていこうかなと思います。
これは先日大学の実験で扱った題材でした。具体的には、雑音が重畳した音声ファイルから雑音を取り除き、目的の音声を取り出して聞き取るという内容です。
今回は、その時に用いたり作成したプログラムについて書いていこうかなと思います。
音声処理
音声処理と一言で言ってもいろいろあります。音量を上げ下げしたり、音の高低を変化させて音色を変化させたり、ボイスチェンジャーのように全く違う声に加工したり。
そういった音声処理を可能にしているのが、俗に言う「デジタルフィルタ」という技術です。
デジタルフィルタ
音声ファイルといっても結局は音のレベルが数値として記録されたデータです。それを、サンプリング周波数などといったパラメータに合わせて処理することにより、音声として再生し、聞き取るわけです。
デジタルフィルタの役目とは、これらの数値データに対して何らかの算術演算をし、その結果音声に変化を付ける、ということなのです。
フィルタの設計
フィルタの設計とはフィルタ係数の設定に他なりません。フィルタ係数を決めるためには、極(音声が増幅される点)と零点(音声が減衰される点)を設定します。
今回は専用のソフトウェアを使いました。平面上で自由に極や零点を配置するとそれに応じたフィルタ係数がわかるというものです。
このソフトウェアは実験の際に与えられたものなので、一応公開は控えておきます。上記のような仕様を持っているということだけ理解してもらえればと思います。
IIRフィルタ
IIRフィルタはインパルス応答が無限の長さを持っているので、Infinite Impulse Responseの略です。
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "inout16.c"
#define HEADER_SIZE 24
#define Fs 44100.0
#define Freq 1000.0
#define PI 3.141592653589793
int main(int argc, char *argv[]){
int i;
double x,y1,y2,y3,y4,y5,y6,y7,y8;
double r1 = -1.79259, s1 = 1.00000, p1 = 1.79959, q1 = -0.81873;
double r2 = -0.00750, s2 = 1.00000, p2 = 1.89007, q2 = -0.89346;
double r3 = -1.93206, s3 = 1.00000, p3 = 1.86878, q3 = -0.87858;
double r4 = -1.26445, s4 = 1.00000, p4 = 1.76765, q4 = -0.78232;
double b0 = 0.00100;
double buf_y1[2]={0.0,0.0};
double buf_y3[2]={0.0,0.0};
double buf_y5[2]={0.0,0.0};
double buf_y7[2]={0.0,0.0};
for(i=0; i<HEADER_SIZE; i++)
putchar(getchar());
while(get16d(stdin,&x)!=EOF){
y1=b0 * x+ p1 * buf_y1[0] + q1 * buf_y1[1];
y2=y1+ r1 * buf_y1[0] + s1 * buf_y1[1];
y3=y2+ p2 * buf_y3[0] + q2 * buf_y3[1];
y4=y3+ r2 * buf_y3[0] + s2 * buf_y3[1];
y5=y4+ p3 * buf_y5[0] + q3 * buf_y5[1];
y6=y5+ r3 * buf_y5[0] + s3 * buf_y5[1];
y7=y6+ p4 * buf_y7[0] + q4 * buf_y7[1];
y8=y7+ r4 * buf_y7[0] + s4 * buf_y7[1];
put16d(stdout,y8);
buf_y1[1]=buf_y1[0];
buf_y1[0]=y1;
buf_y3[1]=buf_y3[0];
buf_y3[0]=y3;
buf_y5[1]=buf_y5[0];
buf_y5[0]=y5;
buf_y7[1]=buf_y7[0];
buf_y7[0]=y7;
}
return 0;
}
inout16.cは実験の時に与えられた入出力用のマクロです。
今回は8次のフィルタを作成したので、フィルタ係数が過去と未来に8個ずつ、計16個用意されます。$b_0$は音声全体の大きさを決めます。音量調整みたいなものですね。
double r1 = -1.79259, s1 = 1.00000, p1 = 1.79959, q1 = -0.81873;
double r2 = -0.00750, s2 = 1.00000, p2 = 1.89007, q2 = -0.89346;
double r3 = -1.93206, s3 = 1.00000, p3 = 1.86878, q3 = -0.87858;
double r4 = -1.26445, s4 = 1.00000, p4 = 1.76765, q4 = -0.78232;
double b0 = 0.00100;
つぎに、信号を記録しておく場所を用意します。
double buf_y1[2]={0.0,0.0};
double buf_y3[2]={0.0,0.0};
double buf_y5[2]={0.0,0.0};
double buf_y7[2]={0.0,0.0};
IIRフィルタでは、過去の信号を用いて次の信号を計算します。そのため、過去の信号を記録しておく必要があります。
次に、フィルタ本体の処理です。
for(i=0; i<HEADER_SIZE; i++)
putchar(getchar());
while(get16d(stdin,&x)!=EOF){
y1=b0 * x+ p1 * buf_y1[0] + q1 * buf_y1[1];
y2=y1+ r1 * buf_y1[0] + s1 * buf_y1[1];
y3=y2+ p2 * buf_y3[0] + q2 * buf_y3[1];
y4=y3+ r2 * buf_y3[0] + s2 * buf_y3[1];
y5=y4+ p3 * buf_y5[0] + q3 * buf_y5[1];
y6=y5+ r3 * buf_y5[0] + s3 * buf_y5[1];
y7=y6+ p4 * buf_y7[0] + q4 * buf_y7[1];
y8=y7+ r4 * buf_y7[0] + s4 * buf_y7[1];
put16d(stdout,y8);
buf_y1[1]=buf_y1[0];
buf_y1[0]=y1;
buf_y3[1]=buf_y3[0];
buf_y3[0]=y3;
buf_y5[1]=buf_y5[0];
buf_y5[0]=y5;
buf_y7[1]=buf_y7[0];
buf_y7[0]=y7;
}
まずはヘッダ部分を読み飛ばします。次に、そして、順次音声データを読み込んで、各々のデータを用いて計算していきます。
y1=b0 * x+ p1 * buf_y1[0] + q1 * buf_y1[1];
y2=y1+ r1 * buf_y1[0] + s1 * buf_y1[1];
y3=y2+ p2 * buf_y3[0] + q2 * buf_y3[1];
y4=y3+ r2 * buf_y3[0] + s2 * buf_y3[1];
y5=y4+ p3 * buf_y5[0] + q3 * buf_y5[1];
y6=y5+ r3 * buf_y5[0] + s3 * buf_y5[1];
y7=y6+ p4 * buf_y7[0] + q4 * buf_y7[1];
y8=y7+ r4 * buf_y7[0] + s4 * buf_y7[1];
ここがデジタルフィルタの伝達関数部分の計算となります。各段階において、過去の信号にフィルた係数を掛け、足し合わせる形です。
そして、作成したデータを出力し、過去の信号を更新します。
put16d(stdout,y8);
buf_y1[1]=buf_y1[0];
buf_y1[0]=y1;
buf_y3[1]=buf_y3[0];
buf_y3[0]=y3;
buf_y5[1]=buf_y5[0];
buf_y5[0]=y5;
buf_y7[1]=buf_y7[0];
buf_y7[0]=y7;
以上がIIRフィルタのプログラムです。これをコンパイルして実行形式を作成し、音声ファイルを入力することで、フィルタを通した音声データを得ることができます。
実際にフィルタを通してみる
実際にフィルタを通してみます。まずは音声データを波形で見てみましょう。
まずはフィルタを通す前。
高周波数帯に酷い雑音が乗っています。もちろん再生したところで雑音以外何も聞こえません。2.5 kHz以下の低周波数帯には音声が含まれているわけですが、ここを取り出してやればいいわけです。ちなみに、寒色系部分ほど音が強いです。なので、この部分を弱め、暖色系部分を強めるようなフィルタを設計すればいいわけです。
ではフィルタをかけてみましょう。上記のCプログラムをコンパイルして実行形式を作成し、音声ファイルを入力します。
雑音がほとんど除去され、暖色系に弱められました。また、低周波帯の音声は増幅され、少し寒色系になっています。
こうして、音声が聞き取れるようになりました。
FIRフィルタ
FIRフィルタというと、インパルス応答が有限の長さを持っているわけです。Finite Impulse Responseの略ですね。
このフィルタを作る場合は、先ほどのようにして設定したフィルタ係数ではなく、要件を満たすようなインパルス応答を何らかの手段で入手し、それを用います。今回の要件は「12.5 kHz以上の高周波をカットする」となっています。
プログラムは以下のような感じです。
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "inout16.c"
#define HEADER_SIZE 24
#define Fs 44100.0
#define Freq 1000.0
#define PI 3.141592653589793
int main(int argc, char *argv[]){
int i,k,n;
int M=12;
double x,y;
double h_plus[13] = {
4.38019092118e-01,
3.09514166759e-01,
5.84022687533e-02,
-8.14914130799e-02,
-4.86437321301e-02,
2.79899005826e-02,
3.52868315211e-02,
-5.00321662753e-03,
-2.16854756524e-02,
-3.89065864425e-03,
1.19315459334e-02,
9.49228987543e-03,
1.10729254269e-03
};
double h_minus[12] = {
3.09514166759e-01,
5.84022687533e-02,
-8.14914130799e-02,
-4.86437321301e-02,
2.79899005826e-02,
3.52868315211e-02,
-5.00321662753e-03,
-2.16854756524e-02,
-3.89065864425e-03,
1.19315459334e-02,
9.49228987543e-03,
1.10729254269e-03,
};
int buf_size = 0;
double *buf_x = NULL;
for(i=0; i<HEADER_SIZE; i++)
putchar(getchar());
while(get16d(stdin,&x)!=EOF){
buf_size++;
buf_x = realloc(buf_x, buf_size * sizeof(double));
buf_x[buf_size - 1] = x;
}
for(n=0; n<buf_size; n++){
y = 0.0;
for(k=0; k<=M; k++){
if(n-k>=0)
y += h_plus[k] * buf_x[n-k];
}
for(k=0; k<M; k++){
if(n+k+1<buf_size)
y += h_minus[k] * buf_x[n+k+1];
}
put16d(stdout,y);
}
return 0;
}
for(i=0; i<HEADER_SIZE; i++)
putchar(getchar());
while(get16d(stdin,&x)!=EOF){
buf_size++;
buf_x = realloc(buf_x, buf_size * sizeof(double));
buf_x[buf_size - 1] = x;
}
この部分で、入力された音声データのヘッダ部分を読み飛ばし、データ部分を読み込んで配列に格納します。
for(n=0; n<buf_size; n++){
y = 0.0;
for(k=0; k<=M; k++){
if(n-k>=0)
y += h_plus[k] * buf_x[n-k];
}
for(k=0; k<M; k++){
if(n+k+1<buf_size)
y += h_minus[k] * buf_x[n+k+1];
}
put16d(stdout,y);
}
ここがフィルタ本体です。以下のFIRフィルタの式に基づいて、インパルス応答から次のデータを計算していきます。
$
y_n= h_0 x_n + \sum_{k=1}^{12}(b_k x_{n-k} + b_{-k} x_{k+n})
$
インパルス応答$h$たちはプログラム冒頭で配列に格納されています。こうした仕様にすることで、FIRフィルタを実装することができました。
このフィルタは高周波をカットするように設計されているので、以下のように雑音を低減します。除去前と比べてみてください。
12.5 kHz以上の高周波がカットされました。
おわりに
実際に雑音を処理することができたときは達成感がありました。フィルタの設計を粘った甲斐があったと思います。
フィルタ内での計算の説明はかなり端折りました。詳しい説明を知りたい人は別途調べてみてください。
それではまた。明日はgaeくんの記事です。よろしくお願いします！