Music 421 Spring 2004-2005 Homework #3 More Windows 46 points Due in one week (4/21/2005)
1. In this problem, we consider the frequencies corresponding to points in the DFT (like we did in problem 5 of the last HW), being, as always, aware of the order output by matlab when we use fft(). In all cases here, assume that the corresponding normalized radian frequencies are defined on [−π, π) (so we should never have an answer with frequnecy in Hz above fs /2). (a) (2 pts) Consider an arbitrary time domain signal of windowed length 999 samples. Say that we zero pad the signal to length 1024 samples, and then take the FFT in matlab. Assuming a sampling frequency of 44100, What frequencies in Hz do each of the 1024 points (also called “bins”) represent, assuming that you do not fftshift the fft? (b) (1 pt) Consider the previous item, but where the windowed signal has length 1000 samples, and is still zero padded to total length 1024 samples. What frequencies in Hz do each of the 1024 bins represent, assuming that you do not fftshift the fft? (c) (1 pt) Consider the first item, but with no zero padding so that we take the 999 point transfrom of the 999 point windowed signal. What frequencies in Hz do each of the 999 bins represent, assuming that you do not fftshift the fft? (d) (2 pts) Say that you take the N point DFT of an M point time-windowed signal (recall that N > M implies the use of zero padding), where the sampling frequency is fs . Find (an) expression(s) that gives the frequency corresponding to bin b in Hz, where b is the MATLAB index from 1 to ceil(N/2). You must consider even or odd DFT lengths. Remember that the first frequency matlab returns is always 0, even though its index is 1. (e) (2 pts) Implement the previous item as a function in matlab. Turn in your code. 2. We compare the spectral features of the Hann, Hann-Poisson, and Poisson windows in this problem. When you write a matlab function, turn in your code. (a) (3 pts) Write a matlab function that takes window-length M and window-parameter α as arguments and returns the Poisson window 1 . The output should be a column vector in “causal” form (as in Matlab windows. You may assume M is an odd integer. 1
http://ccrma.stanford.edu/˜jos/sasp/Poisson Window.html
1
(b) (3 pts) Write a matlab function that takes M as an argument and returns the Hann window as defined by the expression in brackets within the definition of the Hann-Poisson window2 . (This means that an initial and final zero are included in the returned window, so that there will be only M-2 nonzero samples. It is more commonly the convention that only nonzero samples “count” as part of the window, so this could also be called a length M − 2 window.) (c) (3 pts) Write a matlab function that takes M and α as arguments and returns the Hann-Poisson window3 . (d) (3 pts) Overlay plots of the magnitude window transforms for Hann, Poisson, and Hann-Poisson windows, where M = 21 in all cases, α = 2.0, and the zeropadded length is N = 8192 samples. Use the zero-phase zero-padding function from your previous homework to do the zero padding. Label the frequency axis in normalized-radian-frequency units (radians per sample), and normalize the amplitude axis in dB to a peak level of 0 dB. (e) (3 pts) repeat the previous item with α = 4.0. 3. (6 pts) Sketch the window w, corresponding to the window transform ∆
∆
W (ω) = asincM (ω) =
sin(M ω/2) sin(ω/2)
where M is an even integer. (Hint: This is not as trivial as it may first appear.) 4. (5 pts) Find out where the Chebyshev window “breaks down” in Matlab. Let the length be fixed at M = 31, and try various ripple specifications until there is an obvious error in the window obtained. Describe the source of the failure when the ripple specification is (a) too large and (b) too low. The following Matlab code can be used as a starting point: N = 8192; M = 31; w = chebwin(M,rip); W = fft(w,N); % normalize and clip the window transform in dB: Wdb = 20*log10(max(abs(W)/max(abs(W)),... 10^(-rip*1.5/20))); f = [0:N-1]/N - 0.5; plot(f,fftshift(Wdb)); grid; axis([-0.5 0.5 -rip*1.5 0]); hold on; plot([f(1),f(N)],[-rip,-rip],’--k’); text(f(1)+0.02,-rip+rip*1.5/20,... 2 3
http://ccrma.stanford.edu/˜jos/sasp/Hann Poisson Window.html http://ccrma.stanford.edu/˜jos/sasp/Hann Poisson Window.html
2
sprintf(’rip = %0.1f dB’,rip)); xlabel(’Normalized Frequency (cycles/sample)’); ylabel(’Magnitude (dB)’); title(sprintf(’Length %d Chebyshev window’,M)); 5. (5 pts) Sketch the window w(n) corresponding to the window transform 1 1 ∆ W (ω) = − asincM (ω + ω1 ) + asincM (ω) − asincM (ω − ω1 ) 2 2 where ω1 = 2π/M and M is assumed to be odd. 6. Suppose we are going to analyze a sinusoid whose frequency is f0 = 440 Hz using a rectangular window. Assume the window length is M = 255 and the sampling rate is fs = 8192 Hz. (a) (2 pts) What is the relative error (i.e., ∆f /f0 = |f0 − fˆ|/f0 ) between the actual frequency (f0 ) of the sinusoid and the peak frequency (fˆ) when there is no zeropadding? (b) (2 pts) Repeat part (a) with the zero-padding factor of 5. (c) (3 pts) Repeat part (a) with the zero-padding factor of 5 and the parabolic interpolation using the peak and its two neighbors. You don’t need to build algorithms for peak finding or parabolic interpolation for now. Just find the peak and its neighbors graphically.
3