FIR Filter Implementation in MATLAB and in CI. FIR Digital Filter Design and Implementation in MATLAB
The MATLAB function, or “M-file”, shown in Listing 1 below provides an example of how to design and use finite impulse response (FIR) filters using the FIR filter design function fir1(N,Wn), the digital filter response calculating function freqz(A,B), and the digital filtering function filter( B,A,X) which are provided in the “MATLAB Signal Processing Toolbox”, which is described at
(http://www.mathworks.com/).
Make sure that you completely understand this M-file. MATLAB provides the following “help” documentation for these three routines: (For example, you may obtain information on the FIR1 function when running MATLAB by typing “help FIR1 [Enter]” at the MATLAB prompt):
_____________________________________________________________________________________
FIR1 FIR filter design (using the window method) function
B = FIR1(N,Wn) designs an N'th order lowpass FIR digital filter
and returns the filter coefficients in length N+1 vector B.
The cut-off frequency Wn must be between 0 < Wn < 1.0, with 1.0
corresponding to half the sample rate. The filter B is real and
has linear phase, i.e., even symmetric coefficients obeying B(k) =
B(N+2-k), k = 1,2,...,N+1.
If Wn is a two-element vector, Wn = [W1 W2], FIR1 returns an
order N bandpass filter with passband W1 < W < W2.
B = FIR1(N,Wn,'high') designs a highpass filter.
B = FIR1(N,Wn,'stop') is a bandstop filter if Wn = [W1 W2].
If Wn is a multi-element vector,
Wn = [W1 W2 W3 W4 W5 ... WN],
FIR1 returns an order N multiband filter with bands
0 < W < W1, W1 < W < W2, ..., WN < W < 1.
B = FIR1(N,Wn,'DC-1') makes the first band a passband.
B = FIR1(N,Wn,'DC-0') makes the first band a stopband.
For filters with a passband near Fs/2, e.g., highpass
and bandstop filters, N must be even.
By default FIR1 uses a Hamming window. Other available windows,
including Boxcar, Hanning, Bartlett, Blackman, Kaiser and Chebwin
can be specified with an optional trailing argument. For example,
B = FIR1(N,Wn,kaiser(N+1,4)) uses a Kaiser window with beta=4.
B = FIR1(N,Wn,'high',chebwin(N+1,R)) uses a Chebyshev window.
By default, the filter is scaled so the center of the first pass band
has magnitude exactly one after windowing. Use a trailing 'noscale'
argument to prevent this scaling, e.g. B = FIR1(N,Wn,'noscale'),
B = FIR1(N,Wn,'high','noscale'), B = FIR1(N,Wn,wind,'noscale').
FREQZ Z-transform digital filter frequency response function
When N is an integer, [H,W] = FREQZ(B,A,N) returns the N-point frequency
vector W in radians and the N-point complex frequency response vector H
of the filter B/A:
jw -jw -jnbw
jw B(e) b(1) + b(2)e + .... + b(nb+1)e
H(e) = ---- = ----------------------------
jw -jw -jnaw
A(e) a(1) + a(2)e + .... + a(na+1)e
given numerator and denominator coefficients in vectors B and A. The
frequency response is evaluated at N points equally spaced around the
upper half of the unit circle. If N isn't specified, it defaults to 512.
FILTER One-dimensional digital filter function
Y = FILTER(B,A,X) filters the data in vector X with the
filter described by vectors A and B to create the filtered
data Y. The filter is a "Direct Form II Transposed"
implementation of the standard difference equation:
a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)
- a(2)*y(n-1) - ... - a(na+1)*y(n-na)
If a(1) is not equal to 1, FILTER normalizes the filter
coefficients by a(1)
Note that in the example M-file of Listing 1, we must set the pole coefficient array, A, equal to 1, which implies that a1 = 1, and all the other coefficients, a2, a3, a4, .... = 0, since an FIR filter has no poles.
_______________________________________________________________________________________
Listing
1. Example MATLAB M-file illustrating FIR filter design and evaluation.
% Finite Impulse Response filter design example
% found in the MATLAB Signal Processing Toolbox
% using the MATLAB FIR1 function (M-file)
Fs=8e3; %Specify Sampling Frequency
Ts=1/Fs; %Sampling period.
Ns=512; %Nr of time samples to be plotted.
t=[0:Ts:Ts*(Ns-1)]; %Make time array that contains Ns elements
%t = [0, Ts, 2Ts, 3Ts,..., (Ns-1)Ts]
f1=500;
f2=1800;
f3=2000;
f4=3200;
x1=sin(2*pi*f1*t); %create sampled sinusoids at different frequencies
x2=sin(2*pi*f2*t);
x3=sin(2*pi*f3*t);
x4=sin(2*pi*f4*t);
x=x1+x2+x3+x4; %Calculate samples for a 4-tone input signal
grid on;
N=16; %FIR1 requires filter order (N) to be EVEN
%when gain = 1 at Fs/2.
W=[0.4 0.6]; %Specify Bandstop filter with stop band between
%0.4*(Fs/2) and 0.6*(Fs/2)
B=FIR1(N,W,'DC-1'); %Design FIR Filter using default (Hamming window.
B %Leaving off semi-colon causes contents of
%B (the FIR coefficients) to be displayed.
A=1; %FIR filters have no poles, only zeros.
freqz(B,A); %Plot frequency response - both amp and phase response.
pause; %User must hit any key on PC keyboard to go on.
figure; %Create a new figure window, so previous one isn't lost.
subplot(2,1,1); %Two subplots will go on this figure window.
Npts=200;
plot(t(1:Npts),x(1:Npts)) %Plot first Npts of this 4-tone input signal
title('Time Plots of Input and Output');
xlabel('time (s)');
ylabel('Input Sig');
%Now apply this filter to our 4-tone test sequence
y = filter(B,A,x);
subplot(2,1,2); %Now go to bottom subplot.
plot(t(1:Npts),y(1:Npts)); %Plot first Npts of filtered signal.
xlabel('time (s)');
ylabel('Filtered Sig');
pause;
figure; %Create a new figure window, so previous one isn't lost.
subplot(2,1,1);
xfftmag=(abs(fft(x,Ns))); %Compute spectrum of input signal.
xfftmagh=xfftmag(1:length(xfftmag)/2);
%Plot only the first half of FFT, since second half is mirror imag
%the first half represents the useful range of frequencies from
%0 to Fs/2, the Nyquist sampling limit.
f=[1:1:length(xfftmagh)]*Fs/Ns; %Make freq array that varies from
%0 Hz to Fs/2 Hz.
plot(f,xfftmagh); %Plot frequency spectrum of input signal
title('Input and Output Spectra');
xlabel('freq (Hz)');
ylabel('Input Spectrum');
subplot(2,1,2);
yfftmag=(abs(fft(y,Ns)));
yfftmagh=yfftmag(1:length(yfftmag)/2);
%Plot only the first half of FFT, since second half is mirror image
%the first half represents the useful range of frequencies from
%0 to Fs/2, the Nyquist sampling limit.
plot(f,yfftmagh); %Plot frequency spectrum of input signal
xlabel('freq (Hz)');
ylabel('Filt Sig Spectrum');
A. Download the M-file of Listing 1, called lab2.m, from the
afs.rose-hulman.edu\class\ee\hoover\ece581\lab2\lab2.mAFS network “class directory”. Start MATLAB and execute this M-file. Verify that this function does what you expect.
B. Now modify this M-file to obtain the sixteen FIR filter coefficients that correspond to a
1. 16th-order band-pass FIR filter with a pass-band between 800 Hz and 2.4 kHz, and a sampling frequency of 8 kHz.
2. 16th-order high-pass FIR filter with unity-gain passband above 2.0 kHz, and a sampling frequency of 8 kHz.
Include a printout of each of these sets of FIR coefficients as Attachment A in your lab report. II. Real-Time, Floating Point FIR Digital Filter Implementation Study the digital filtering program shown in Listing 2. This FIR filtering program was written to be easily understood, and is therefore not very efficient. For example, the “for loop” in the ISR that updates the input sample storage array, x[ ], by shifting the newly converted sample into x[0], what was in x[0] goes into x[1], what was in x[1] goes into x[2], etc., is certainly not efficient, though it makes the convolution calculation very straightforward. The use of a circular input buffer to hold the last N input samples leads to much more efficient code, though the subscript variable management becomes trickier. In this case, a reference pointer (that points to the most recent input sample in the buffer) is simply rotated around the circular buffer, allowing the previous inputs in the buffer remain in their original position within the buffer.Listing 2. C-Language FIR Real-time Digital Filtering Program
_______________________________________________________________________
/* Floating Point FIR Digital Filter Implementation */
/* Digital Signal Processing Laboratory */
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
float coeff[20],x[20]; /* coeff array holds FIR filter coeffs */
int Norder; /* x array holds past input samples */
void hookint(void);
interrupt void McBSPRcvISR(void);
int main()
{
Mcbsp_dev dev;
Mcbsp_config mcbspConfig;
int sampleRate,Actual_Sampling_Rate;
/* Band-Stop Filter Coefficients Cut From MATLAB
B =
Columns 1 through 7
-0.0038 0.0000 0.0218 0.0000 -0.0821 0.0000 0.1625
Columns 8 through 14
0.0000 0.8031 0.0000 0.1625 0.0000 -0.0821 0.0000
Columns 15 through 17
0.0218 0.0000 -0.0038
*/
coeff[0] = -0.0038; /* Here are the 17 coefficients corresponding */
coeff[1] = 0.0000 ; /* to a 16th-order FIR Band Stop FIR filter */
coeff[2] = 0.0218; /* FIR filter obtained using MATLAB FIR1 */
coeff[3] = 0.0000; /* Fs = 8 kHz, stop band between 1.6 kHz and */
coeff[4] = -0.0821; /* 2.4 kHz. */
coeff[5] = 0.0000;
coeff[6] = 0.1625;
coeff[7] = 0.0000;
coeff[8] = 0.8031;
coeff[9] = 0.0000;
coeff[10] = 0.1625;
coeff[11] = 0.0000;
coeff[12] =-0.0821;
coeff[13] = 0.0000;
coeff[14] = 0.0218;
coeff[15] = 0.0000;
coeff[16] =-0.0038;
/******************************************************/
/* Initialize EVM */
/******************************************************/
printf("Initializing EVM board\n");
evm_init();
printf("Done initializing EVM\n");
/******************************************************/
/* Open MCBSP */
/******************************************************/
mcbsp_drv_init(); /* initialize McBSP driver, allocates memory for
the device handles */
dev=mcbsp_open(0); /* dev is the handle to control the McBSP */
if(dev==NULL)
{
printf("Error Opening MCBSP 0\n");
return(ERROR);
}
/******************************************************/
/* Configure McBSP */
/******************************************************/
memset(&mcbspConfig,0,sizeof(mcbspConfig));
mcbspConfig.loopback =FALSE;
mcbspConfig.tx.update =TRUE;
mcbspConfig.tx.clock_mode =CLK_MODE_EXT;
mcbspConfig.tx.frame_length1 =0;
mcbspConfig.tx.word_length1 =WORD_LENGTH_32;
mcbspConfig.rx.update =TRUE;
mcbspConfig.rx.clock_mode =CLK_MODE_EXT;
mcbspConfig.rx.frame_length1 =0;
mcbspConfig.rx.word_length1 =WORD_LENGTH_32;
mcbsp_config(dev,&mcbspConfig); /* configuration adjustments */
MCBSP_ENABLE(0,MCBSP_BOTH); /* McBSP is activated */
/******************************************************/
/* Configure CODEC */
/******************************************************/
codec_init();
/* A/D 0.0dB (min) gain, turn OFF 20dB mic gain, sel(L/R)LINE input */
/* Set codec for stereo mode of operation */
codec_adc_control(LEFT,MIN_ADC_INPUT_GAIN,FALSE,LINE_SEL);
codec_adc_control(RIGHT,MIN_ADC_INPUT_GAIN,FALSE,LINE_SEL);
/* MUTE(L/R)LINE input to “DSP Bypass Path” by setting mute switch to TRUE */
codec_line_in_control(LEFT,MIN_AUX_LINE_GAIN,TRUE);
codec_line_in_control(RIGHT,MIN_AUX_LINE_GAIN,TRUE);
/* D/A 0.0dB attenuation, do not mute DAC outputs */
codec_dac_control(LEFT,0.0,FALSE);
codec_dac_control(RIGHT,0.0,FALSE);
sampleRate=8000;
Actual_Sampling_Rate = codec_change_sample_rate(sampleRate,TRUE);
/* set to the closest allowed rate */
printf("The actual sampling rate is = %d\n",Actual_Sampling_Rate);
hookint();
codec_interrupt_enable();
/* codec generates interrupts when data are received in the DRR */
/******************************************************/
/* Main Loop, wait for Interrupt */
/******************************************************/
Norder = 16;
while(1)
{
}
}
void hookint()
{
/* an interrupt is assigned to DRR event of the serial port
then, the interrupt will branch to ISP */
intr_init(); /* initialize ISTP with the address of vec_table.
Placing the base address of the vector table in ISTP */
intr_map(CPU_INT15,ISN_RINT0); /* map CPU_INT15 to DRR interrupt */
intr_hook(McBSPRcvISR,CPU_INT15); /* connect ISR to the CPU_INT15 */
INTR_ENABLE(15);
INTR_GLOBAL_ENABLE();
return;
}
interrupt void McBSPRcvISR(void)
{
/* ISR for the DRR interrupt */
/*
This routine convolves the present and the Norder previous input
samples with the “Norder+1” FIR coefficients h(0) through h(Norder):
y(n) = h(0)*x(n) + h(1)*x(n-1) + ... + h(Norder)*x(n-Norder)
*/
int intsamp,i;
float floatsamp,sum;
intsamp=MCBSP_READ(0); /* read from CODEC’s data receive register (DRR) */
intsamp = (intsamp >>16); /* Shift right channel data down to bottom 16 bits */
if(intsamp & 0x8000)
intsamp = intsamp 0xffff0000; /* Sign extend right channel data */
floatsamp = (float) intsamp; /* Convert right channel to floating point */
for(i=Norder; i >= 0;i--) /* Update past input sample array, where
x[0] holds present sample, x[1] holds
sample from 1 sample period ago, x[N]
holds sample from N sampling periods ago*/
{ /* This time-consuming loop could be */
x[i]=x[i-1]; /* eliminated if circular buffering were*/
} /* to be employed. */
x[0] = floatsamp;
sum = 0;
for(i=0;i<=Norder;i++) /* Perform FIR filtering (convolution) */
{
sum=sum+x[i]*coeff[i];
}
intsamp = (int) sum; /* Convert result back to integer form */
intsamp = intsamp << 16; /* Shift result back into Right Channel position
while zeroing the Left Channel position */
MCBSP_WRITE(0,intsamp); /* Send to CODEC (right channel only) */
}
This real-time FIR filtering program implements the digital band-stop filter that corresponds to the MATLAB M-file of Listing 1. Recall that this was a stop-band filter with a stop band between 1.6 kHz and 2.4 kHz at a sampling frequency of 8 kHz.
A. Download this program from the previously cited AFS network class directory. It is named “lab2a.c”. Place this file in a subdirectory, along with the other relevant linker command, assembly, and library files, and follow the steps outlined in Lab 1 to build the project. To make sure that the C67x C compiler produces well-optimized code, make sure that compiler optimization is enabled compiler optimization in Code Composer Studio 2 by clicking on Project – Build Options – Basic. Make sure that “Speed Most Critical” option is selected in the “Opt Speed Vs. Size” box. Also, in the “Opt Level” box, choose “Register –o0”. The “Generate Debug Info” box should normally have “Full Symbolic Debug” selected, when you are in the process of debugging your code, though I have found that program execution is substantially speeded up by changing this to “No Debug” after you code is debugged, allowing you to go to a higher sampling rate.Run the program. Verify, using a function generator and an oscilloscope, that this filter attenuates audio signals roughly between 1.6 kHz and 2.4 kHz, as expected. Replace the oscilloscope with a loudspeaker. Note that the attenuation in the stop band is not as apparent with the loudspeaker, since your hearing responds logarithmically rather than linearly.B. Now insert the seventeen filter coefficients you obtained in Part II (B) above for the bandpass filter. Run the program. Listen to the output of this filter when you speak into the microphone. Your voice should have that characteristic band-limited “telephone sound”!
C. Use the Microsoft EXCEL spreadsheet to plot the observed filter magnitude response of this band-pass filter in decibels versus frequency, where AvdB = 20log(Voutpeak/Vinpeak) Vary the frequency using a function generator over a range of 100 Hz – 3.7 kHz. Be sure that the amplitude of the function generator is not turned up too high, since we do not want the filter output to become distorted at any frequency within the specified range. Use an oscilloscope to measure the gain at (at least) eleven evenly-spaced frequencies of 100 Hz, 500 Hz, 900 Hz, 1.3 kHz, 1.7 kHz, 2.1 kHz, 2.5 kHz, 2.9 kHz, 3.3 kHz, 3.7 kHz. You may want to take additional measurements at frequencies where the response of the filter is changing dramatically. You may have trouble estimating the amplitude of the output sinusoid at the higher frequencies, since there are fewer samples per cycle, although the CODEC data sheet claims to provide the proper reconstruction low-pass filtering, with a break frequency of Fs/2. Compare the observed frequency response with that predicted by MATLAB for this band-pass filter. Include both of these plots in your report memo as Attachment B. The experimentally measured frequency response plot of gain (in dB) vs. frequency should have the same shape as the MATLAB-predicted plot, however, it may differ by an additive constant, since the gain of the mixer box is not known (unless the mixer box is bypassed.)
D. Now change the coefficients to those of Part II (B) above (for the high-pass filter). (Isn’t it amazing that changing the coefficient values can so thoroughly change the “personality” of the filter!) Then repeat Parts C and D. Include both the observed and the MATLAB predicted frequency response curves for this high-pass filter in your report memo as Attachment C.
III. Fixed-Point FIR Filter Implementation
Modify the floating point FIR band-pass filtering program that you obtained in Part I (B) so it can be run on the (cheaper and faster) C62x DSP chip, which executes the (fixed-point) subset of the C67 instruction set. The C62x DSP has the same instructions and architecture (even the same pinout!) as the C67x DSP chip. Now the “rules of the game” have changed: you may no longer use any floating point variables or operations in your C program, since floating point operations cannot be efficiently executed by the C62.Obviously, you will have to translate your filter coefficients into integer form. In order to take advantage of the full dynamic range of the C62’s 32-bit integer representation, you should multiply these coefficients by 0x7FFF. This will convert (scale up) these original floating point coefficients into signed 16-bit integer quantities (of type “short int”), since the FIR coefficients of a unity gain filter range between (-1.0, 1.0). These (short int) coefficient values should be entered into your new program, and declared to be of the “short int” data type. In your filtering program, you should multiply by the signed 16-bit input sample in the upper half of the 32-bit integer coming in from the CODEC, corresponding to the right CODEC channel. This may be done without downshifting by 16 bits, and then sign extending, as was done in the floating point version of the program (See Listing 2). Instead, you should use a “C compiler instrinsic function” which forces the use of a specific C62x/67x machine instruction that has no direct counterpart in the C language. The intrinsic C function that is “just what the doctor ordered” in this particular case is “_mpyhl” (multiply high-low), and it may be invoked by a line of C code similar to this: temp=(_mpyhl(sample,coeff));This line of C code forces the in-line insertion of the “MPYHL” C62x/C67x DSP chip’s “16 X 16” integer multiply instruction. This instruction multiplies the 16 MSBs of the first (32-bit) integer argument (sample) by the 16 LSBs of second (32-bit) integer argument (coeff), and returns a 32-bit signed integer product to the integer variable “temp”. The upper 16 bits of this 32-bit variable “temp” holds the16-bit signed output sample. Note that taking this upper 16-bit portion of the 32-bit product in “temp” as the result, corresponds to dividing the 32-bit result by 0x10000, since it lies in the upper 16 bits of the 32-bit integer. However, we want to divide by 0x8000, which almost corrects for the multiplication of the coefficients by 0x7FFF that we had to perform earlier, in order to scale the coefficients for use in this fixed-point (integer) version of the FIR filtering program. Therefore, one final left shift must be applied to our result to make the upper 16-bits of “temp” to be equivalent to the result divided by 0x8000. temp = temp << 1;The 16-bit result is now in the proper position (upper-most 16 bits) to go out to the right channel of the CODEC. However, before sending the result (temp) out to the CODEC, the bottom 16 bits of temp should be masked out, to ensure that no residual noise will go out the CODEC on the left channel.
Demonstrate that your fixed-point FIR band-pass filter program passes frequencies with least attenuation in the range 800 Hz – 2.4 kHz . Obtain the instructor’s validating signature on the program listing. Include the listing of your modified “fixed-point” FIR filter program listing as Attachment D.IV. More Efficient FIR Filter ImplementationWrite, and then demonstrate, the proper operation of a more efficient FIR floating-point band-pass filtering program. As suggested earlier, a circular buffer should be used to replace the input sample storage array x[ ]. This modification would eliminate the need for the time-wasting “for loop” that updates x[ ] each time the interrupt service routine. (Of course, the most efficient FIR implementations would require use of C6x assembly-language coding.) Try running this more efficient filtering program at higher clock rates. Can you get this program to filter properly at 48 kHz, which is the highest sampling frequency supported by the CODEC?
Include the listing of your more efficient filter implementation as Attachment E, and on this listing, indicate the highest sampling rate you were able to attain.
By the way, please note that the band-pass filter’s break frequencies will be scaled up with increasing sampling frequency. Recall that the bandpass filter was designed to break (provide unity gain) between the “normalized frequencies” of 0.1 and 0.3. Thus the filter passband was found to lie between the break frequencies of 0.1*8 kHz = 0.8 kHz = 800 Hz and 0.3*8 kHz = 2.4 kHz, when fs was set to 8 kHz. However, if fs is changed, the break frequencies will lie between 0.1*fs and 0.3*fs.
A simple way to tell if the filter is operating properly (and not missing one or more sampling interrupts) is to gradually increase the frequency of the function generator from 0 Hz up to fs while listening in the loudspeaker. Verify that no aliasing can be hear over this range. Aliasing is heard when the audible frequency starts to go back down, even while the input frequency is being increased. The reason you would expect to hear aliasing is that the CODEC sets its switched-capacitor LPF to ½ of the current sampling rate, but if the interrupt processing is taking too long and every other sample is missed because the previous sample is still being processed, the effective sampling rate is only ½ of what you thought it was, and thus the anti-aliasing filter’s break frequency is set to twice as high as it should be, and so aliasing will be heard. (I will demonstrate this effect in the laboratory!)
.
