logo资料库

FIR低通滤波器C语言实现.pdf

第1页 / 共9页
第2页 / 共9页
第3页 / 共9页
第4页 / 共9页
第5页 / 共9页
第6页 / 共9页
第7页 / 共9页
第8页 / 共9页
资料共9页,剩余部分请下载后查看
文文: C:\Documents and Settings\user\桌桌\FirAlgs.c 2009-3-27 27, 15:53:51 /************************************************************* *************** * * Name: FirAlgs.c * * Synopsis: FIR filter algorithms for use in C * * Description: * * This module provides functions to implement Finite Impulse Response (FIR) * filters using a variety of algorithms. * * These functions use most or all of the following input parameters: * * input - the input sample data * ntaps - the number of taps in the filter * h[] - the FIR coefficient array * z[] - the FIR delay line array * *p_state - the "state" of the filter; initialize this to 0 before * the first call * * The functions fir_basic, fir_shuffle, and fir_circular are not very * efficient in C, and are provided here mainly to illustrate FIR * algorithmic concepts. However, the functions fir_split, fir_double_z * and fir_double_h are all fairly efficient ways to implement FIR filters * in C. * * Copyright 2000 by Grant R. Griffin * * Thanks go to contributors of comp.dsp for teaching me some of these * techniques, and to Jim Thomas for his review and great suggestions. * * The Wide Open License (WOL) * * Permission to use, copy, modify, distribute and sell this software and its * documentation for any purpose is hereby granted without fee, provided that * the above copyright notice and this license appear in all 第 1 页
文文: C:\Documents and Settings\user\桌桌\FirAlgs.c 2009-3-27 27, 15:53:51 source copies. * THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT EXPRESS OR IMPLIED WARRANTY OF * ANY KIND. See http://www.dspguru.com/wol.htm for more information. * ************************************************************** ***************/ #include #ifdef WIN32 #include #endif #define SAMPLE double /* define the type used for data samples */ /************************************************************* ***************/ /* clear: zeroize a FIR delay line */ /************************************************************* ***************/ void clear(int ntaps, SAMPLE z[]) { int ii; for (ii = 0; ii < ntaps; ii++) { z[ii] = 0; } } /************************************************************* *************** * fir_basic: Does the basic FIR algorithm: store input sample, calculate * output sample, move delay line ************************************************************** ***************/ 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; 第 2 页
文文: C:\Documents and Settings\user\桌桌\FirAlgs.c 2009-3-27 27, 15:53:51 /* 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; } /************************************************************* *************** * fir_circular: This function illustrates the use of "circular" buffers * in FIR implementations. The advantage of circular buffers is that they * alleviate the need to move data samples around in the delay line (as * was done in all the functions above). Most DSP microprocessors implement * circular buffers in hardware, which allows a single FIR tap can be * calculated in a single instruction. That works fine when programming in * assembly, but since C doesn't have any constructs to represent circular * buffers, you need to "fake" them in C by adding an extra "if" statement * inside the FIR calculation, even if the DSP processor provides hardware to * implement them without overhead. ************************************************************** ***************/ 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 */ 第 3 页
文文: C:\Documents and Settings\user\桌桌\FirAlgs.c 2009-3-27 27, 15:53:51 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; } /************************************************************* *************** * fir_shuffle: This is like fir_basic, except that data is shuffled by * moving it _inside_ the calculation loop. This is similar to the MACD * instruction on TI's fixed-point processors ************************************************************** ***************/ 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; 第 4 页
文文: C:\Documents and Settings\user\桌桌\FirAlgs.c 2009-3-27 27, 15:53:51 } /************************************************************* *************** * fir_split: This splits the calculation into two parts so the circular * buffer logic doesn't have to be done inside the calculation loop ************************************************************** ***************/ 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; } /************************************************************* *************** 第 5 页
文文: C:\Documents and Settings\user\桌桌\FirAlgs.c 2009-3-27 27, 15:53:51 * fir_double_z: double the delay line so the FIR calculation always * operates on a flat buffer ************************************************************** ***************/ 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; } /************************************************************* *************** * fir_double_h: uses doubled coefficients (supplied by caller) so that the * filter calculation always operates on a flat buffer. ************************************************************** ***************/ 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; 第 6 页
文文: C:\Documents and Settings\user\桌桌\FirAlgs.c 2009-3-27 27, 15:53:51 /* 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; } /************************************************************* *************** * main: This provides a simple test suite for the functions above. An * impulse is fed into each filter implementation, so the output should be * the "impulse response", that is, the coefficients of the filter. You * should see some zeroes, followed by the "coefficents" below, followed by a * few more zeroes. ************************************************************** ***************/ 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); 第 7 页
文文: C:\Documents and Settings\user\桌桌\FirAlgs.c 2009-3-27 27, 15:53:51 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++) { 第 8 页
分享到:
收藏