C语言实现fir1函数?

Question

matlab里面可以用fir1函数求得滤波系数，因为我要设计的滤波器是中心频率可变的，所以不同频率要求不同的系统，就需要在我的C程序里面有求系数的过程，也就是用C语言实现fir1谁能帮帮我啊~~~感激不尽FIR滤波器！

Answer 1

#include <stdio.h>
#ifdef WIN32
#include <conio.h>
#endif

#define SAMPLE double /* define the type used for data samples */

void clear(int ntaps, SAMPLE z[])
{
int ii;
for (ii = 0; ii < ntaps; ii++) {
z[ii] = 0;
}
}

SAMPLE fir_basic(SAMPLE input, int ntaps, const SAMPLE h[], SAMPLE z[])
{
int ii;
SAMPLE accum;

/* store input at the beginning of the delay line */
z[0] = input;

/* calc FIR */
accum = 0;
for (ii = 0; ii < ntaps; ii++) {
accum += h[ii] * z[ii];
}

/* shift delay line */
for (ii = ntaps - 2; ii >= 0; ii--) {
z[ii + 1] = z[ii];
}

return accum;
}

SAMPLE fir_circular(SAMPLE input, int ntaps, const SAMPLE h[], SAMPLE z[],
int *p_state)
{
int ii, state;
SAMPLE accum;

state = *p_state; /* copy the filter's state to a local */

/* store input at the beginning of the delay line */
z[state] = input;
if (++state >= ntaps) { /* incr state and check for wrap */
state = 0;
}

/* calc FIR and shift data */
accum = 0;
for (ii = ntaps - 1; ii >= 0; ii--) {
accum += h[ii] * z[state];
if (++state >= ntaps) { /* incr state and check for wrap */
state = 0;
}
}

*p_state = state; /* return new state to caller */

return accum;
}

SAMPLE fir_shuffle(SAMPLE input, int ntaps, const SAMPLE h[], SAMPLE z[])
{
int ii;
SAMPLE accum;

/* store input at the beginning of the delay line */
z[0] = input;

/* calc FIR and shift data */
accum = h[ntaps - 1] * z[ntaps - 1];
for (ii = ntaps - 2; ii >= 0; ii--) {
accum += h[ii] * z[ii];
z[ii + 1] = z[ii];
}

return accum;
}

SAMPLE fir_split(SAMPLE input, int ntaps, const SAMPLE h[], SAMPLE z[],
int *p_state)
{
int ii, end_ntaps, state = *p_state;
SAMPLE accum;
SAMPLE const *p_h;
SAMPLE *p_z;

/* setup the filter */
accum = 0;
p_h = h;

/* calculate the end part */
p_z = z + state;
*p_z = input;
end_ntaps = ntaps - state;
for (ii = 0; ii < end_ntaps; ii++) {
accum += *p_h++ * *p_z++;
}

/* calculate the beginning part */
p_z = z;
for (ii = 0; ii < state; ii++) {
accum += *p_h++ * *p_z++;
}

/* decrement the state, wrapping if below zero */
if (--state < 0) {
state += ntaps;
}
*p_state = state; /* return new state to caller */

return accum;
}

SAMPLE fir_double_z(SAMPLE input, int ntaps, const SAMPLE h[], SAMPLE z[],
int *p_state)
{
SAMPLE accum;
int ii, state = *p_state;
SAMPLE const *p_h, *p_z;

/* store input at the beginning of the delay line as well as ntaps more */
z[state] = z[state + ntaps] = input;

/* calculate the filter */
p_h = h;
p_z = z + state;
accum = 0;
for (ii = 0; ii < ntaps; ii++) {
accum += *p_h++ * *p_z++;
}

/* decrement state, wrapping if below zero */
if (--state < 0) {
state += ntaps;
}
*p_state = state; /* return new state to caller */

return accum;
}

SAMPLE fir_double_h(SAMPLE input, int ntaps, const SAMPLE h[], SAMPLE z[],
int *p_state)
{
SAMPLE accum;
int ii, state = *p_state;
SAMPLE const *p_h, *p_z;

/* store input at the beginning of the delay line */
z[state] = input;

/* calculate the filter */
p_h = h + ntaps - state;
p_z = z;
accum = 0;
for (ii = 0; ii < ntaps; ii++) {
accum += *p_h++ * *p_z++;
}

/* decrement state, wrapping if below zero */
if (--state < 0) {
state += ntaps;
}
*p_state = state; /* return new state to caller */

return accum;
}

int main(void)
{
#define NTAPS 6
static const SAMPLE h[NTAPS] = { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 };
static SAMPLE h2[2 * NTAPS];
static SAMPLE z[2 * NTAPS];
#define IMP_SIZE (3 * NTAPS)
static SAMPLE imp[IMP_SIZE];
SAMPLE output;
int ii, state;

/* make impulse input signal */
clear(IMP_SIZE, imp);
imp[5] = 1.0;

/* create a SAMPLEd h */
for (ii = 0; ii < NTAPS; ii++) {
h2[ii] = h2[ii + NTAPS] = h[ii];
}

/* test FIR algorithms */

printf("Testing fir_basic:\n ");
clear(NTAPS, z);
for (ii = 0; ii < IMP_SIZE; ii++) {
output = fir_basic(imp[ii], NTAPS, h, z);
printf("%3.1lf ", (double) output);
}
printf("\n\n");

printf("Testing fir_shuffle:\n ");
clear(NTAPS, z);
state = 0;
for (ii = 0; ii < IMP_SIZE; ii++) {
output = fir_shuffle(imp[ii], NTAPS, h, z);
printf("%3.1lf ", (double) output);
}
printf("\n\n");

printf("Testing fir_circular:\n ");
clear(NTAPS, z);
state = 0;
for (ii = 0; ii < IMP_SIZE; ii++) {
output = fir_circular(imp[ii], NTAPS, h, z, &state);
printf("%3.1lf ", (double) output);
}
printf("\n\n");

printf("Testing fir_split:\n ");
clear(NTAPS, z);
state = 0;
for (ii = 0; ii < IMP_SIZE; ii++) {
output = fir_split(imp[ii], NTAPS, h, z, &state);
printf("%3.1lf ", (double) output);
}
printf("\n\n");

printf("Testing fir_double_z:\n ");
clear(2 * NTAPS, z);
state = 0;
for (ii = 0; ii < IMP_SIZE; ii++) {
output = fir_double_z(imp[ii], NTAPS, h, z, &state);
printf("%3.1lf ", (double) output);
}
printf("\n\n");

printf("Testing fir_double_h:\n ");
clear(NTAPS, z);
state = 0;
for (ii = 0; ii < IMP_SIZE; ii++) {
output = fir_double_h(imp[ii], NTAPS, h2, z, &state);
printf("%3.1lf ", (double) output);
}

#ifdef WIN32
printf("\n\nHit any key to continue.");
getch();
#endif

return 0;
}

1. fir_basic: 实现基本的FIR滤波器

2. fir_circular: 说明环行buffer是如何实现FIR的。

3. fir_shuffle: 一些TI的处理器上使用的shuffle down技巧

4. fir_split: 把FIR滤波器展开为两块，避免使用环行缓存。

5. fir_double_z: 使用双精度的延迟线，使可以使用一个flat buffer。

6. fir_double_h: 使用双精度的系数，使可以使用一个flat buffer。