|
代做DFT留学生作业、代写Matlab实验作业、Matlab程序语言作业调试、FFT课程作业代做
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.
因为专业，所以值得信赖。如有需要，请加QQ：99515681或邮箱：99515681@qq.com 微信：codehelp
|
