EE265 Lab 4: FFT and Spectrum Analysis

Winter 2008-2009
Instructor: Teresa Meng

I. Introduction:

In this Lab you will learn how to use an efficient FFT routine for useful DSP applications. For example, you will use the FFT routine to perform filtering (linear convolution in the time domain) using the techniques of overlap-and-add and overlap-and-save. You will also be able to experiment with some FFT properties; e.g. a circular delay in the time domain means a linear phase in the frequency domain, or the duality of symmetry and real and imaginary components between the two domains. It is also essential that you understand the difference between the frequency response of a finite-duration signal and its DFT.

II. Readings:

Discrete-Time Signal Processing, Second Edition (text)

Section 9.3 - Decimation-in-time FFT algorithms.

Section 9.4 - Decimation-in-frequency FFT algorithms.

Section 8.7 - Linear Convolution Using the Discrete Fourier Transform. Focus on the discussion on overlap-and-save and overlap-and-add methods.

C55x DSP Library Programmer's Guide

Refer to the DSP Library Reference Guide for descriptions of available library functions. (On-line Programming Reference->C55x DSP Library)

Sample code that may be useful is also available in the folder C:\CCStudio_v3.1\c5500\dsplib\examples.

III. Lab exercises:

In the first part of the Lab, you will implement frequency domain filtering in MATLAB. You will then go through a DSP implementation of a small FFT to become familiar with using the DSP library functions and to get a basic understanding of their implementation and performance. After you have done this, you will then go through frequency domain filtering on the DSP, using your knowledge of DSP programming to optimize the filtering algorithm. You will then compare the performance of the frequency-domain filter to the time-domain filter. And finally, you will experiment with some properties of the FFT.

A. Filtering in the Frequency Domain (MATLAB Implementation)

To start, you will first build the framework for filtering in the frequency domain through simulation in MATLAB. This will help you test your algorithm as well as specify the necessary characteristics of the filtering process.

Use the directory lab4/fftfilt/matlab/ for your MATLAB implementation.
Download h_bp40.mat and h_bp80.mat for filter coefficients.
Algorithm overview:
The algorithm for filtering in the frequency domain can be broken down into five parts:

Record the signal and perform the necessary padding for FFT.
Perform FFT.
Multiply the frequency-domain filter coefficients with the FFT of the data signal.
Perform IFFT to obtained the time-domain filtered result.
Post- process the time-domain filtered result to extract useful data. This is your filtered output.

Design Considerations:
There are a number of things to consider when coding this algorithm:

The size of the FFT, the size of the signal buffer, and the size of the filter coefficients are somewhat related and constrained by the following:

FFT sizes are limited to powers of 2 on the DSP.
Multiplication in the frequency domain is equivalent to a circular convolution in the time domain (so steps must be taken to make the result the same as linear convolution).
If the data signal is of length L and the filter is of length P, then the size of the filter circularly convolved with the data signal is L + P - 1.

In the MATLAB code, you should parameterize the size of the data set, the size of the filter, and the FFT size so that you can change any one of these variables without having to go through major changes in the MATLAB code.
It may be useful to structure your MATLAB code to behave like a DSP implementation (for example, only getting one new frame of inputs each time, etc.)
Remember that the time-domain signal is always purely real.
You need to assume that the frequency domain filter coefficients are complex. You must perform complex multiplication of the frequency-domain filter coefficients and the FFT of the data signal for the filtering to work. MATLAB automatically performs complex multiplication when the signals are complex, but you will have to do this explicitly when you implement the algorithm for the DSP.
After you have completed IFFT, you will need to post-process the result to extract useful filtered data. There are two methods that you can employ to do this (see Readings for more information):

Overlap-and-add - where the overlapping section from the circular convolution is preserved. This is then summed with the corresponding overlapping section of the next filtered segment - See figure 8.23 in Discrete-Time Signal Processing, Oppenheim & Schafer, Second Edition.
Overlap-and-save - where the overlapping section is discarded and only the useful (un-corrupted) filtered outputs are retained - see figure 8.24 in O & S.

You are to choose one of the two methods that you think is the most efficient for the DSP. You will need to implement only one method in MATLAB and on the DSP.

MATLAB Warning: If you need to take the transpose of a vector/matrix of complex numbers in MATLAB, make sure that you use transpose operator ( .' ) instead of the conjugate transpose operator ( ' ).

In order to take the simple transpose of a matrix, you will need to use the " .' " operator (example: y_t = y.').
The transpose operator commonly used is " ' " (example: y_ct = y'), which gives you the complex-conjugate transpose of y.
Something to think about: if you take the the complex-conjugate transpose of a signal's frequency domain coefficients, what does this correspond to in the time domain?

In MATLAB, you are to do the following experiments:

Experiments - implement frequency-domain filtering with the following fiter sizes and FFT sizes:

Filter size = 40 taps, FFT size = 128 points
Filter size = 80 taps, FFT size = 128 points
Filter size = 40 taps, FFT size = 256 points
Filter size = 80 taps, FFT size = 256 points

Filter coefficients:

Sample bandpass filter coefficients have been generated and provided to you in the lab4/fftfilt/matlab/ directory. The file h_bp40.mat and the file h_bp80.mat contain the 40-point bandpass filter coefficient array and the 80-point bandpass filter coefficient array, respectively. You can load these into MATLAB using the load command. For example, load h_bp40 will load the variable, h_bp, contained in h_bp40.mat. These coefficients have unity bandpass spectrum. To use these coefficients in a DSP code, you need to multiply these coefficients by 2^15.
The sample bandpass filters are generated using Kaiser windows and have the following characteristics:

40-point:

passband = 1000 - 2500 Hz
passband ripple = 0.01%
stopband ripple = 0.1%
Transition width = 450 Hz assuming Fs = 8000 Hz

80-point:

passband = 1000 - 2500 Hz
passband ripple = 0.002%
stopband ripple = 0.003%
Width = 322 Hz assuming Fs = 8000 Hz

Verifying your design:

if your frequency-domain filtering code is working correctly, then the output it produces should be equivalent to the result of performing linear convolution of the input data and the filter impulse response. Use Matlab's filter(h,1,x) to compare, where h is your filter and x is your data input. This will truncate the last P-1 values of the full linear convolution, which is what you want.

Estimating Performance:

You are looking at two tradeoffs, performance in terms of number of operations per data sample and quality of filtered result in terms of number of filter taps. How does doubling the number of taps affect the performance in terms of CPU time? How does doubling the FFT size affect the performance in terms of CPU time?
Use the tic and toc commands to calculate the relative cpu_time of filtering 100,000 samples with your overlap-add or overlap-save implementation (type "help tic" for info in MATLAB). This will give you an idea of the relative amount of cpu time needed to perform the filtering operations with different parameter settings (filter length, FFT size, etc.).

For best results, put tic right before your main loop and toc right after your main loop.
Note: the actual number of samples is not critical. The idea is to get the tic/toc times to be on the order of seconds so that we can observe changes.
If your code is taking more than 10 seconds to process 100,000 samples, feel free to reduce the number of samples until the tic/toc times reported are in the 1-10 second range.
(In the past, MATLAB provided a FLOPS command that was more precise for this purpose, but it is now obsolete, hence the need for this cruder measurement with a large data size.)

MATLAB Speed-up: Depending on how you implement your MATLAB code, the majority of your algorithm's running time may be used for memory allocation. For example, the common practice of concatenating two vectors ( yy = [yy next_frame] ) forces MATLAB to perform memory allocation each time the size of the vector grows. This is OK, but it makes the running time of your MATLAB script dominated by this memory allocation instead of the FFTs involved. This can be prevented with a little bit of extra work:

Since we know the final size of all of our vectors, we can pre-allocate the vectors to be their final size:

yy = zeros(1, final_size);

The new data can be placed into this vector by using indexing:

yy(1:frame_size) = next_frame;

In general, this can be done in a for-loop with something like this, where ii is a loop index starting at 1:

yy( (ii-1)*frame_size + 1: ii*frame_size) = next_frame;

This can be very important if we are trying to estimate the effect that the FFT size has on the performance of the algorithm. Otherwise, our running time is dominated simply by the number of vector concatenations that must be done and the effect of the FFT calculation will not be noticed directly.

Your numbers probably won't match the expected results. For example, you may get that P=40, N=256 is twice as fast as P=40, N=128, and the speedup should theoretically be something else. Try calling fftw('wisdom',[]); after every fft/ifft call. Matlab likes to somehow cache its fft/ifft calculations, so this command clears the fftw database (basically its fft cache) so that subsequent fft calculations take as long as they're supposed to.

B. FFT and Intro to the DSP Library

In this section, you will implement a 64-point FFT using the C55x DSP Library. The DSP Library contains many useful functions for DSP algorithms. Examples of such functions are FFT routines and trigonometric functions. TI has made this library available so that programmers do not need to implement these functions themselves and can use these functions as building blocks for implementing more complex DSP algorithms. The following exercise will help you become familiar with using the DSP Library in a project. For information on the DSP Library, refer to the DSP Library Reference Guide. (On-line Programming Reference->C55x DSP Library)

Coding Requirements:
Implement a 64-point FFT followed by a 64-point Inverse FFT on the EVM using the DSP library.

Download and decompress fft64.zip to your account under labs/lab4/fft64/.
Open the project fft.pjt and resolve any compiler/linker errors.

We use two libraries; runtime support library, which is standard C library, and dsp library.
NOTE: CCS may complain that it cannot find the library file 55xdspx.lib. This is because the project file was created on a machine whose CCS installation was on a particular drive letter (C:\). If you are on a machine for which CCS was intalled on a different drive letter, then you may have to change the drive letter for the location of the actual DSP library .lib file. The .lib file should be located at either C:\CCStudio_v3.1\c5500\dsplib\55xdspx.lib or D:\CCStudio_v3.1\c5500\dsplib\55xdspx.lib, depending on the machine that you are using.
55xdsp.lib is for small memory model project, whereas 55xdspx.lib is for large memory model project. Note that fft.pjt is configured as large memory model project.
Similarly, CCS may complain that it cannot find a header (.h) file. In this case, you will have to change the drive letter for the location of the DSP library "include" directory. Under the Compiler tab under Build Options, (Project->Build Options), change the drive letter of the include path from C:\CCStudio_v3.1\c5500\dsplib\include to D:\CCStudio_v3.1\c5500\include or vice versa to ensure that the compiler knows where to find the proper header (.h) files to include.
Runtime support library and header files can be found in C:\CCStudio_v3.1\c5500\cgtools\lib and C:\CCStudio_v3.1\c5500\cgtools\include.
Similarly, rts55.lib is for small memory model, and rts55x.lib is for large memory model.

Familiarize yourself with this project template. You will be using it to integrate DSP library with your real-time audio code later on so that you can hear the output of your DSP code. ( Refer to the DSP Library Reference (On-line Programming Reference->C55x DSP Library) for information on available functions/macros.) There are several noteworthy points about this template:

It serves as an example of running C programs without the use of DSP/BIOS.
Note that the library file 55xdspx.lib and rts55x.lib have been included in the Libraries folder of the project. This allows the linker to access the DSP Library and runtime support functions.
Also note that the header files for the DSP Library and runtime support functions have been included in the Compiler tab under Build Options (Project->Build Options). The include path is specified with the "-i" option (as either C:\CCStudio_v3.1\c5500\dsplib\include or D:\CCStudio_v3.1\c5500\dsplib\include, and C:\CCStudio_v3.1\c5500\cgtools\include or D:\CCStudio_v3.1\c5500\cgtools\include).

Edit main.c by adding the appropriate DSP library functions with the appropriate arguments. All the DSP library setup (i.e. memory sections, variables, etc.) are already implemented in this template.

You should either implement a cfft followed by a cifft, or implement an rfft followed by rifft.
See below for a description of how to use the DSP Library properly.

Profile the 64-point FFT routine. Compare the Profiled result from the DSP code and the documentation on the fft function in the DSP library. What are the discrepencies between the actual EVM code and the DSP library performance table? Would you have easily been able to write a 64-point FFT code that out performs the DSP library routine?

Proper use of the DSP library:

As specified in the DSP library user manual, you need to include the .twiddle section in your linker command file. In addition, you must allocate this section onto the internal memory of DSP with the alignment of 2048. Use .twiddle: {} > MEMORY_SECTION PAGE 0, align (2048) in linker command file.
Most FFT functions (cfft, rfft) requires the elements of their input arrays to be in bit-reversed order. cbrev is a complex bit reversal routine that will perform this re-ordering for complex data (with a real word and an imaginary word). cbrev function expects that the input array of this function to be aligned with even word boundary.
Pay attention to the scale factor presented in the DSP library. If scale == 1, scale factor = nx, else scale factor == 1. Scaling determines whether or not the result will be divided by N (the FFT size).
Specifically _cfft_NOSCALE and _cifft_NOSCALE function don't take advantages of DSP's automatic saturation feature. Test cfft / cifft pair with full range sine wave. If you want to enable saturation when using these functions, use asm("\tBSET SATD"); prior to calling a function. Be sure to reset the processor status after that function by using asm("\tBCLR SATD"); Note that these functions are called when using cfft, cifft, rfft, rifft functions with NOSCALE option.
Several known issues with the rfft/rifft function (DO NOT IGNORE THESE!)

1. rfft and rifft are, in fact, macros. See dsplib.h file for the definitions; RFFT is composed of cfft / cbrev / unpack, and RIFFT comprises unpacki / cifft / cbrev.

2. More RIFFT Issues: The rifft function also suffers from an undocumented shift of two bits (factor of 4) such that taking the rfft followed by an rifft will result in an attenuation by 4 (approximately, not exactly 4). This affects the quality of waveform slightly, and it seems that the accuracy at the edge of block is degraded quite a bit. You may or may not need to compensate for this attenuation. It depends on the required precision of your application.

3. RFFT: The documentation for the rfft has some slight typos in the description of the packing scheme used for the frequency domain output. First of all, the second word of the output array is the Nyquist term, which should be purely real for real input values. The documentation mistakenly refers to this as an imaginary term (y(nx/2)Im instead of y(nx/2)Re). Secondly, the last two words of the output array will be the real and imaginary terms for the frequency bin before the Nyquist frequency, so the corresponding indices should be listed as (nx/2-1) instead of (nx/2).

4. Scaled function (cfft_SCALE, cifft_SCALE) introduces some noise, and if the fft spectrum, which is generated by cfft_SCALE, is inverse transformed, the the time-domain response would have quite a larger error on the first two samples than on the other samples. This will be a critical problem with overlap-and-addition because this method doesn't throw away the first two samples.

There may be inconsistencies in the DSP Library documentation. If you run into any strange behavior with the DSP Library you are strongly encouraged to isolate and test each library call in your code. MATLAB can be very useful for this; writing a script to inspect the actual input and output data from DSP routines can save you a lot of debugging time.
In order to suppress noise level sufficiently small in 16bit audio processing, it is strongly recommended to use 32bit version of functions such as cfft32, cifft32, rfft32, and rifft32. Put your 16bit input data on the upper 16bit of 32bit input data to transform, then do the internal processing, whatever it may be, maintaining data width to 32bit accuracy. After you do the inverse transform, the upper 16bit data of the final 32bit data can be extracted.

FFT Scaling:

One of the arguments to the FFT/IFFT routines is scale, which specifies whether or not the DFT summation will be divided by N, the FFT size.
Traditionally, the (forward) DFT does not scale by 1/N, while the (inverse) IDFT does scale by 1/N. (Take a look at the definitions).
However, on the DSP, we typically scale by 1/N in the forward direction to avoid overflow. Here is an example where overflow will occur:

If you have an input array of all 1.0's (fractional), the DC term of an N-point FFT should then be N. Since this cannot be properly stored as a 16-bit fraction (Q15 format) the result will overflow into the guard bits of the accumulator, and the 16-bit result (taken from AH or BH) will be wrong.
If we scale by 1/N for this same input, our output is guaranteed to be no greater than 1.0, so it will not be affected by overflow.

Note that if we scale by 1/N in the forward direction, we do not need to scale by N in the reverse direction, so the two operations will still be inverses of each other.
Scaling by 1/N in the forward direction can cause interesting problems, since the forward FFT is not normally scaled by 1/N. For example:

Say we take the forward FFT of two signals, x and y, and scale by 1/N while taking the FFT.
If we now multiply the frequency domain signals together, we get a result that is scaled by 1/(N*N).
This may cause problems, depending on what you are doing, so you will want to keep this in mind.

MATLAB Verification:
You will be expected to write a script to verify your code. This will prove that your routines work as expected and it will also help you fully understand the DSP Library routines that you will be using. In particular, it will force you to look at the scaling done by the routines and the output format of the FFT routines. The work you do here should make your job easier for the second part of this lab, as well as for the final project.
Here are the requirements.

Verify your FFT routine (either cfft32 or rfft32). This should involve several things:

Overwrite the input to the FFT with some data before performing the FFT. Using random data proves that the routine works in the general case, but using a signal such as a sine wave may give you more insight.
Perform the FFT (and other necessary routines such as cbrev32) and verify that your routines produce the proper output for this input.
Use MATLAB to perform the FFT on the same input data.
Compare your DSP result to the result generated by MATLAB. There are a couple of considerations here:

The output of your FFT routine (cfft32 or rfft32) will be in some packed form (such as {real, imaginary, real, imaginary, ...} or something more complicated). You can deal with this in one of two ways:

Either expand the packed data out to an array of complex numbers (the form MATLAB uses), or
Take the MATLAB result and pack it into the same form as the DSP's FFT output

You will have to account for scaling differences due to the use of "fractional" integers, the scaling parameters passed to the FFT routine, and also any scaling inherent in the FFT routines themselves (rifft32). By going through this, you will be forced to understand how the scaling works before moving on to debug larger problems.

Verify your Inverse FFT (IFFT) routine (either cifft32 or rifft32). This should involve similar issues to the FFT verification described above. In addition:

You will now be providing a complex, frequency-domain input to the IFFT routine. This data is in some packed form (the same as the output from the cifft32 or rifft32).
Use the output from your FFT routine as the input to the IFFT routine. However, note the following:

This does not effectively isolate the routines (the IFFT output will depend on the FFT routine, and we would like to verify the IFFT separately).
If your FFT is already verified, then passing its output to the IFFT should be OK and it will allow you to see the result of a signal passing through the FFT and IFFT in succession. This will expose you to any possible scaling issues.
To test the IFFT in isolation, you would have to generate some input data and make sure that it is packed into the expected form. Once again, when verifying the result in MATLAB, you will need a single array of complex numbers (for taking the IFFT), so the array of data that goes into your DSP IFFT will not be the same as the array that MATLAB uses. This can be resolved in a couple of ways:

Either create random 32-bit data to be passed to the DSP and then un-pack this into a single complex array for use by MATLAB, or
Create a single array of random complex numbers for MATLAB and then pack this into the proper form before writing to the DSP

Once again, you will need to be aware of scaling issues and get a full understanding of the scaling caused by IFFT routines. As mentioned above, rifft has a scaling quirk that you may want to be aware of.

Recommendation: Create MATLAB functions

It is recommended that you write some MATLAB functions to convert the packed output from DSP FFT routines to unpacked arrays for use by MATLAB. This will allow you to make meaningful MATLAB plots of the frequency-domain data produced by your FFT routines and also to use this data in MATLAB computations.
It would also be very useful to write a function that does the reverse (converts from a complex array to a packed array).
rfft consideration: If you are writing a routine to expand the output of rfft32 you will have to remember that only half of the spectrum is represented in the array and the other half will have to be recreated as the complex conjugate of the first half.
Writing these functions now will make your MATLAB code cleaner and will allow you to reuse them later on (this lab, and the final project).
Unpacking example: Here is a simple example of converting a "packed" cfft32 output to a complex MATLAB vector:

Assume that we have some input that is packed as (re, im, re, im, ...). We will use the numbers 1 through 8 for this example. We can extract the "real" components using indexing starting at index 1, jumping by 2, and going to the end of x. Similarly, we can extract the "imaginary" components starting at index 2. These real and imaginary vectors can be combined into a complex vector in MATLAB as follows :

x = 1:8; (this is the same as x = [1 2 3 4 5 6 7 8])
real_x = x(1:2:end); (this gives us real_x = [1 3 5 7] )
imag_x = x(2:2:end); (this gives us imag_x = [2 4 6 8] )
unpacked_x = real_x + j * imag_x; (this gives us unpacked_x = [(1 + j*2) (3+j*4) (5+j*6) (7+j*8)] which is what we wanted.)

The rfft output can be unpacked using something slightly more complicated. This may require flipping a vector, which can be done by specifying a decreasing index range or by using the flip left-right or up-down functions (fliplr/flipud):

x_flipped = x(end:-1:1); (this gives us x_flipped = [8 7 6 5 4 3 2 1] )
x_flipped2 = fliplr(x); (this gives us x_flipped2 = [8 7 6 5 4 3 2 1] )

Caution: Overflow

It is now possible to experience overflow problems in your code. Here is an example:

If you have an input array of all 1.0's (fractional) and the DC term of an N-point FFT should then be N. Since this cannot be properly stored as a fraction (Q15 or Q31 format) the result will overflow.

This means that if you have random data that is large (uses up the full range of the 16-bits), you may get overflow after taking the FFT.
You can check for overflow by using MATLAB to determine if the ideal result is greater than 32,767 (after accounting for all other scaling: fractional, and otherwise). This may take some thinking, but it is important to be aware of it.

C. Filtering in the Frequency Domain (DSP Implementation)

In a previous section, you have gone through the frequency-domain filtering algorithm in MATLAB. Here, you will port it to the DSP. You will use the DSP library to perform the FFT portion of this algorithm. Your job is to integrate this algorithm into the real-time processing environment that was used in earlier labs so that you will be able to hear the output of your frequency domain filter.

1. DSP Implementation

The working directory for this section is lab4/fftfilt/asm/.
For this section, there are two things you need to do. Implement frequency-domain filtering algorithm using the DSP library and integrate this algorithm into a real-time code. There are two approaches depending on what code template you start with. Choose one of these two approaches for implementation:

Recommended approach: If you started with a real-time code, then you need to integrate DSP library into it. This will ensure that once you are successful, you are done. However, you may encounter real-time issues that can sidetrack you in implementing the algorithm. If you decide to implement this way, then:

Copy the source codes from a working sinewave generator or audio recorder that you did in earlier labs to lab4/fftfilt/asm/.
Modify the project to include the DSP Library .lib file and the proper include directory for the library function header files (as in the earlier part of this lab).
(DON'T include rts55x.lib. That's not needed for DSP/BIOS application)
Implement the frequency-domain filter using DSP library routines directly into the real-time code.

You can start by testing the algorithm using a code that you know works with DSP library. Once you have verified your algorithm, you can then try to integrate this into a working real-time code. In this case, the real-time debugging would be done during the integration instead of during algorithm porting. If you decide to implement this way, then:

Copy the source codes from lab4/fft64/ (a working code that uses the DSP library) that you have completed in Section III.B. to lab4/fftfilt/asm/dsplib/.
Implement the frequency-domain filter using DSP library routines and make sure it works in this code first.
When you have verified that your frequency-domain filter works, integrate it into a real-time code. Copy the source codes from a working sinewave generator or audio recorder that you did in earlier labs to lab4/fftfilt/asm/.

Specifications:

Filter size = 40 taps. You may choose to use the sample filter coefficients provided to you in lab4/fftfilt/matlab/h_bp40.mat. However, you will need to convert the filter coefficients into the appropriate format.
FFT size = 128 points.
From the MATLAB code, determine the number of new input samples that can be processed for each pass of the algorithm. Remember that this is a function of the filter size and FFT size.
For real-time integration:

Set the sampling rate to 8 KHz in the real-time code.
You are to choose the size of your audio buffer.

You will be graded on the functionality and performance of this code. (Note that this does not suggest hand-optimizing the DSP library functions yourself.)
You will also be required to perform verification with MATLAB, as described below. The MATLAB verification capabilities can be very useful for verifying individual parts of your code (complex multiplication functions, etc.) so it is recommended that you use this during your coding process. (More below...).
Be sure to adjust the buffer size appropriately in audio.h ( #define BUFFSIZE N )
Also, you should filter both channels of audio

Algorithmic Comparison:
Profile the core filtering code and compare the performance to the time domain filter in Lab 3 Section II.B.1. In order to to make a fair comparison, both implementations should be based on the same number of filter taps. You will need to estimate the execution time that your lab 3 implementation would have taken for a 40 tap filter.

Based on the MATLAB exercise in section B.1, can you determine an estimate of the size of the FFT used for the frequency-domain filter such that the frequency-domain filter is more efficient than the time-domain filter?
Does such a size exist? If not, can you explain why?
Here are the factors involved in comparing the two algorithms:

The parameters: Filter length is P, Number of input samples is L, FFT size is N
In terms of algorithmic complexity, linear convolution that takes in L input samples and produces L new output samples is on the order of L*P operations: O(L*P). This means that the number of operations/cycles can be estimated as k1*L*P, for some constant k1.
For the frequency-domain filter, the algorithmic complexity is on the order of N*log2(N) operations O(N*log2(N)). This means that the number of operations can be estimated as k2*N*log2(N), for some constant k2.
In order to directly compare the two methods, we must relate the parameters. This is typically done by fixing some of the parameters to be constants or by making one parameter a function of another one. Here are a couple possibilities:

We could fix P and make N a function of the input size L and the constant P (since N will be related to the filter length and the number of input samples)
We could let P vary as a function of L or N and then reduce the problem down to a single variable. Example: let P = N/2 and L = ? (constrained by N and P). This gives us everything in terms of a single parameter N.

If we relate both filtering methods by way of a single parameter, we still cannot directly compare the two until we know the constants used (k1 and k2). We really only need an estimate of these constants, and we can get these estimates from our knowledge of the DSP. For example, for linear convolution, we know that k1 will be approximately 1 if k1*L*P represents the number of cycles to perform linear convolution (as we saw in Lab 3). The constant k2 can be estimated from the cycle count required for the FFT functions in the DSP library.
Once we have an estimate of the constants involved (and both equations are related by a single parameter), then we can estimate which algorithm will be more efficient as a function of the single parameter. The result will depend heavily on the choices made (fixed P or P=N/2, etc.) so any conclusion should be qualified by stating these assumptions.

2. MATLAB Verification

As in Lab 3, you will be expected to write a script to verify your code. Here are the requirements.

At the very least, you will be required to verify that the output of your frequency domain algorithm is equivalent to the ideal (MATLAB) filtered output, as described in Lab 3. This should have already been done by your Lab 3 script, so this shouldn't require too much work.
Debugging Suggestion:
Although you are only required to verify the output of the entire algorithm, you are strongly encouraged to use the MATLAB verification capabilities to test individual components of your design. A lot of time can be saved by debugging smaller pieces of a design, instead of attempting to debug the entire design.

For example, your frequency-domain filtering algorithm will need to have some form of complex multiplication routine for frequency-domain samples. MATLAB can be used to isolate and test this routine, since it can overwrite the input before performing the operation and can then verify that the result of the operation was correct for the given input.
In this lab and in future labs, the algorithms we implement will typically consist of several different stages that are performed sequentially (example: FFT -> complex multiply -> inverse FFT). Each of these stages has an input and an output (which is usually fed to the input of the next stage). MATLAB verification can be used to debug each of these stages separately, instead of treating the completed system as a black box and trying to verify its output from its input.

3. "Ideal" Filter

In order to better understand filter design, you will now filter the frequency response of the input by multiplying with an ideal filter. This ideal filter has a unity passband and zero stopband (as specified in the frequency domain).

The working directory for this section is lab4/fftfilt/ideal/.
Copy the working frequency-domain filter code from the working directory of the previous section into this working directory.
Plot the magnitude and phase of the 128-point FFT of the filter coefficient stored in h_bp40.mat. Note that the phase is linear but the magnitude has small ripples.
What would happen if rather than dealing with non-flat magnitude response in the frequency domain, you create an ideal filter by simply setting the passband filter frequency coefficients to unity (zero phase) and the stopband filter frequency coefficients to zero?
To see the effect of the above question, replace the filter coefficients in your DSP implementation with the filter coefficients of an ideal filter. Set the passband to be 1 KHz to 2.5 KHz (same as the sample filter coefficients) (assuming an 8kHz sampling rate).
For the purposes of your DSP code, you can specify the frequency-domain coefficients directly, without computing the time-domain sequence. (Note: we will eventually ask you to compute the time-domain impulse response, so it may be easier to compute this first.) If you do specify the frequency domain values for the DSP directly, make sure that the format of these samples is compatible with your code (like the output of a cfft32 or rfft32 call, most likely).
Qualitatively compare the output of this ideal filter to the output of the 40-tap sample filter.
Can you explain why this ideal filter does not work? There are three distinct reasons, and the following should help you see why:

Provide a MATLAB plot of the time-domain impulse response of this filter. Note the length, and comment on how this effects the algorithm.
Also provide a plot of the actual magnitude spectrum of this filter (using zero-padding to get higher resolution in the frequency domain). Observe the behavior of the filter for frequencies in-between the frequency-domain samples that were forced to 0 or 1).
How well does this filter satisfy the requirements of your filtering algorithm?

Last modified: 9:00 am 3/07/2007