今は多くの学生が夏休みですね。
「することが無くて暇だ！」って人は信号処理に挑戦してみてはどうでしょうか？
道具(一部)を置いておきますので。
8/14 : コードの間違いを修正
・IIRでコンパイルできないバグがあったので修正
IIR
コード:
// IIR差分方程式
void IIR_Filtering(double y[], double x[], int y_length, double a[], double b[])
{
int n;
for (n = 2; n < y_length; n++) {
y[n] = b[0] * x[n] + b[1] * x[n - 1] + b[2] * x[n - 2] - a[1] * y[n - 1] - a[2] * y[n - 2];
}
}
// 全域通過フィルタ
void IIR_AllPassFilter(double fc, double fs, double Q, double a[], double b[])
{
double w0 = 2.0 * M_PI * fc / fs;
double alpha = sin(w0) / (2.0 * Q);
a[0] = 1.0 + alpha;
a[1] = (-2.0 * cos(w0)) / a[0];
a[2] = (1.0 - alpha) / a[0];
b[0] = (1.0 - alpha) / a[0];
b[1] = (-2.0 * cos(w0)) / a[0];
b[2] = (1.0 + alpha) / a[0];
a[0] = 1.0;
}
// 低域通過フィルタ
void IIR_LowPassFilter(double fc, double fs, double Q, double a[], double b[])
{
double w0 = 2.0 * M_PI * fc / fs;
double alpha = sin(w0) / (2.0 * Q);
a[0] = 1.0 + alpha;
a[1] = (-2.0 * cos(w0)) / a[0];
a[2] = (1.0 - alpha) / a[0];
b[0] = ((1.0 - cos(w0)) / 2.0) / a[0];
b[1] = (1.0 - cos(w0)) / a[0];
b[2] = ((1.0 - cos(w0)) / 2.0) / a[0];
a[0] = 1.0;
}
// 高域通過フィルタ
void IIR_HighPassFilter(double fc, double fs, double Q, double a[], double b[])
{
double w0 = 2.0 * M_PI * fc / fs;
double alpha = sin(w0) / (2.0 * Q);
a[0] = 1.0 + alpha;
a[1] = (-2.0 * cos(w0)) / a[0];
a[2] = (1.0 - alpha) / a[0];
b[0] = ((1.0 + cos(w0)) / 2.0) / a[0];
b[1] = (-(1.0 + cos(w0))) / a[0];
b[2] = ((1.0 + cos(w0)) / 2.0) / a[0];
a[0] = 1.0;
}
#ifdef USE_BPF_CONSTANT_SKIRT_GAIN
// 帯域通過フィルタ（constant skirt gain, peak gain = Q）
void IIR_BandPassFilter1(double fc, double fs, double Q, double a[], double b[])
{
double w0 = 2.0 * M_PI * fc / fs;
double alpha = sin(w0) / (2.0 * Q);
a[0] = 1.0 + alpha;
a[1] = (-2.0 * cos(w0)) / a[0];
a[2] = (1.0 - alpha) / a[0];
b[0] = (sin(w0) / 2.0) / a[0]; // = (Q * alpha) / a[0];
b[1] = (0.0) / a[0];
b[2] = (-sin(w0) / 2.0) / a[0]; // = (-Q * alpha) / a[0];
a[0] = 1.0;
}
// 帯域通過フィルタ（constant skirt gain, peak gain = Q）
void IIR_BandPassFilter2(double fc1, double fc2, double fs, double a[], double b[])
{
double fc = sqrt(fc1 * fc2);
double Q = fc / (fc2 - fc1);
double w0 = 2.0 * M_PI * fc / fs;
double alpha = sin(w0) / (2.0 * Q);
a[0] = 1.0 + alpha;
a[1] = (-2.0 * cos(w0)) / a[0];
a[2] = (1.0 - alpha) / a[0];
b[0] = (sin(w0) / 2.0) / a[0]; // = (Q * alpha) / a[0];
b[1] = (0.0) / a[0];
b[2] = (-sin(w0) / 2.0) / a[0]; // = (-Q * alpha) / a[0];
a[0] = 1.0;
}
#else
// 帯域通過フィルタ（constant 0 dB peak gain）
void IIR_BandPassFilter1(double fc, double fs, double Q, double a[], double b[])
{
double w0 = 2.0 * M_PI * fc / fs;
double alpha = sin(w0) / (2.0 * Q);
a[0] = 1.0 + alpha;
a[1] = (-2.0 * cos(w0)) / a[0];
a[2] = (1.0 - alpha) / a[0];
b[0] = (alpha) / a[0];
b[1] = (0.0) / a[0];
b[2] = (-alpha) / a[0];
a[0] = 1.0;
}
// 帯域通過フィルタ（constant 0 dB peak gain）
void IIR_BandPassFilter2(double fc1, double fc2, double fs, double a[], double b[])
{
double fc = sqrt(fc1 * fc2);
double Q = fc / (fc2 - fc1);
double w0 = 2.0 * M_PI * fc / fs;
double alpha = sin(w0) / (2.0 * Q);
a[0] = 1.0 + alpha;
a[1] = (-2.0 * cos(w0)) / a[0];
a[2] = (1.0 - alpha) / a[0];
b[0] = (alpha) / a[0];
b[1] = (0.0) / a[0];
b[2] = (-alpha) / a[0];
a[0] = 1.0;
}
#endif
// 帯域阻止フィルタ
void IIR_NotchFilter(double fc, double fs, double Q, double a[], double b[])
{
double w0 = 2.0 * M_PI * fc / fs;
double alpha = sin(w0) / (2.0 * Q);
a[0] = 1.0 + alpha;
a[1] = (-2.0 * cos(w0)) / a[0];
a[2] = (1.0 - alpha) / a[0];
b[0] = (1.0) / a[0];
b[1] = (-2.0 * cos(w0)) / a[0];
b[2] = (1.0) / a[0];
a[0] = 1.0;
}
// 低域増幅フィルタ
void IIR_LowShelvingFilter(double fc, double fs, double Q, double g, double a[], double b[])
{
double A = pow(10.0, g / 40.0);
double w0 = 2.0 * M_PI * fc / fs;
double alpha = sin(w0) / (2.0 * Q);
a[0] = (A + 1.0) + (A - 1.0) * cos(w0) + 2.0 * sqrt(A) * alpha;
a[1] = (-2.0 * ((A - 1.0) + (A + 1.0) * cos(w0))) / a[0];
a[2] = ((A + 1.0) + (A - 1.0) * cos(w0) - 2.0 * sqrt(A) * alpha) / a[0];
b[0] = (A * ((A + 1.0) - (A - 1.0) * cos(w0) + 2.0 * sqrt(A) * alpha)) / a[0];
b[1] = (2.0 * A * ((A - 1.0) - (A + 1.0) * cos(w0))) / a[0];
b[2] = (A * ((A + 1.0) - (A - 1.0) * cos(w0) - 2.0 * sqrt(A) * alpha)) / a[0];
a[0] = 1.0;
}
// 高域増幅フィルタ
void IIR_HighShelvingFilter(double fc, double fs, double Q, double g, double a[], double b[])
{
double A = pow(10.0, g / 40.0);
double w0 = 2.0 * M_PI * fc / fs;
double alpha = sin(w0) / (2.0 * Q);
a[0] = (A + 1.0) - (A - 1.0) * cos(w0) + 2.0 * sqrt(A) * alpha;
a[1] = (2.0 * ((A - 1.0) - (A + 1.0) * cos(w0))) / a[0];
a[2] = ((A + 1.0) - (A - 1.0) * cos(w0) - 2.0 * sqrt(A) * alpha) / a[0];
b[0] = (A * ((A + 1.0) + (A - 1.0) * cos(w0) + 2.0 * sqrt(A) * alpha)) / a[0];
b[1] = (-2.0 * A * ((A - 1.0) + (A + 1.0) * cos(w0))) / a[0];
b[2] = (A * ((A + 1.0) + (A - 1.0) * cos(w0) - 2.0 * sqrt(A) * alpha)) / a[0];
a[0] = 1.0;
}
// 帯域増幅フィルタ
void IIR_PeakingFilter(double fc, double fs, double Q, double g, double a[], double b[])
{
double A = pow(10.0, g / 40.0);
double w0 = 2.0 * M_PI * fc / fs;
double alpha = sin(w0) / (2.0 * Q);
a[0] = 1.0 + alpha / A;
a[1] = (-2.0 * cos(w0)) / a[0];
a[2] = (1.0 - alpha / A) / a[0];
b[0] = (1.0 + alpha * A) / a[0];
b[1] = (-2.0 * cos(w0)) / a[0];
b[2] = (1.0 - alpha * A) / a[0];
a[0] = 1.0;
}
Q(クオリティファクタ) = 1.0 / sqrt(2.0);
でバターワース特性(遮断周波数で-3.0dB)です。
FIR ※(J + 1 - 1)は私にとってわかりやすいようにしてあるだけです。
コード:
// sinc(x)
double sinc(double x)
{
if (x == 0.0) {
return 1.0;
} else {
return (sin(x) / x);
}
}
// FIR差分方程式
void FIR_Filtering(double y[], double x[], int y_length, double b[], int J)
{
int n, m;
for (n = 0; n < y_length; n++) {
y[n] = 0.0;
for (m = 0; m <= J; m++) {
y[n] += b[m] * x[J + n - m];
}
}
}
// FIR遅延器数を取得 (delta:遷移帯域幅[Hz], fs:標本化周波数[Hz])
int FIR_GetTheNumberOfDelayUnits(double delta, double fs)
{
int J = (int)(3.1 / (delta / fs) + 0.5) - 1; // 遅延器数
if (J % 2 == 1) {
J++; // 偶数にする
}
return J;
}
// FIR遅延配列を確保
double *FIR_MemoryAllocationOfDelayUnits(int J)
{
int i;
double *b;
b = (double *)calloc(J + 1, sizeof(double));
for (i = 0; i <= J; i++) { // memset使えよってツッコミはなしで
b[i] = 0.0;
}
return b;
}
// 低域通過フィルタ
void FIR_LowPassFilter(double fe, double fs, int J, double b[])
{
int m;
int offset;
fe = fe / fs;
offset = J / 2;
for (m = -J / 2; m <= J / 2; m++) {
b[offset + m] = 2.0 * fe * sinc(2.0 * M_PI * fe * m);
}
for (m = 0; m < J + 1; m++) {
b[m] *= 0.5 - 0.5 * cos(2.0 * M_PI * (double)m / (double)(J + 1 - 1));
}
}
// 高域通過フィルタ
void FIR_HighPassFilter(double fe, double fs, int J, double b[])
{
int m;
int offset;
fe = fe / fs;
offset = J / 2;
for (m = -J / 2; m <= J / 2; m++) {
b[offset + m] = sinc(M_PI * m) - 2.0 * fe * sinc(2.0 * M_PI * fe * m);
}
for (m = 0; m < J + 1; m++) {
b[m] *= 0.5 - 0.5 * cos(2.0 * M_PI * (double)m / (double)(J + 1 - 1));
}
}
// 帯域通過フィルタ
void FIR_BandPassFilter(double fe1, double fe2, double fs, int J, double b[])
{
int m;
int offset;
fe1 = fe1 / fs;
fe2 = fe2 / fs;
offset = J / 2;
for (m = -J / 2; m <= J / 2; m++) {
b[offset + m] = 2.0 * fe2 * sinc(2.0 * M_PI * fe2 * m) - 2.0 * fe1 * sinc(2.0 * M_PI * fe1 * m);
}
for (m = 0; m < J + 1; m++) {
b[m] *= 0.5 - 0.5 * cos(2.0 * M_PI * (double)m / (double)(J + 1 - 1));
}
}
// 帯域阻止フィルタ
void FIR_BandEliminationFilter(double fe1, double fe2, double fs, int J, double b[])
{
int m;
int offset;
fe1 = fe1 / fs;
fe2 = fe2 / fs;
offset = J / 2;
for (m = -J / 2; m <= J / 2; m++) {
b[offset + m] = sinc(M_PI * m) - 2.0 * fe2 * sinc(2.0 * M_PI * fe2 * m) + 2.0 * fe1 * sinc(2.0 * M_PI * fe1 * m);
}
for (m = 0; m < J + 1; m++) {
b[m] *= 0.5 - 0.5 * cos(2.0 * M_PI * (double)m / (double)(J + 1 - 1));
}
}
// 低域増幅フィルタ
void FIR_LowShelvingFilter(double fe, double fs, int J, double b[], double g)
{
int m;
int offset;
g = pow(10.0, g / 20.0);
fe = fe / fs;
offset = J / 2;
for (m = -J / 2; m <= J / 2; m++) {
b[offset + m] = (2.0 * fe * sinc(2.0 * M_PI * fe * m)) * g;
}
for (m = 0; m < J + 1; m++) {
b[m] *= 0.5 - 0.5 * cos(2.0 * M_PI * (double)m / (double)(J + 1 - 1));
}
}
// 高域増幅フィルタ
void FIR_HighShelvingFilter(double fe, double fs, int J, double b[], double g)
{
int m;
int offset;
g = pow(10.0, g / 20.0);
fe = fe / fs;
offset = J / 2;
for (m = -J / 2; m <= J / 2; m++) {
b[offset + m] = (sinc(M_PI * m) - 2.0 * fe * sinc(2.0 * M_PI * fe * m)) * g;
}
for (m = 0; m < J + 1; m++) {
b[m] *= 0.5 - 0.5 * cos(2.0 * M_PI * (double)m / (double)(J + 1 - 1));
}
}
// 帯域増幅フィルタ
void FIR_PeakingFilter(double fe1, double fe2, double fs, int J, double b[], double g)
{
int m;
int offset;
g = pow(10.0, g / 20.0);
fe1 = fe1 / fs;
fe2 = fe2 / fs;
offset = J / 2;
for (m = -J / 2; m <= J / 2; m++) {
b[offset + m] =
2.0 * fe1 * sinc(2.0 * M_PI * fe1 * m)
+ (2.0 * fe2 * sinc(2.0 * M_PI * fe2 * m) - 2.0 * fe1 * sinc(2.0 * M_PI * fe1 * m)) * g
+ sinc(M_PI * m) - 2.0 * fe2 * sinc(2.0 * M_PI * fe2 * m);
}
for (m = 0; m < J + 1; m++) {
b[m] *= 0.5 - 0.5 * cos(2.0 * M_PI * (double)m / (double)(J + 1 - 1));
}
}
基本的にIIRの使い方とFIRの使い方は同じです。
配列の確保
・IIR
y[2 + 波形長] ：出力
x[2 + 波形長] ：入力
・FIR
y[波形長] ：出力
x[J(=FIR_GetTheNumberOfDelayUnits) + 波形長] ：入力
使用例 (偶に追加します)
・walkmanのカラオケ機能:IIR,FIR
普通は(左-右)でカラオケにすると音像が中央部に位置する低音が消えるのですが、元の波形に左右別々にローパスフィルタを通して加算合成してあげると、walkmanのようなカラオケが作れます。
※カラオケの計算方法は(右-左)でもいい。
※walkmanでも(左-右)波形はモノラルになっている。
※完全に低域だけを通過させるわけではないのでボーカルが少し残りますが、walkmanでも共通です。
「左カラオケ = 左元波形をローパスに通した波形 + 左 - 右」
「右カラオケ = 右元波形をローパスに通した波形 + 左 - 右」
・walkmanのVPT(サラウンド)機能：FIR
公式サイトを見ると、インパルス応答を高速フーリエ変換して波形に周波数領域で乗算しているようです。
インパルス応答をFIR_Filtering()で畳み込んでも同じです。
インパルス応答の畳み込みに関するサイト「
http://www.ari-web.com/service/soft/reverb-4.htm」
・イコライザ:IIR,FIR
ピーキングフィルタをバンド数分用意して順に通してください。
・FFTを使わない周波数解析:IIR,FIR
バンドパスフィルタを沢山使えば実現できます。サンプル数が少なくても周波数分解能は一定なので周波数解析で楽譜を作りたいならこちらの方がよいかもしれません。
2013/9/28 追加
※もし間違っていたら指摘していただけるとありがたいです
・フィルタ特性
IIR振幅特性
コード:
// IIR振幅特性（対数軸）
double IIR_FrequencyCharacteristics_Logarithm(double fc, double fs, double a[], double b[])
{
double w = 2.0 * M_PI * fc / fs;
double process_a = b[0] + b[1] * cos(-w) + b[2] * cos(-2.0 * w);
double process_b = b[1] * sin(-w) + b[2] * sin(-2.0 * w);
double process_c = 1.0 + a[1] * cos(-w) + a[2] * cos(-2.0 * w);
double process_d = a[1] * sin(-w) + a[2] * sin(-2.0 * w);
double real = (process_a * process_c + process_b * process_d) / (process_c * process_c + process_d * process_d);
double imag = (process_b * process_c - process_a * process_d) / (process_c * process_c + process_d * process_d);
return (10.0 * log10(real * real + imag * imag));
}
IIR位相特性
コード:
// IIR位相特性
double IIR_PhaseCharacteristics(double fc, double fs, double a[], double b[])
{
double w = 2.0 * M_PI * fc / fs;
double process_a = b[0] + b[1] * cos(-w) + b[2] * cos(-2.0 * w);
double process_b = b[1] * sin(-w) + b[2] * sin(-2.0 * w);
double process_c = 1.0 + a[1] * cos(-w) + a[2] * cos(-2.0 * w);
double process_d = a[1] * sin(-w) + a[2] * sin(-2.0 * w);
double real = (process_a * process_c + process_b * process_d) / (process_c * process_c + process_d * process_d);
double imag = (process_b * process_c - process_a * process_d) / (process_c * process_c + process_d * process_d);
return atan2(imag, real);
}
FIR振幅特性
コード:
// FIR振幅特性（対数軸）
double FIR_FrequencyCharacteristics_Logarithm(double fc, double fs, double b[], int J)
{
double w = 2.0 * M_PI * fc / fs;
double real = 0.0;
double imag = 0.0;
for (int i = 0; i <= J; i++) {
real += b[i] * cos(-i * w);
imag += b[i] * sin(-i * w);
}
return (10.0 * log10(real * real + imag * imag));
}
FIR位相特性
コード:
// FIR位相特性
double FIR_PhaseCharacteristics(double fc, double fs, double b[], int J)
{
double w = 2.0 * M_PI * fc / fs;
double real = 0.0;
double imag = 0.0;
for (int i = 0; i <= J; i++) {
real += b[i] * cos(-i * w);
imag += b[i] * sin(-i * w);
}
return atan2(imag, real);
}
・IIRを例に使い方の説明
標本化周波数44100Hzなので、対象となる周波数は22050Hzまでです。
コード:
// フィルタ係数の周波数特性を対数値で取得
for (double fc = 0.0; fc < 44100.0 / 2; fc++) {
IIR_FrequencyCharacteristics_Logarithm(fc, 44100.0, a, b);
}
次に時間が空いて気が乗ったらタイムストレッチの方法でも書きたいと思います。
使用アルゴリズム：PSOLA法（自己相関関数を使用）
ざっくりとやり方を。
自己相関関数で波形の周期を算出し、周期分の波形を足したり引いたりして、つなぎ目が自然になるように重畳加算する。
・重畳加算
｜＼ ＋ ／｜ ＝ ｜￣｜
・早送り
｜周期｜｜周期｜｜周期｜
｜＼ ／｜
↓
｜周期｜｜周期｜
｜￣｜
・スロー再生
｜周期｜｜周期｜
／｜ ｜＼
↓
｜周期｜｜周期｜｜周期｜
｜￣｜
※一応ソースコードは作りましたが、あまりにも汚いのでいつアップロードできるかわかりません。一言声をかけていただければお送りしますが。
調整めんどくさいんでそのままアップします
2chに対応させるにはひと工夫必要です
読み込む量が左右で一定ではないので
・標本化周波数44100Hzのwavを標本化周波数を変えずに88200Hzで再生する方法
1秒間に44100個のデータが必要だったのが、1秒間に88200個に変わるということです。
しかし標本化周波数44100Hzなのでデータ数は変えられませんので無理やり88200個詰めてください。
例えば
short data_44100[44100] = {0}, data_88200[88200] = {0};
for (int i = 0; i < 44100; i++) {
data_44100
= data_88200[i * 2];
}
現在、ステレオのままカラオケにする方法を思いついたので実験中・・・
→停滞中
----------------------------------------------------------------------------------------
以前私がどこかで質問した、
「周波数領域で振幅を加工したらIFFTしたときに両端が0にならなかった」
原因は巡回畳み込みにあります。
理想周波数特性のフィルタを掛けてたせいです。
インパルス応答のデータ数が無限長（めっちゃでかい）になってしまい、サンプル波形長＋インパルス応答長＞FFTサンプル数 となってしまいました。
-----------------------------------------------------------------------------------------
2014/1/3
ついに念願の立体音響の実装に成功しました！！
頭部伝達関数を使っています。
Bienという立体音響化プログラムと同じアルゴリズムです。
メニュー
カレンダー
Clock
最新メンバー
投稿数
統計データ
Link to us
カウンター