代做DFT留学生、代写Matlab代写、Matlab程序语言调试、FFT代做

日期：2019-11-22 10:23

Introduction

This lab is a revision of the Discrete Fourier Transform (DFT), and the

Fast Fourier Transform (FFT), and an introduction to the Short-Time Fourier Transform (STFT) and the

spectrogram.

The outcomes from the lab are to be handed in as a “folder” of results, showing that

you have completed the steps of the lab successfully. A box in the right-hand

margin indicates where an outcome is expected – like this:

The Matlab programs you write will be short: you can print them out if you wish, but

hand-written listings are OK too. Sketch graphs are also OK.

1. Getting Started

In your home directory, create the subdirectories “EBU6018” and “EBU6018/lab1”.

Start Matlab. Use “cd <directory>” to get into the directory “lab1” you have just

created.

2. Discrete Fourier Transform

In Matlab, type “edit” to start the Matlab editor.

Create a Matlab function in the file “dft.m” to calculate the Discrete Fourier

Transform of a signal. Recall that the DFT is given by [Qian, eqn (2.34)]

[NB: The “j” is missing in Qian’s definition of WN on p33.]

Hints:

• Start your function with

function sw = dft(st)

so “st” is the time waveform vector, and “sw” is the frequency waveform vector

• Matlab vectors (e.g. st and sw) start from 1, not zero, so use “n-1” and “m-1” to

refer to the appropriate element

• Assume that N=M, and use the “length(st)” to find the value to use for these.

An example outline for your Matlab function is provided below.

EBU6018 Advanced Transform Methods

Lab 1: DFT, FFT and STFT

Department of Electronic Engineering and Computer Science

Example DFT function outline in Matlab

Generate some waveforms to test your function. Test your dft on at least the

following signals:

• Uniform function: “s=ones(1,64);”

• Delta function: “s = ((1:64)= =1);”

[NB: “1:64” generates the vector (1 2 … 64) ].

• Sine wave: “s = sin(((1:64)-1)*2*pi*w/100)” for various values of w.

Why do we need to use “(1:64)-1”?

What values of w give the cleanest dft?

What happens if we use “cos”?

• Symmetrical rectangular pulse: “s = [0:31 32:-1:1]<T” for various values of T.

(NB: Why doesn’t this “look” symmetrical? Remember that the DFT repeats, so

the time interval 32 .. 63 is “the same as” the interval -31 .. -1).

The following function may be useful to display your results:

If you want zero frequency (or time) to appear in the middle of your plot, use

“fftshift”, e.g. “stem4(fftshift(dft(s)));”

Explain your results in terms of what you know about the Fourier Transform.

function stem4(s)

% STEM4 - View complex signal as real, imag, abs and angle

subplot(4,1,1); stem(real(s)); title('Real');

subplot(4,1,2); stem(imag(s)); title('Imag');

subplot(4,1,3); stem(abs(s)); title('Abs');

subplot(4,1,4); stem(angle(s)); title('Angle');

end

function sw = dft(st)

% DFT - Discrete Fourier Transform

M = length(st);

N = M;

WN = exp(2*pi*j/N);

% Main loop

for n=0:N-1

temp = 0;

for m=0:M-1

[** Do something useful here **]

end

sw(n+1) = temp;

end

3. Comparison with Matlab’s FFT function

Matlab has a built-in Fast Fourier Transform, “fft”.

Compare the results of your dft against the built-in fft. Are the results the same? If

so, why: if not, why not?

Find out the complexity of your dft and the built-in fft, i.e. how long they take to

perform their calculation for various lengths of s. Use “tic” and “toc” to measure the

time taken to perform the operation, so e.g.

tic; dft(ones(1,4)); toc % No “;” for final expression

will report how long a 4-point DFT took to calculate.

Hint: You may find your dft is too fast for tic/toc to measure any useful difference. If

so, run it several times, e.g.

tic; for (i=1:1e4) dft(ones(1,4)); end; toc

(Of course, remember to divide your measure by the number of times round the loop!)

Make a log-log plot (using “loglog”) showing the time increase with the size n of s.

On your plot, show that the DFT takes O(n

2

) time, while the FFT takes O(n log n).

Hint: Use “hold on” if you want to add a second “loglog” plot to an existing plot.

Explain what this tells you about the DFT compared to the FFT in real applications,

i.e. as n gets larger.

*[OMIT]3.1 DIY-FFT [Optional, but highly recommended] [OMIT]

Write a Matlab function (called e.g. “my_fft”) to calculate the FFT of a signal. If

you like, you could write this as a recursive function (one that calls itself) – see the

outline below.

Plot and compare its speed to the DFT, showing that your “my_fft” function takes

O(n log n) time rather than O(n2

) time

Derivation of the FFT:odd the of FFT point- the is and

the is where samples, even the ofFFT point

:get weodd for and even for using

even odd

Notes:

(1) The above only works if N is a power of 2 (64, 128, 1024, etc), so your program

may not work if you use other lengths of s (you could check this, if you like!)

(2) Note that a 1-point FFT of a signal is the signal itself, so the 1-point FFT is easy

(to be sure of this, check the DFT formula with N=1).

(3) Remember that Matlab vectors start at 1 (not zero), so go from 1...N not 0…N-1

Example outline of Matlab function to calculate FFT:

function sw = my_fft(st);

% Recursive Implementation of Fast Fourer Transform

N = length(st);

% check length of N is 2^k

if (rem(log(N),log(2)))

disp('slow_fft: N must be an exact power of 2')

return

end

WN = exp(2*pi*j/N);

% split st into even and odd samples

st_even = st(1:2:end-1);

st_odd = st(2:2:end);

% implement recursion here...

if (N==2)

g = st_even; % = st(0+1)

h = st_odd; % = st(2)

gg = [g g];

hh = [h -h];

else

g = [** Something useful here **];

h = [** Something useful here **];

gg = [g g];

hh = WN.^(-[0:N-1]).*[h h];

end

sw = gg+hh;

4. Single Windowed Fourier Transform

Save one of the audio files on the course details page at

https://www.student.elec.qmul.ac.uk/courseinfo/EBU6018/

into your “lab1” directory.

Read into Matlab, using “s = wavread('file.wav')”.

Where ‘file.wav’ could be ‘dbarrett2.wav’

Plot the magnitude (“abs”) of the FFT of the waveform. (“plot” is probably better

than “stem” for these longer signals). Explain what this tells you about the waveform.

We will now construct a function that will allow you to “zoom in” on a short section

of the signal. To smooth out end effects, we will use a “Hanning” window to multiply

the segment that we select. You can show the Hanning window of length 256 in

Matlab using “plot(hanning(256))”.

Construct a Matlab function in the file “wft.m” that will select a section from a file

and window it. The function “wft” is to be called as follows:

y = wft(s, t, n);

where s is the signal, t is the time in the middle of the window, and n is a window

length. You might use the following steps:

1) Select the desired section from the signal, for example using

s(floor(t-n/2)+(1:n));

(if you don’t see how this works, try “help colon”).

2) Multiply elementwise with a Hanning window of length n, using “.*”

3) Use the built-in Matlab fft function to calculate the DFT.

Plot the magnitude of this single windowed Fourier transform of your signal for

various values of t and n (note that values of t near the beginning and end of s may

cause an error, depending on how clever you were at step (1)). Try also plotting with a

log y-scale. Explain the difference between these results

*[OMIT][Optional]: Make a matlab m-file that loops through different values for t in steps of

e.g. 50, using “pause” between each step.

5. STFT and Spectrogram

Now we will construct a “spectrogram” to visualize the time-frequency information in

a signal on one image.

Read the Matlab documentation for the Matlab “specgram” function (try “help

specgram” for information).

Using specgram, investigate the audio files on the course details page at

https://www.student.elec.qmul.ac.uk/courseinfo/EBU6018/

Try different window sizes (“NFFT”) to see the effect. For fastest results on long

files, use powers of 2 (Why?). Record what values of window size give best

visualization results for different files, and suggest why.

5.1 Analysis of Piccolo sound

From the course webpage download ‘piccolo.wav’ and load it into Matlab using:

[x fs] = wavread(‘piccolo.wav’); % fs = sampling frequency

Record the sampling frequency, fs.

If you have headphones, try listening to the signal, using

soundsc(x,fs); %fs is the sampling frequency of x

Plot a spectrogram of x, using the ‘specgram’ function.

From the spectram plot, estimate the fundamental frequencies (f0) of the 3 notes in

the sample, giving your answers in Hz.

Repeat your estimates for different window sizes.

Notes: You will need to use your window size, (NFFT) and the value for the sampling

frequency (fs) in your calculation. Figure 1 is given as a guide to help you.

Make your calculation in 2 ways:

(1) by calculating the frequency range displayed by specgram, and

(2) by supplying specgram with the correct value fs when you call it.

Check that both of these methods agree.

Figure 1: Angular frequency representation for f0 estimation

Explain what happens to the accuracy of your f0 values as you vary the window size.

For further experimentation, try visualizing other “wav” files available on the internet

using your spectrogram.

*[OMIT]5.1 DIY STFT and Spectrogram [Optional, but highly

recommended]

Construct a Matlab function “sg(s,N)” in a file called “sg.m” to compute a

spectrogram of a waveform s with window size N (NFFT in Matlab’s specgram).

To do this, your function should

i) divide the signal “s” into sections of length N,

ii) multiply s by a Hanning window

iii) perform an FFT of each section, and

iv) construct a matrix where each column is the absolute value of one FFT

Hints:

• You can select the k-th segment of length N using “s( ((1:N)+(k-1)*N) )”

• You can get a Hanning window of length N by using the Matlab function

“W=hanning(N)”. Multiply by a segment s1 using “s1.*W” (dot-star).

• Since the signal is real, you know the FFT result will be Hermitian symmetric, so

you can discard one half of the vector of results.

• You can set the n-th column of a matrix to be a 1xN vector y by using

M(:,n) = y'

Plot using

imagesc(log10(abs(B))); axis xy;

where B is the spectrogram (“axis xy” restores the origin to the bottom).

How should you call “specgram” to get the most similar results to your function “sg”?

Modify your function “sg” so that it overlaps its windows in the same way as the

default operation of “specgram”.

6. Handing In

Compile the answers to the exercises, including the answers to specific questions,

program listings (including comments), and plots from experiments, into a “folder” of

results showing that you have completed the lab, and submit electronically. You do

not need to write a formal report.

IMPORTANT: Plagiarism (copying from other students, or copying the work of

others without proper referencing) is cheating, and will not be tolerated.

IF TWO “FOLDERS” ARE FOUND TO CONTAIN IDENTICAL MATERIAL,

BOTH WILL BE GIVEN A MARK OF ZERO.

Updated by MPD, MEPD

Modified ARW for EBU6018.