Cプログラミング言語を使用してFIRフィルターを実装する方法を教えてください。
3 に答える
FIR フィルターの設計は単純なトピックではありませんが、既に設計されたフィルターを実装することは (FIR 係数が既にあると仮定して) それほど悪くはありません。このアルゴリズムは畳み込みと呼ばれます。これが単純な実装です...
void convolve (double *p_coeffs, int p_coeffs_n,
double *p_in, double *p_out, int n)
{
int i, j, k;
double tmp;
for (k = 0; k < n; k++) // position in output
{
tmp = 0;
for (i = 0; i < p_coeffs_n; i++) // position in coefficients array
{
j = k - i; // position in input
if (j >= 0) // bounds check for input buffer
{
tmp += p_coeffs [i] * p_in [j];
}
}
p_out [k] = tmp;
}
}
基本的に、畳み込みは入力信号の加重移動平均を行います。重みはフィルター係数であり、合計すると 1.0 になると想定されます。重みの合計が 1.0 以外になる場合は、増幅/減衰とフィルタリングが行われます。
ところで - この関数は逆方向の係数配列を持っている可能性があります - 私は再確認していません。
特定のフィルターの FIR 係数を計算する方法については、その背後にかなりの量の数学があります。デジタル信号処理に関する優れた本が本当に必要です。これはPDFで無料で入手できますが、それがどれほど優れているかはわかりません. 私はRorabaughとOrfandisを持っています。どちらも 90 年代半ばに出版されましたが、これらは実際に時代遅れになることはありません。
複数のフィルタを組み合わせるには:
単位インパルス (最初の位置に 1 があり、それ以外は 0 の信号) から始めます。最初のフィルターを適用します。2 番目のフィルターを適用します。すべてのフィルターが適用されるまで続行します。結果は、結合されたフィルターが単位インパルスをどのように畳み込むかを示します (データが失われないように配列が十分に長い場合)。したがって、その中の値は、他のフィルターの合成である 1 つのフィルターの係数です。
サンプルコードは次のとおりです。
#include <stdio.h>
#include <string.h>
#define NumberOf(a) (sizeof (a) / sizeof *(a))
/* Convolve Signal with Filter.
Signal must contain OutputLength + FilterLength - 1 elements. Conversely,
if there are N elements in Signal, OutputLength may be at most
N+1-FilterLength.
*/
static void convolve(
float *Signal,
float *Filter, size_t FilterLength,
float *Output, size_t OutputLength)
{
for (size_t i = 0; i < OutputLength; ++i)
{
double sum = 0;
for (size_t j = 0; j < FilterLength; ++j)
sum += Signal[i+j] * Filter[FilterLength - 1 - j];
Output[i] = sum;
}
}
int main(void)
{
// Define a length for buffers that is long enough for this demonstration.
#define LongEnough 128
// Define some sample filters.
float Filter0[] = { 1, 2, -1 };
float Filter1[] = { 1, 5, 7, 5, 1 };
size_t Filter0Length = NumberOf(Filter0);
size_t Filter1Length = NumberOf(Filter1);
// Define a unit impulse positioned so it captures all of the filters.
size_t UnitImpulsePosition = Filter0Length - 1 + Filter1Length - 1;
float UnitImpulse[LongEnough];
memset(UnitImpulse, 0, sizeof UnitImpulse);
UnitImpulse[UnitImpulsePosition] = 1;
// Calculate a filter that is Filter0 and Filter1 combined.
float CombinedFilter[LongEnough];
// Set N to number of inputs that must be used.
size_t N = UnitImpulsePosition + 1 + Filter0Length - 1 + Filter1Length - 1;
// Subtract to find number of outputs of first convolution, then convolve.
N -= Filter0Length - 1;
convolve(UnitImpulse, Filter0, Filter0Length, CombinedFilter, N);
// Subtract to find number of outputs of second convolution, then convolve.
N -= Filter1Length - 1;
convolve(CombinedFilter, Filter1, Filter1Length, CombinedFilter, N);
// Remember size of resulting filter.
size_t CombinedFilterLength = N;
// Display filter.
for (size_t i = 0; i < CombinedFilterLength; ++i)
printf("CombinedFilter[%zu] = %g.\n", i, CombinedFilter[i]);
// Define two identical signals.
float Buffer0[LongEnough];
float Buffer1[LongEnough];
for (size_t i = 0; i < LongEnough; ++i)
{
Buffer0[i] = i;
Buffer1[i] = i;
}
// Convolve Buffer0 by using the two filters separately.
// Start with buffer length.
N = LongEnough;
// Subtract to find number of outputs of first convolution, then convolve.
N -= Filter0Length - 1;
convolve(Buffer0, Filter0, Filter0Length, Buffer0, N);
// Subtract to find number of outputs of second convolution, then convolve.
N -= Filter1Length - 1;
convolve(Buffer0, Filter1, Filter1Length, Buffer0, N);
// Remember the length of the result.
size_t ResultLength = N;
// Convolve Buffer1 with the combined filter.
convolve(Buffer1, CombinedFilter, CombinedFilterLength, Buffer1, ResultLength);
// Show the contents of Buffer0 and Buffer1, and their differences.
for (size_t i = 0; i < ResultLength; ++i)
{
printf("Buffer0[%zu] = %g. Buffer1[%zu] = %g. Difference = %g.\n",
i, Buffer0[i], i, Buffer1[i], Buffer0[i] - Buffer1[i]);
}
return 0;
}
このコード スニペットが機能しないことがわかりました (Visual Studio 2005)。
最終的に、畳み込みの質問に素晴らしい答えがあることがわかりました。
知らない人のために説明すると、畳み込みは FIR フィルタリングとまったく同じ操作です。「カーネル」は FIR フィルターのインパルス応答であり、信号は入力信号です。
これが、FIR コードを探していた貧弱な人に役立つことを願っています :-)