讲解：DFT、Matlab、Matlab、FFTSQL|Matlab

IntroductionThis lab is a revision of the Discrete Fourier Transform (DFT), and theFast Fourier Transform (FFT), and an introduction to the Short-Time Fourier Transform (STFT) and thespectrogram.The outcomes from the lab are to be handed in as a “folder” of results, showing thatyou have completed the steps of the lab successfully. A box in the right-handmargin indicates where an outcome is expected – like this:The Matlab programs you write will be short: you can print them out if you wish, buthand-written listings are OK too. Sketch graphs are also OK.1. Getting StartedIn your home directory, create the subdirectories “EBU6018” and “EBU6018/lab1”.Start Matlab. Use “cd ” to get into the directory “lab1” you have justcreated.2. Discrete Fourier TransformIn Matlab, type “edit” to start the Matlab editor.Create a Matlab function in the file “dft.m” to calculate the Discrete FourierTransform 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” torefer 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 MethodsLab 1: DFT, FFT and STFTDepartment of Electronic Engineering and Computer ScienceExample DFT function outline in MatlabGenerate some waveforms to test your function. Test your dft on at least thefollowing 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](NB: Why doesn’t this “look” symmetrical? Remember that the DFT repeats, sothe 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 anglesubplot(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);endfunction sw = dft(st)% DFT - Discrete Fourier TransformM = length(st);N = M;WN = exp(2*pi*j/N);% Main loopfor 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 functionMatlab has a built-in Fast Fourier Transform, “fft”.Compare the results of your dft against the built-in fft. Are the results the same? Ifso, why: if not, why not?Find out the complexity of your dft and the built-in fft, i.e. how long they take toperform their calculation for various lengths of s. Use “tic” and “toc” to measure thetime taken to perform the operation, so e.g. tic; dft(ones(1,4)); toc % No “;” for final expressionwill 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. Ifso, 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(n2) 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. Ifyou like, you could write this as a recursive function (one that calls itself) – see theoutline below.Plot and compare its speed to the DFT, showing that your “my_fft” function takesO(n log n) time rather than O(n2) timeDerivation 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 usingeven oddNotes:(1) The above only works if N is a power of 2 (64, 128, 1024, etc), so your programmay 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-1Example outline of Matlab function to calculate FFT:function sw = my_fft(st);% Recursive Implementation of Fast Fourer TransformN = length(st);% check length of N is 2^kif (rem(log(N),log(2))) disp(slow_fft: N must be an exact power of 2) returnendWN = exp(2*pi*j/N);% split st into even and odd samplesst_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];endsw = gg+hh; 4. Single Windowed Fourier TransformSave one of the audio files on the course details page athttps://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 betterthan “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 sectionof the signal. To smooth out end effects, we will use a “Hanning” window to multiplythe segment that we select. You can show the Hanning window of length 256 inMatlab using “plot(hanning(256))”.Construct a Matlab function in the file “wft.m” that will select a section from a fileand 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 windowlength. You might use the following steps:1) Select the desired section from the signal, for example usings(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 forvarious values of t and n (note that values of t near the beginning and end of s maycause an error, depending on how clever you were at step (1)). Try also plotting with alog y-scale. Explain the difference between these results*[OMIT][Optional]: Make a matlab m-file that loops through different values for t in steps ofe.g. 50, using “pause” between each step.5. STFT and SpectrogramNow we will construct a “spectrogram” to visualize the time-frequency information ina signal on one image.Read the Matlab documentation for the Matlab “specgram” function (try “helpspecgram” for information).Using specgram, investigate the audio files on the course details page athttps://www.student.elec.qmul.ac.uk/courseinfo/EBU6018/Try different window sizes (“NFFT”) to see the effect. For fastest results on longfiles, use powers of 2 (Why?). Record what values of window size give bestvisualization results for different files, and suggest why. 5.1 Analysis of Piccolo soundFrom the course webpage download ‘piccolo.wav’ and load it into Matlab using:[x fs] = wavread(‘piccolo.wav’); % fs = sampling frequencyRecord the sampling frequency, fs.If you have headphones, try listening to the signal, usingsoundsc(x,fs); %fs is the sampling frequency of xPlot a spectrogram of x, using the ‘specgram’ function.From the spectram plot, estimate the fundamental frequencies (f0) of the 3 notes inthe 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 samplingfrequency (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 estimationExplain 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 internetusing your spectrogram. *[OMIT]5.1 DIY STFT and Spectrogram [Optional, but highlyrecommended]Construct a Matlab function “sg(s,N)” in a file called “sg.m” to compute aspectrogram of a waveform s with window size N (NFFT in Matlab’s specgram).To do this, your function shouldi) divide the signal “s” into sections of length N,ii) multiply s by a Hanning windowiii) perform an FFT of each section, andiv) construct a matrix where each column is the absolute value of one FFTHints:• 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, soyou 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 usingM(:,n) = yPlot usingimagesc(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 thedefault operation of “specgram”.6. Handing InCompile the answers to the exercises, including the answers to specific questions,program listings (including comments), and plots from experiments, into a “folder” ofresults showing that you have completed the lab, and submit electronically. You donot need to write a formal report.IMPORTANT: Plagiarism (copying from other students, or copying the work ofothers 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, MEPDModified ARW for EBU6018.转自：http://www.daixie0.com/contents/12/4346.html