Lab 1. 9/8/98 Due 9/22/98.

Exercises to be done using MATLAB. Two PCs will be available in 87 MSC.

Problem 1. In order to familiarize everyone with the use of the discrete Fourier transform:
generate y(t) = 10*sin(2*pi*f*t) where f = 100Hz. Make the signal 1 second long. Use a power of 2 for the number of samples. Note: the Nyquist theorem states you have to sample a signal at least twice as fast as the highest frequency present. It is often better to sample at least 5Xs faster, especially, if you want a nice looking waveform to plot.

Plot the waveform, and its Fourier magnitude and phase.

There will be an MATLAB m-file, L1P1, available as a starter for those who have difficulty with this exercise. It can be used as a template for several of the following exercises. You can display any m-file in MATLAB by typing dbtype FileName. This can also be printed by selecting the lines on the screen you want to type with the mouse and then use the File command: followed by print selection.

Problem 2. Clip the waveform from Problem 1 at the following levels: 100%, 90%, 75%, 50%, 25% and 10 % of the maximum & minimum signal level, i.e., this is a symmetric clipping function. What amount of the total signal energy is present in the each of the first three harmonics for each three condition. Make a table and/or graph of the results. Note that as the clipping level gets lower that the signal is essentially a square wave.

Problem 3. Distortion of signals results in the generation of components that are at frequencies that are not present in the input signal. One can represent the distortion process with a polynomial function. If x(t) is the input and y(t) the output then we can write y = aF x + bF x² +cF x³ +…

set a = 0, b= 1, and c=0. Let the input consist of two sinusoids with frequencies of 1000 and 1100Hz.

What frequency components are present in the output signal?

Note that you can do this analytically. It just trigonometry.

Remember: Phythagorus: the Square of the hypothenous = the sum of the squares of the sides of a right triangle.

sin²(x) + cos²(x) = 1.

sin(A ± B) = sin(A)cos(B) ± sin(B)cos(A)

cos(A ± B) = cos(A)cos(B) -+ sin(B) sin(A)

and if you set A = B = x you can derive:

cos²(x) = (1 + cos(2x))/2. You can verify your result using MATLAB.

Problem. 4. With the same frequency components as in Problem 3, set a = 0, b = 0 and c = 1.

What frequency components are present in the output. Analytically or with MATLAB FFT.

Problem 5. Same as Problem 4 but set a = 1 and c = .05.

What is the spectrum? Plot it.

Problem 6. A special case occurs when the ratio of the two frequency components is set to 1.5 (f2/f1 = 1.5, f1 = 1000Hz). Set a = 1, b = c = 0.05. Vary the phase angle between the two primary frequencies and plot the magnitude of the difference frequency, f2-f1 and the best known distortion product, 2f1-f2, often labeled the cubic difference tone (CDT). We will discuss these distortion products and their generation in the auditory system more later. Use MATLAB to obtain the results.

Introduction to Matlab 5.2

Double click on the Matlab 5.2 icon.

At the Matlab prompt >>

Type demo

Click on matrices (this is all Matlab really understands, i.e., matrices).

Then Run Basic matrix.

Continue through the slides.

When done click on close.

Then select numerics.

Select Fast fourier transform

Run fast fourier transform

Sequence the slides.

At end close.

Your now ready to enter a set of instructions for Matlab to execute.

Click on File on the toolbar.

Select New and then m-file.

You are now in an editor/debugger.

Type in a series of instructions.

e.g.

x=1:2:20; % says set the x variable to 1,3,5,7,… 19

y=x .^2; % says make y= each element of x squared.

plot(x,y)

Pause; % to advance beyond this point hit any key

plot(x,y,’*r’)

%save this file with a name such as temp in your directory

%then close the editor; this returns you to Matlab

type cd c:\neuro730\yourname

temp %this will run the series of instruction that you just saved as an m-file.

%L1P1c.m, W. Rhode, 8/29/98
%generate a sinewave and fourier analyze it.
%time resolution: tresol determines the sample rate
clear all;
tresol=1/512; %note a power of two was chosen for the denominator
%the signal will be 'sampled' at time =t, range of 1 second
t=0:tresol:1-tresol; %we want exactly 512 samples in the 1 second period.
tmax = 1;
ampl=10;
%generate the sinewave with a frequency of 100Hz
f=100;
y=ampl*sin(2*pi*f*t +pi/2);%note the phase is set to pi/2.
n=1/tresol;
fy = fft(y,n);
%the FFT results are in complex form, a + jb.
% you can obtain the magnitude by taking the absolute value of the output
magfy= abs(fy);
%setup the frequency vector for the spectral analysis components.
freqres = 1/tmax;%since tmax = 1sec, the resolution of the resulting spectra is 1Hz.
freq = 0:1/freqres:(n/2-1)/freqres;%this give exactly 256 spectral components
figure
hold on;
subplot(1,3,1);%the subplot function allows several plots per page
%(this says 1row; 3 cols); subplot(4,2,i) says(the ith subplot of 4 rows, 2 cols);
plot(t,y)
axis([0 .1 -10 10]);%plot only .1 of the waveform. 0<x<.1; -10<y<10
subplot(1,3,2);
plot(freq,2*magfy(1:256)/n);% the magfy had to be normalized by the number of sample points.
% it has to be multiplied by 2
%since we are discarding the upper 256 points
xlabel 'frequency in Hz';
ylabel 'magnitude';
title 'fourier analysis of a sinewave';
yphase = angle(fy) + pi/2;%The phase is returned in cosine, change to sine by adding pi/2
yinsig = find(magfy < .00001); %locate the frequency components that are insignificant
yphase(yinsig) = yphase(yinsig)*0;
subplot(1,3,3);
plot(freq,yphase(1:n/2));
xlabel 'frequency in Hz';
ylabel 'phase in radians';