イコライザ アルファ版

自分で作ったゲームや動画、面白いネタをみんなに宣伝しましょう!
また、気軽に作品の感想をコメントで残してあげて下さい。
tlnlniri
記事: 17
登録日時: 7年前

イコライザ アルファ版

#1

投稿記事 by tlnlniri » 7年前

日記からここに移動しました。

このページを見てくださった方に性能をテストしてほしいのです。
イコライザを起動して60FPS以上の状態でどれくらいCPUを使うのか調べて欲しいです。
タスクマネージャを使っただいたいの数値でかまいません。
環境は特に書いていただかなくてかまいません。

ご協力いただきありがとうございました。

※イコライザは同フォルダ内のwave.wavだけを起動時に読み込みます。
※FFT表示にも変更できません。

最終版で改善
※ステレオ(2ch)のみ対応してます。

イコライザ初期設定数値
ゲイン値すべて0dB

コード:

double eq_frequency[EQ_BANDS] = {31.25, 62.5, 125.0, 250.0, 500.0, 1000.0, 2000.0, 4000.0, 8000.0, 16000.0};	// 周波数テーブル
double eq_quality_factor[EQ_BANDS] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};	 // クオリティファクタ
double eq_dBgain[EQ_BANDS] = {15.0, 25.0, 24.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};	 // ゲイン[dB]
・バグ(正確には違う)
ヘッダのみのwavファイルを読み込むとCPU使用率が90%近くになる。

・5/29
アルファ2に更新 (飽きやすい性格のため、ベータにはいかないかと。)
グラフィックイコライザ追加
クオリティファクタ設定追加
最大振幅メーター追加(クリッピング検出用)
アルファ1の名残で起動時はwave.wavを読み込むが、ドラッグ&ドロップで別の.wavファイルを再生

・5/30
サンプル画像変更

・6/1
最終版アップロード(一番下にあります)
最後に編集したユーザー tlnlniri on 2013年8月06日(火) 07:56 [ 編集 10 回目 ]

nil
記事: 428
登録日時: 8年前

Re: イコライザ アルファ版

#2

投稿記事 by nil » 7年前

OS: Windows 7
CPU: Intel Pentium P6200 2.13GHz
問題なく60FPSを維持できています。
CPU使用率はおよそ20%です。

tlnlniri
記事: 17
登録日時: 7年前

Re: イコライザ アルファ版

#3

投稿記事 by tlnlniri » 7年前

>>涼雅さん
テストありがとうございます。
このような事に時間を使わせてしまって申し訳ありません。
お礼といっては何ですが、IIRの計算部分のソースコードでもそのうち載せようかと思っています。

アバター
みけCAT
記事: 6275
登録日時: 9年前
住所: 千葉県
連絡を取る:

Re: イコライザ アルファ版

#4

投稿記事 by みけCAT » 7年前

OS:Windows Vista SP2 32ビット
CPU: Intel(R) Core(TM)2 Duo T8100 @ 2.10GHz 2.10GHz
RAM: 4.00GB
CPU使用率:7~12%

UTAUで作成した音声をAudacityで限界まで増幅したデータと、単純なサイン波のデータ(1000Hz、0dB)を試しましたが、
どちらも線のグラフがほとんど動かず(水平のまま)、サイン波の音量は-42dB程度を示していました。
形式は44100Hz、16ビットです。
サイン波のデータを添付します。
添付ファイル
sign.zip
サイン波wavファイル
(11.55 KiB) ダウンロード数: 187 回
複雑な問題?マシンの性能を上げてOpenMPで殴ればいい!(死亡フラグ)

tlnlniri
記事: 17
登録日時: 7年前

Re: イコライザ アルファ版

#5

投稿記事 by tlnlniri » 7年前

>>みけCATさん

テストありがとうございます。また、貴重なデータありがとうございます。

・設定音量100%
・10バンドイコライザ全0dB
で再生させてもらいました。
(
 添付イコライザはクリッピングを防ぐため設定音量1%です。
 言い換えると1/100倍です。0.01倍です。
 そして http://www.geocities.jp/fkmtf928/dB_sound.html で見ると、
 0.01倍 → -40dB となっています。
 が、やはり正しい値と少しずれますね...。
)

-3dBと-4dBを行き来しました。
0dBに近いには近いですが、これは...すみませんわかりません。
 →離散化誤差、量子化誤差が関係してる?

また、同じ形の波形が続いているのでメーターは一定でいいと思います。
添付ファイル
WS000000.png
WS000000.png (50.26 KiB) 閲覧数: 6916 回
最後に編集したユーザー tlnlniri on 2013年5月02日(木) 18:55 [ 編集 4 回目 ]

tlnlniri
記事: 17
登録日時: 7年前

Re: イコライザ アルファ版

#6

投稿記事 by tlnlniri » 7年前

では、前に書いた通り、ソースを張ろうと思います。

ピーキングフィルタとかは
http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt
http://www.g200kg.com/jp/docs/makingvst/04.html
を参照してください。

コード:

/********************************************************************************
f0:中心周波数(遮断周波数)
fs:標本化周波数
Q:クオリティファクタ
g:ゲイン[dB]
a,b:フィルタ係数
********************************************************************************/
// ピーキングフィルタ
void BiquadPeakingFilter(double f0, double fs, double Q, double g, double a[], double b[])
{
	double A = pow(10.0, g / 40.0);
	double w0 = 2.0 * M_PI * f0 / fs;
	double alpha = sin(w0) / (2.0 * Q);
	b[0] = 1.0 + alpha * A;
	b[1] = -2.0 * cos(w0);
	b[2] = 1.0 - alpha * A;
	a[0] = 1.0 + alpha / A;
	a[1] = -2.0 * cos(w0);
	a[2] = 1.0 - alpha / A;
}
フィルタ係数の設定は

コード:

// フィルタ係数設定 (eq_sampling_rate = 44100.0)
for (int e = 0; e < EQ_BANDS; e++) {
	BiquadPeakingFilter(eq_frequency[e], eq_sampling_rate, eq_quality_factor[e], eq_dBgain[e], a[e], b[e]);
}
下記はコードの断片です。
『 参考程度 』に見てください。『 参考程度 』に。
もっと効率のいい方法があるはずです多分。
IIRを勉強すれば何やってるかわかると思います。

コード:

#define WAVE_CHANNELS				2		// ステレオのみ
#define SPK_BUFFER_LENGTH			512		// Left + Right
#define EQ_BANDS					10		// イコライザバンド数
#define IIR_WAVEFORM_LENGTH			(SPK_BUFFER_LENGTH / 2 + 2)	// IIR用波形長 ( / 2 = チャンネル分離)

double a[EQ_BANDS][3], b[EQ_BANDS][3];			// フィルタ係数
double y[2][EQ_BANDS + 1][IIR_WAVEFORM_LENGTH];	// 元波形(y[c][0]) & フィルタリング波形(y[c][1~EQ_BANDS]) (2 = Stereo)

// 遅延波形更新
for (int c = 0; c < WAVE_CHANNELS; c++) {
	for (int e = 0; e <= EQ_BANDS; e++) {
		for (int n = 0; n < 2; n++) {
			y[c][e][n] = y[c][e][IIR_WAVEFORM_LENGTH - 2 + n];
		}
	}
}

// 元波形更新
for (int c = 0; c < WAVE_CHANNELS; c++) {
	for (int n = 0; n < SPK_BUFFER_LENGTH / 2; n++) {
		y[c][0][2 + n] = (double)waveout_buffer[spk_during_playback][2 * n + c] * (volume / 100.0);		// チャンネル分離
	}
}

// フィルタリング (IIR)
for (int c = 0; c < WAVE_CHANNELS; c++) {
	for (int e = 0; e < EQ_BANDS; e++) {
		for (int n = 2; n < IIR_WAVEFORM_LENGTH; n++) {
			y[c][e + 1][n] = (b[e][0] / a[e][0]) * y[c][e][n] + (b[e][1] / a[e][0]) * y[c][e][n - 1] + (b[e][2] / a[e][0]) * y[c][e][n - 2] - (a[e][1] / a[e][0]) * y[c][e + 1][n - 1] - (a[e][2] / a[e][0]) * y[c][e + 1][n - 2];	// 遅延器数2で固定
		}
	}
}

// 再生バッファに書き込み
for (int c = 0; c < WAVE_CHANNELS; c++) {
	for (int n = 0; n < SPK_BUFFER_LENGTH / 2; n++) {
		waveout_buffer[spk_during_playback][2 * n + c] = (short)Limiter(y[c][EQ_BANDS][2 + n], 32767.0, -32768.0);		// ハードクリップ(無理やりshortの範囲に)
	}
}
何やってるか適当に書いておきます。
[0](元波形)を使って[1]を作成
[1]を使って[2]を作成
...
[EQ_BANDS - 1]を使って[EQ_BANDS]を作成
EQ_BANDSバンド数でいイコライズした波形が[EQ_BANDS] → 再生

IIR波形データ更新方法も適当に書いておきます。
[0]:
[■■]□□□□□□[□□]     最初の[■■]に[□□]をコピー
[□□]□□□□□□[□□]
[□□][新しい波形データ]  [□□]の後ろに新しい波形データ

[1 ~ EQ_BANDS]:
[■■]□□□□□□[□□]     最初の[■■]に[□□]をコピー
[□□]□□□□□□[□□]     後ろのデータは上書きされるから何もしなくていい



Q.何で全ソースコードを提示しないの?
A.オブジェクト指向じゃないし汚いし恥ずかしいから。
最後に編集したユーザー tlnlniri on 2013年5月06日(月) 22:59 [ 編集 1 回目 ]

tlnlniri
記事: 17
登録日時: 7年前

Re: イコライザ アルファ版

#7

投稿記事 by tlnlniri » 7年前

タイトルが変更できないのでここに書いておきます。

このイコライザは Biquad Equalizer と呼ばれる IIRフィルタ を用いたイコライザです。
位相ずれがあるようですが人間にはわからないようです。(私にはわかりません)

FFTを使ってもイコライザは作れますがそっちは Linear Phase Equalizer と呼ばれます。
数値制御するためには自由曲線を使うようです。
(lilith、ulilithがこちらに当てはまります。)
最後に編集したユーザー tlnlniri on 2013年7月02日(火) 01:19 [ 編集 1 回目 ]

tlnlniri
記事: 17
登録日時: 7年前

Re: イコライザ アルファ版

#8

投稿記事 by tlnlniri » 7年前

とりあえず更新はこれで終わりにします。

説明は面倒くさいので少し省略します。
右クリックで機能切り替えです。(FFTによってCPU使用率が跳ね上がります)
下にコントロール表示があります。
.wavファイルをD&Dして再生。

※初回起動時は必ずエラー表示が出ます。「設定ファイルがない」って。問題はありません。

動作にそれほど影響はないですが、低スペックPC(だけ?)だと読み込み遅延がまれに発生します。
なので みけCATさん の偉大なるプログラムを参考にさせてもらい、いつかは解消する予定です。
添付ファイル
WS000001.png
WS000000.png
WS000000.png (43.1 KiB) 閲覧数: 6627 回
Biquad Equalizer.zip
最終版
(980.5 KiB) ダウンロード数: 243 回

tlnlniri
記事: 17
登録日時: 7年前

Re: イコライザ アルファ版

#9

投稿記事 by tlnlniri » 7年前

今は多くの学生が夏休みですね。
「することが無くて暇だ!」って人は信号処理に挑戦してみてはどうでしょうか?
道具(一部)を置いておきますので。


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という立体音響化プログラムと同じアルゴリズムです。
添付ファイル
タイムストレッチ.c
(4.65 KiB) ダウンロード数: 206 回

返信

“作品お披露目掲示板” へ戻る