DSP Projects: Lab Project -2

FIR Filter Implementation in MATLAB and in C


I. 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\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!)

Sitting for the RHCE

By Ken Barber on April 06, 2004 (8:00:00 AM)

I've taken some pretty tough tests in my life, and passed every one of them on the first try. And since I'm teaching Linux system administration at my local community college, I thought I would be hot stuff when I signed up to take the Red Hat Certified Engineer (RHCE) exam. I don't think that way anymore.

I've always done well on tests. Mensa once offered me a membership in their organization because of my SAT scores. In a former life I was a refrigeration mechanic and earned the coveted CM certification from the Refrigeration Service Engineers Society -- a cert well known in that industry for requiring multiple attempts to pass -- on the first try. I even passed the written exam for a private pilot's license on my first try.

Then I started a new career in IT. I breezed through my MCSE a few years ago, including the notoriously difficult Exchange Server exam. Then I earned a GSEC security certification from the SANS Institute, requiring two open-book exams. I took the first one closed-book, and passed both. When I learned that the Linux Professional Institute was offering its exams for free at Linuxworld last year, I walked in late to one of them on a lark, totally unprepared. I barely squeaked through, but I passed.

RHCE was different. It was the second-toughest exam I've ever taken.

Red Hat vs. other certifications

The day after I earned my MCSE the Dilbert comic strip featured a character in tights and a cape saying, "Step away from that network server! I'm certified!", but he couldn't fix anything when he sat down in front of it. My friends and co-workers thought it was really funny, and made sure I had clippings of it coming in from all over for the next few days.

But the problem of people holding "paper" certifications is not a laughing matter to employers, and the number of "paper CNEs" and "paper MCSEs" has become the stuff of legends. I do not believe there is any such thing as a "paper RHCE."

To my knowledge, there are only two IT industry certifications that require a candidate to set up and repair an actual running system. Red Hat's is one of them; the other is a Cisco exam. There are no multiple-choice questions to answer; you spend the entire session repairing a broken system and then building a new one from scratch. At the end of the day, the things you've been asked to do either work, or they do not -- and you pass or fail on that basis alone.

It's not as easy as it sounds. The failure rate hovers around 40%.

Fortunately, all is not necessarily lost for those unlucky 40% who fail to meet the lofty requirements for an RHCE. For a little more than a year now, Red Hat has been issuing a Red Hat Certified Technician (RHCT) certificate to those who demonstrate competence in the portions of the test that deal with workstation (as opposed to server) administration.

Preparing for the exam

What's on the test? I'm not allowed to tell you, but Red Hat provides a list of everything you need to know. It is safe to assume that you will be tested on every item in that list; if you're weak in two or three of the list's items, don't take the exam until you've done some more preparation.

The choices for "more preparation" are somewhat limited: you can take a class from Red Hat, or you can take a different class from Red Hat! Red Hat does provide a very nice set of online pre-assessment tests (free, but registration required) to help you choose which class is right for you. The exam prep guides being sold by various booksellers can help you prepare for a class, but are useless in preparing for the exam itself. Believe me -- I own two of them.

The main value of the classes is to learn how do the exam's tasks quickly, the way Red Hat wants them done. But I can tell you from experience that even Red Hat's classes are not enough to fully prepare you. If you want to have any hope of passing, you will have to have been installing, repairing, and configuring Linux, in all of the areas mentioned in the list above, often enough and long enough to do most of it without referring to a man page. Yes, you're allowed to use man pages during the exam, but if your system won't even boot when you walk into the room, you had better know your stuff cold.

I took the RH 300 course, the one for people who supposedly already know what they're doing. It runs for the four days prior to test day.

Exam day

On a recent Friday morning eight of us filed into the classroom we'd been calling home for the last week to find our systems re-imaged. After signing non-disclosure agreements, we were given a list of 10 things that weren't working and 2-½ hours to fix them. Five of the items had to be fixed in the first hour.

You must get at least eight of the 10 items fixed to earn an RHCE. If you get only the first five items, in the first hour, you can earn an RHCT.

One of our number didn't make it through the first hour. No one gloated when he shook the proctor's hand and left early. The pervasive feeling in the room was that any of us could have met the same fate. The mood afterward, as we ate our catered lunch, was somber. One of us had already been eliminated. How many more of us wouldn't make it through the three hours to come?

We returned to the test room to find that our hard drives had been erased. We were given a boot CD with instructions to build a server with an unbelievably huge list of requirements. As I looked over the list, my heart sank. "There is easily two days' worth of work here," I said to myself, "and I have to have it all done in three hours?"

The list was divided into RHCT tasks and RHCE tasks. You must score at least 70% in each area to earn an RHCE, and a high score in one won't help in the other.

A few of the items on the list were easy, but many were not. So I set myself to the task, starting with the things that I already knew how to do and plugged along at my usual snail's pace. For the things that I'd never done before taking the class, the man pages and online guides weren't enough help. There is simply not enough time to read them, and doing a task once in a class exercise isn't enough to remember how to do it cold.

I actually got through the entire list of items about five minutes before the time was up. Only one other candidate had finished; from the sound of feverish key-tapping in the rest of the room I guessed that most were still trying desperately to get as much done as possible before the bell. All of the things I had tested worked, but I didn't test everything. And there were enough of those to sink my boat if they didn't actually work.

I decided there was no use trying to test any more stuff because if they didn't work there wouldn't be time to fix them anyway. So I rebooted to make sure the system came back up with the services running that were needed, and then just shut the thing down and leaned back in my chair. "I'm done," I said to no one in particular. "Either I passed or I didn't, but I'm not doing any more to this machine." I shook the proctor's hand and left about two minutes before the bell rang.

The next morning my head still hurt, and the pain went down into my shoulders. My sweetheart tried to massage my neck, but said that my muscles there were so tight that they felt like bones. In the end, I had to take a muscle relaxer to calm down. Still, I was cautiously optimistic. While I knew I hadn't aced the exam, I thought I had a good shot at getting the double 70%.

A few days later, I received my scores:

SECTION I: TROUBLESHOOTING AND SYSTEM MAINTENANCE
Overall Section I score: 100%

SECTION II: INSTALLATION AND CONFIGURATION
RHCT components score: 100.0%
RHCE components score: 67.9%

RHCT Certification: PASS
RHCE Certification: NO PASS

It's a brutal exam.

DSP Lab-3

Laboratory 3:

Low-pass FIR Filter Design

I. Introduction

In the present project, students are required to implement and simulate a low-pass FIR filter using DSP Builder in the Simulink environment. The design has to be downloaded to the FPGA device on the Stratix EP1S25 DSP development board to perform hardware simulation and verification.

II. Theory

An FIR (Finite Impulse Response) filter, oppositely to IIR filters, has a finite response to impulse signals, which is explained because it does not have feedback. This way, FIR filters define a class of filter that has only zeros in the z-transform (the poles are in the origin z=0).

In addition, FIR filters have other characteristics that become these filters very attractive to many applications, such as linearity of phase, stability in the frequency response and constant group delay.

Equation describes an FIR filter of length K:

Where:

x and y represent the input and transformed data, respectively.

a_k is the set of constant coefficients of the filter.

(K-1) is the order of the FIR filter.

An FIR filter can also be characterized by its number of taps (K), which is the order incremented by one. The transfer function A(z) of the FIR filter is expressed as follows:

Given equation, an FIR filter is also called an all-zero filter because the frequency response is only determined by the zeros in the z-transform.

In general, FIR filters are preferred because of its linear phase characteristic and stability. However, IIR filters can be used in applications that require sharp cut-off or narrow band filters and where linear phase is not a requirement. That is because FIR filters require much higher order implementations than IIR filters for a similar performance.

III. Background about DSP Builder and the Altera EP1S25 DSP development board

Capabilities of DSP Builder and Simulink were introduced with an example design in the Tutorial: “A/D and D/A Conversion on Altera Stratix EP1S25 Development Board using Simulink and DSP Builder” [1], which contain a design targeting the 12-bit 125MHz A/D and 14-bit 165 MHz D/A converters onboard. For further information about the Stratix EP1S25 DSP development board, the student is referred to [3] and [4].

Similarly to the previous lab design, the present project will include multirate capabilities by using the PLL blocks in the Altera FPGA. The student is referred to Laboratory 1 [2] for information about working with different sampling rates in the same circuit design.

IV. Requirements

For this laboratory, students are required to accomplish the next requirements:

• Define the floating-point coefficients for a low-pass FIR filter with the next specifications:

Direct-form FIR filter

Cut-off frequency: 0.1178 x F_N

Stop-band frequency: 0.2958 x F_N

Max. attenuation pass-band: 3dB

Min. attenuation stop-band: 97dB

Use fdatool from Matlab to design the filter under the previous specifications and determine the number of coefficients. To open this tool, type ‘fdatool’ in the Matlab command window. Set the Design Method of the FIR filter to Equiripple and the Design Factor to 20. It is important to note that we are using normalized frequency F_N with a range from 0.0 to 1.0. The maximum value (1.0) corresponds to the half of the sampling rate.

• The coefficients of the designed filter using fdatool are in floating-point format. You can plot very easily the frequency response of the filter by using one of the available options in fdatool (Analysis>Magnitude Response). Check if the filter requirements are fulfilled.

• Export the coefficients from fdatool to the Matlab Workspace (File>Export…). By default, coefficients are stored in the variable ‘Num’. Then, convert coefficients to 16-bit fixed-point with the next command in Matlab:

Num1 = round (Num*2^16-1)

• Check if the new fixed-point coefficients fulfill the filter specifications. File > Import Filter from Workspace permits to analyze the fixed-point design. First, select ‘Direct-Form FIR’ from the drop down menu in Filter Structure and then type Num1 in the ‘Numerator’ field, which contains the 16-bit coefficients. In the ‘Units’ field, ‘Normalized (0 to 1)’ has to be select from the drop down menu below Sampling Frequency. Then, click Import Filter, and now you can check the frequency response of the fixed-point FIR filter. Analyze if it adjusts to the requirements and compare it with the frequency response obtained from the floating-point design.

• Implement the direct-form FIR filter in Simulink (DSP Builder) using the 16-bit fixed-point coefficients (stored in variable ‘Num1’) and the A/D and D/A converters onboard. You can use the design given in the Tutorial [1] as a template. The A/D converter must convert a sinusoidal signal to digital, filter it and finally convert it back to analog through the D/A converter. Place a SignalTap II node at the output of the filter to acquire the signal later (see Figure 1, Tutorial [1]).

NOTE: filter symmetry property should be used for efficiency.

• Simulate your filter in Simulink with different sinusoidal signals from 100KHz to 5MHz, observing at which frequencies the amplitude of the signal drops significantly.

• Introduce two additional sampling rates to your design: 2.5MHz and 10MHz, similarly to Laboratory 1. Make two additional copies of the FIR filter working at 80MHz and make the necessary modifications to have them working at 2.5MHz and 10MHz, respectively. A PLL and Tsamp blocks are required to modify the sampling rate in these two additional circuits. At the output, use a Multiplexer and two switches from SW3 onboard to select in real-time one out of the three FIR filter outputs (default at 80MHz, 10MHz and 2.5MHz).

NOTE: the Down Sampling block would not be required since the FIR filter is not combinatorial. Also, to simplify the complexity of the circuit, it is suggested to use HDL Subsystem blocks from the Altera DSPBuilder library. Each of these blocks would contain one of the three FIR filters working at different sampling rates.

• Simulate your filter in Simulink with different sinusoidal signals from 100KHz to 5MHz, observing at which frequencies the amplitude of the signal is attenuated significantly. Observe the effects of the different sampling rates, and particularly, observe the aliasing effect when working with the filter sampled at 2.5MHz.

• Download your design to the Stratix EP1S25 DSP development board and test it by using a signal generator and oscilloscope for the different sampling rate options. Observe the effect of the different sampling rates, and again, observe the aliasing effect when working with the circuit sampled at 2.5MHz.

• Capture and analyze the data using SignalTap II Analysis. Show when the amplitude of the signals drops significantly. Also, using the data acquired through SignalTap II, show the frequency response using ‘fft’ command, choosing one input signal frequency in pass-band and another in stop-band.

NOTE: you can use similar commands to those given in the Tutorial [1], “Importing the data acquired from the board in Matlab workspace”, step2.

• What happened if the bit-precision of the coefficients is fixed at 8 bits? Modify your design by replacing 16-bit coefficients by 8-bit coefficients (convert your floating-point coefficients to 8 bits), and observe the behavior. Acquire the data through SignalTap II and show the frequency response using ‘fft’ command, choosing one input signal frequency in pass-band and another in stop-band.

V. Submission

You must submit:

· The Simulink model (.mdl) developed using DSP Builder blocks with three FIR filters working at three different sampling rates: 2.5, 10MHz and 80MHz, respectively.

· A report describing the general procedure to accomplish the project requirements, and answering questions regarding filtering and aliasing effects as detailed in the requirements. The report must contain graphs showing the filtered signals and frequency response using ‘fft’ command for different sampling rates. Choose one input signal frequency in pass-band and another in stop-band. Also, discuss and show the effects of reducing the coefficient precision from 16 to 8 bits.

VI. References

[1] Tutorial: “A/D and D/A Conversion on Altera Stratix EP1S25 Development Board using Simulink and DSP Builder”.

[2] Laboratory 1: Real-Time Implementation for Observing Quantization Effects

[3] DSP Builder User Guide, ver. 5.1.0, Altera, 2005. Local copy at: c:\altera\61\DSPBuilder\Doc\ug_dspbuilder.pdf

[4] Stratix EP1S25 DSP Development Board Data Sheet, ver. 1.6, Altera, 2004. Available online at: http://www.altera.com/literature/ds/ds_stratix_dsp-board-starter.pdf or local copy at: c:\altera\61\kits\stratix_dsp_kit-v1.3.0\Docs\ds_stratix_dsp_bd.pdf

[5] DSP Development Kit Stratix & Stratix Professional Edition (Getting Started User Guide), ver. 1.3.0 rev. 1, Altera, 2004. Available online at: http://www.altera.com/literature/ug/ug_stratix_dsp_kit.pdf or local copy at: c:\altera\61\kits\stratix_dsp_kit-v1.3.0\Docs\ug_stratix_dsp_kit.pdf

DSP LAB-1

LAB 1. Signals in Matlab

Introduction

This lab will describe how to use Matlab for some basic signal representation and manipulation:

• Creating and importing signals

• Sampling and resampling

• Signal visualization

• Modeling noise

• Modulation

Discrete Signals

Time base: t = [0.0 0.1 0.2 0.3]

Signal data: x = [1.0 3.2 2.0 8.5]

The central data construct in Matlab is the numeric array, an ordered collection of real or complex numeric data with one or more dimensions. The basic data objects of signal processing (one-dimensional signals or sequences, multichannel signals, and two-dimensional signals) are all naturally suited to array representation. Matlab represents ordinary one-dimensional sampled data signals, or sequences, as vectors. Vectors are 1-by-n or n-by-1 arrays, where n is the number of samples in the sequence. One way to introduce a sequence into Matlab is to enter it as a list of elements at the command

prompt. The statement

x = [1 2 3 4 5]

creates a simple five-element real sequence in a row vector. It can be converted to a column vector by taking the transpose:

x = [1 2 3 4 5]’

Column vectors extend naturally to the multichannel case, where each channel is represented by a column of an array.

c 2006GM

Another method for creating vector data is to use the colon operator. Consider a 1-second

signal sampled at 1000 Hz. An appropriate time vector would be

t = 0:1e-3:1;

where the colon operator creates a 1001-element row vector representing time from zero to one second in steps of one millisecond.

You can also use linspace to create vector data:

t = linspace(0,1,1e3);

creates a vector of 1000 linearly spaced points between 0 and 1.

Try:

t1 = [0 .1 .2 .3];

t2 = 0:0.1:0.3;

t3 = linspace(0, 0.3, 4);

T = [t1’ t2’ t3’];

X = sin(T)

Q: What does this code show?

Sampling Signals

Analog signal sources include electromagnetic, audio, sonar, biomedical and others. Analog signals must be sampled in order to be processed digitally.

Sampling

x(n) = xa(nTs)

x is a discrete signal sampled from the analog signal xa with a sample period of Ts and a

sample frequency of Fs = 1/Ts.

Try:

Fs = 100;

N = 1000;

stoptime = 9.99;

t1 = (0:N-1)/Fs;

t2 = 0:1/Fs:stoptime;

x1 = sin(2*pi*2*t1);

x2 = sin(2*pi*3*t2);

plot(x1)

figure, plot(x2)

An alternative to creating signals is to use a toolbox function. A variety of toolbox functions generate waveforms. Each of them requires that you begin with a vector representing a time base. Some of these functions will be described later in this lab.

Aliasing

Digital signals are often derived by sampling a continuous-time signal with an analog-to digital (A/D) converter. If the continuous signal, xa(t), is bandlimited, meaning that it does not contain any frequencies higher than a maximum frequency fM, the Shannon sampling theorem says that it can be completely recovered from a set of samples if the sampling frequency fs is greater than two times the maximum frequency of the signal to be sampled:

Fs > 2fM

This maximum frequency fM is known as the Nyquist frequency. If the sampling frequency is

not greater than two times the Nyquist frequency, the continuous signal cannot be uniquely recovered and aliasing occurs. (You heard examples of aliased signals in Homework No.1).

fs > 2fM: Original signal and sampled signal have the same frequency.

fs _ 2fM: Sampled signal is aliased to half the original frequency.

Try:

t = 0:0.001:2;

xa = sin(2*pi*5*t);

plot(t,xa)

hold on

fs = 15;

ts = 0:1/fs:2;

xs1 = sin(2*pi*5*ts);

plot(ts,xs1,’ro-’)

fs = 7.5;

ts = 0:1/fs:2;

xs2 = sin(2*pi*5*ts);

plot(ts,xs2,’ro-’)

hold off

Q: What is the frequency of xs2?// (There is aliasing here. We need sampling theory.

However can use the fft function on the signal to determine the frequency).

Signal Visualization

• View signal amplitude vs. time index

• Functions: plot, stem, stairs, strips

• Listen to data: sound

Note: the sound and soundsc commands will not work if your computer hardware isn’t set up. If that is the case, view the signals instead of listening to them.

Try:

t = [0.1 0.2 0.3 0.4];

x = [1.0 8.0 4.5 9.7];

plot(t,x)

figure, stem(t,x)

figure, stairs(t,x)

fs = 1000;

ts = 0:1/fs:2;

f = 250 + 240*sin(2*pi*ts);

x = sin(2*pi*f.*ts);

strips(x,0.25,fs)

sound(x,fs)

plot(ts,x)

plot(ts(1:200),x(1:200))

Q: What does the strips command do? (See ’help strips’.)

Q: What does the .* operator do?

Signal Processing Tool

The Signal Processing Toolbox application, SPTool, provides a rich graphical environment for signal viewing, filter design, and spectral analysis. You can use SPTool to analyze signals, design filters, analyze filters, filter signals, and analyze signal spectra. You can accomplish these tasks using four GUIs that you access from within SPTool:

• The Signal Browser is for analyzing signals. You can also play portions of signals using

your computer’s audio hardware.

• The Filter Designer is for designing or editing FIR and IIR digital filters. Note that the

FDATool is the preferred GUI to use for filter designs. FDATool is discussed in later labs.

• The Filter Viewer is for analyzing filter characteristics.

• The Spectrum Viewer is for spectral analysis.

Open SPTool by typing sptool at the command prompt.

Try:

sptool

Look at the train signal, FIRbp filter, and trainse spectrum. (You see 3 panes - Signals, Filters, Spectra. Filter Designer is available through File 7! Preferences. You can play sounds using the LS icon. When viewing spectra, note that many methods of determining spectra, including the fft, are available.)

Importing a Signal

You can use SPTool to analyze the signals, filters, or spectra that you create at the Matlab

command line. You can import signals, filters, or spectra from the Matlab workspace into the SPTool workspace using the Import item under the File menu.

Try:

fs = 1000;

ts = 0:1/fs:0.5;

f = 250 + 240*sin(2*pi*ts);

x = sin(2*pi*f.*ts);

Import these signals (f and x) into the SPTool and use the tool to examine them.

Q: What are the icons to use for horizontal zoom?

Try zooming in using the mouse.

Signal Browser

The Signal Browser tool is an interactive signal exploration environment. It provides a

graphical view of the signal object(s) currently selected in teh Signals list of SPTool.

Using the Signal Browser you can

• View and compare vector/array signals

• Zoom in on a range of signal data to examine it more closely

• Measure a variety of characteristics of signal data

• Play signal data on audio hardware

To open/activate the Signal Browser for the SPTool,

• Click one or more signals (use the Shift key for multiple selections) in the Signals list of

SPTool.

• Click the View button in the Signals list of SPTool.

Changing Sample Rates

To change the sample rate of a signal in SPTool,

1. Click a signal in the Signals list in SPTool.

2. Select the Sampling frequency item in the Edit menu.

3. Enter the desired sampling frequency and cliick OK.

Try changing the sampling rate of the imported signal.

Signal Generation

Signals

• Create a time base vector

t = [0:0.1:2];

• Create a signal as a function of time

x = sin(pi*t/2);

plot(t,x)

Useful Matlab functions

• Nonperiodic functions

ones, zeros

• Periodic functions

sin, cos, square, sawtooth

Nonperiodic Signals

t = linspace(0,1,11)

• Step:

y = ones(1,11);

stem(y)

• Impulse:

y = [1 zeros(1,10)];

stem(y)

• Ramp:

y = 2*t;

plot(y)

Useful Matlab functions

step, impulse, gensig

Try:

Step function:

fs = 10;

ts = [0:1/fs:5 5:1/fs:10];

x = [zeros(1,51) ones(1,51)];

stairs(ts,x)

Impulse function with width w:

fs = 10;

w = 0.1;

ts = [-1:1/fs:-w 0 w:1/fs:1];

x = [zeros(1,10) 1 zeros(1,10)];

plot(ts,x)

Delta function:

ts = 0:0.5:5;

x = [1 zeros(1,length(ts)-1)];

stem(ts,x)

axis([-1 6 0 2])

Sinusoids

Sinusoid parameters

• Amplitude, A

• Frequency, f

• Phase shift, _

• Vertical offset, B

The general form of a sine wave is

y = Asin(2_ft + _) + B

Example: generate a sine wave given the following specifications:

• A = 5

• f = 2 Hz

• _ = _/8 radians

t = linspace(0,1,1001);

A = 5;

f = 2;

p = pi/8;

sinewave = A*sin(2*pi*f*t + p);

plot(t, sinewave)

Try:

edit sine_wave

sine_wave

edit sinfun

[A T] = sinfun(1,2,3,4)

Square Waves

Square wave generation is like sine wave generation, but you specify a duty cycle, which is the percentage of the time over one period that the amplitude is high.

Example:

• duty cycle is 50% (the Matlab default)

• frequency is 4 Hz.

t = linspace(0,1,1001);

sqw1 = square(2*pi*4*t);

plot(t,sqw1)

axis([-0.1 1.1 -1.1 1.1])

Example:

• duty cycle is 75%

• frequency is 4 Hz.

t = linspace(0,1,1001);

sqw2 = square(2*pi*4*t,75);

plot(t,sqw2)

axis([-0.1 1.1 -1.1 1.1])

Sawtooth Waves

Sawtooth waves are like square waves except that instead of specifying a duty cycle, you

specify the location of the peak of the sawtooth.

Example:

• peak at the end of the period (the Matlab default)

• frequency is 3 Hz.

t = linspace(0,1,1001);

saw1 = sawtooth(2*pi*3*t);

plot(t,saw1)

Example:

• peak is halfway through the period

• frequency is 3 Hz.

t = linspace(0,1,1001);

saw2 = sawtooth(2*pi*3*t,1/2);

plot(t,saw2)

Complex Signals

Periodic signals can be represented by complex exponentials:

x(t) = ej2_ft = cos(2_ft) + jsin(2_ft) = cos(t) + jsin(t)

If t is measured in seconds, then f will have units of sec−1, and will have units of radians/ second. In signal processing, we associate the unit circle with one sampling cycle, so that a sampling frequency of Fs is associated with 2_ radians, and the Nyquist frequency Fs/2 is associated with _ radians. Values of in the upper half-plane, in units of Hz, then correspond to frequencies within the sampled signal. In Matlab, type:

x = exp(2*pi*j*f*t);

plot(x)

Matlab recognizes either j or i as the square root of -1, unless you have defined variables j

or i with different values. Useful Matlab functions real, imag, abs, angle

Try:

edit zsig

zsig(5)

Look at both figures and describe what you see.

Importing Data

An important component of the Matlab environment is the ability to read and write data

from/to external sources. Matlab has extensive capabilities for interfacing directly with data from external programs and instrumentation.

In this lab, we concentrate on reading and writing data that has already been stored in external files. Files come in a variety of standard formats, and Matlab has specialized routines for working with each of them. To see a list of supported file formats, type:

help fileformats To see a list of associated I/O functions, type:

help iofun

Matlab provides a graphical user interface, the Import Wizard, to the various I/O functions. You access the Wizard by choosing File! Import Data or by typing:

Uiimport The Matlab command importdata is a programmatic version of the Wizard, accepting all of the default choices without opening the graphical user interface. You can use importdata in M-files to read in data from any of the supported file formats. Matlab also has a large selection of low-level file I/O functions, modeled after those in the C

programming language. These allow you to work with unsupported formats by instructing Matlab to open a file in memory, position itself within the file, read or write specific formatted data, and then close the file.

Try:

help fileformats

help iofun

jan = textread(’all_temps.txt’,’%*u%u%*[^\n]’,’headerlines’,4);

[data text] = xlsread(’stockdata.xls’);

plot(data(:,2))

legend(text{1,3})

Explain how the colon operator works in the preceding plot command.

I = importdata(’eli.jpg’);

image(I)

which theme.wav

uiimport

Browse for:

theme.wav

soundsc(data,fs)

Save and Load

Two data I/O functions are especially useful when working with Matlab variables.

• The save command writes workspace variables to a binary Matlab data file (MAT-file)

with a .mat extension. The file is placed in the current directory.

• The load command reads variables from a MAT-file back into the Matlab workspace.

Although quite specialized, save and load can be used for day-to-day management of your

Matlab computations.

Try:

doc save

doc load

t = 0:0.1:10;

x1 = sin(t);

x2 = sin(2*t);

x3 = sin(3*t);

save myvars

clear

load myvars t x3

Note the list of variables in the workspace tab in the upper left of the Matlab window.

Modeling Noise

To model signals in space, in the atmosphere, in sea water, or in any communications channel, it is necessary to model noise.

Matlab has two functions for generating random numbers, which can be added to signals to model noise.

Uniform random numbers

A = rand(m,n);

generates an mxn array of random numbers from the uniform distribution on the interval

[0,1]. To generate uniformly distributed random numbers from the interval [a,b], shift and

stretch:

A = a + (b-a)*rand(m,n);

Gaussian random numbers

A = randn(m,n);

generates an mxn array of random numbers from the standard normal distribution with

mean 0 and standard deviation 1. To generate random numbers from a normal distribution with mean mu and standard deviation sigma, shift and stretch:

A = mu + sigma*rand(m,n);

Random numbers from other distributions

Random numbers from other distributions can be generated using the uniform random number generator and knowledge of the distribution’s inverse cumulative distribution function. Random number generators for several dozen common distributions are available in the Statistics Toolbox.

Adding Noise to a Signal

noisy signal = signal + noise

y1 = x + rand(size(x)) % uniform noise

y2 = x + randn(size(x)) % Gaussian noise

Example:

Add Gaussian noise to middle C.

fs = 1e4;

t = 0:1/fs:5;

sw = sin(2*pi*262.62*t); % middle C

n = 0.1*randnsize(sw);

swn = sw + n:

Try:

edit noisyC

noisyC

strips(swn, .1,1e4)

Zoom in on the strips plot. (Note: you might have to cut and paste from the noisyC script

to generate swn.)

Pseudorandomness

This number:

0.95012928514718

is the first number produced by the Matlab uniform random number generator with its

default settings. Start up Matlab, set format long, type rand, and you get the number.

If all Matlab users, all around the world, all on different computers, keep getting this same number, is it really “random”? No, it isn’t. Computers are deterministic machines and should not exhibit random behavior. If your computer doesn’t access some external device, like a gamma ray counter or a clock, then it must really be computing pseudorandom numbers.

A working definition of randomness was given in 1951 by Berkeley professor D. H. Lehmer, a pioneer in computing and, especially, computational number theory:

A random sequence is a vague notion ... in which each term is unpredictable

to the uninitiated and whose digits pass a certain number of tests traditional with

statisticians ... Random number generators proceed deterministically from their current state. To view the current state of rand, type:

s = rand(’state’)

This returns a 35-element vector containing the current state.

To change the state of rand:

rand(’state’,s) Sets the state to s.

rand(’state’,0) Resets the generator to its initial state.

rand(’state’, sum(100*clock)) Sets to a new state each time.

Commands for randn are analogous.

Try:

s = rand(’state’)

format long

rand

rand(’state’,sum(100*clock))

s = rand(’state’)

format long

rand

Resampling

The Signal Processing Toolbox provides a number of functions that resample a signal at a

higher or lower rate.

y = downsample(x,n)

decreases the effective sampling rate of x by keeping every nth sample starting with the first sample. x can be a vector or a matrix. If x is a matrix, each column is considered a separate

sequence.

y = upsample(x,n)

Increases the effective sampling rate of x by inserting n−1 zeros between samples. x can be a vector or a matrix. If x is a matrix, each column is considered a separate sequence. The upsampled y has x _ n samples.

y = resample(x,p,q)

resamples the sequence in vector x at p/q times the original sampling rate, using a polyphase filter implementation. p and q must be positive integers. The length of y is equal to ceil(length(x) _ p/q). If x is a matrix, resample works down the columns of x.

y = interp(x,r)

increases the sampling rate of x by a factor of r. The interpolated vector y is r times longer than the original input x.

y = decimate(x,r)

reduces the sampling rate of x by a factor of r. The decimated vector y is r times shorter

in length than the input vector x. By default, decimate employs an eighth-order lowpass

Chebyshev Type I filter. It filters the input sequence in both the forward and reverse

directions to remove all phase distortion, effectively doubling the filter order.

Try:

load mtlb

sound(mtlb,Fs)

mtlb4 = downsample(mtlb,4)

mtlb8 = downsample(mtlb,8)

sound(mtlb8,fs/8)

What are the sizes of mtlb, mtlb4, and mtlb8?

(If sound doesn’t work, plot the signals.)

t = 0:0.00025:1;

x = sin(2*pi*30*t) + sin(2*pi*60*t);

y = decimate(x,4);

subplot(211), stem(x(1:120))

axis([0 120 -2 2])

title(’Original Signal’)

subplot(212), stem(y(1:30))

title(’Decimated Signal’)

Modulation and Demodulation

Modulation varies the amplitude, phase, or frequency of a carrier signal with reference to a message signal.

The Matlab modulate function modulates a message signal with a specified modulation

method. The syntax is

y = modulate(x,fc,fs,’method’)

where:

• x is the message signal.

• fc is the carrier frequency.

• fs is the sampling frequency.

• method is a flag for the desired modulation method (see table below).

Method Description

amdsb-sc or am Amplitude modulation, double side-band, suppressed carrier

amdsb-tc Amplitude modulation, double side-band, transmitted carrier

amssb Amplitude modulation, single side-band

fm Frequency modulation

pm Phase modulation

ppm Pulse position modulation

pwm Pulse width modulation

qam Quadrature amplitude modulation

The demod function performs demodulation, that is, it obtains the original message signal

from the modulated signal. The syntax is:

x = demod(y,fs,fs,’method’)

demod uses any of the methods shown for modulate. The signal x is attenuated relative to

y because demodulation uses lowpass filtering.

Exercise: High and Low

1. Create a signal equal to the sum of two sine waves with the following characteristics:

• 3-second duration

• Sampling frequency = 2 kHz

• Sinusoid 1: frequency 50 Hz (low), amplitude 10, phase = 0

• Sinusoid 2: frequency 950 Hz (high), amplitude 1, phase = 0

2. View and listen to the signal using an M-file

3. Import the signal into SPTool and view it. Listen to the signal.

HAND IN:

ANSWERS TO ALL QUESTIONS STARTING WITH THE LETTER Q IN

FRONT OF IT.

(The material in this lab handout was put together by Paul Beliveau and derives principally from the MathWorks training document “MATLAB for Signal Processing”, 2006.)

Try out these examples