[music-dsp] oversampled Fourier Transform
given a N-point window h(n), how do I create an oversampled Fourier transform H out of h(n)? I need to have very high resolution of the H so that I can find out the magnitude of arbitrary frequency with little frequency error. This is what I did: to oversample h(n) with factor 2, I first sampled the window h with N * 2 points, and then simply FFT the 'new h' with N*2-point. But the H with N*2 points looks very different from H with N points. What am I doing wrong? Here is detailed description with some sample output of the original and oversampled fft This is my windowing function: *function blackman(N)* a0 = 0.42659; a1 = 0.49656; a2 = 0.076849; w = zeros(N); for n = 1:N w[n] = a0 - a1 * cos(2*pi*n/(N-1)) + a2 * cos(4*pi*n/(N-1)); end return w; end I assume an oversampled fft of the window function by a factor of 2 would be *fft(blackman(N*2))* However, the fft values of the original and the oversampled value aren't increasing/decreasing in a similar way (I expect the oversampled version to takes twice as much samples to increase and decrease, in case the over-sample factor is 2). Here are some the fft sample values from 1 - 10. the original is FFT(blackman(512)), the oversampled is FFT(blackman(512 * 2)) k = 1 original: 217.9943833003818 + 0.0im oversampled: 873.2366098911153 + 0.0im k = 2 original: -127.11893711697411 - 2.340240508514901im oversampled: -508.4935935475258 - 2.3400747304536136im k = 3 original: 19.845020222085132 + 0.7309352010595286im oversampled: 78.87486492124337 + 0.7259756614936985im k = 4 original: 0.00044320678301414774 + 2.423663020980352e-5im oversampled: 0.0005516543659794568 + 7.61242746175689e-6im k = 5 original: 0.007485288714345478 + 0.0005517961643643162im oversampled: 0.00748836292444639 + 0.00013785416101184825im k = 6 original: 0.0060558634779024605 + 0.000558513387226844im oversampled: 0.00605535370050697 + 0.00013934970191932598im k = 7 original: 0.004581346143077173 + 0.0005075321354533686im oversampled: 0.004584081103264545 + 0.0001265980340830028im k = 8 original: 0.0035101418651299794 + 0.0004541979779443242im oversampled: 0.003516027383715841 + 0.0001132935089731223im k = 9 original: 0.0027517589243924367 + 0.00040747551324277084im oversampled: 0.0027600729094881945 + 0.00010164856218561597im k = 10 original: 0.002205670972923391 + 0.0003679932983093956im oversampled: 0.0022158009449307643 + 9.181309340287595e-5im What's wrong with my approach? Thanks in advance! -- dupswapdrop -- the music-dsp mailing list and website: subscription info, FAQ, source code archive, list archive, book reviews, dsp links http://music.columbia.edu/cmc/music-dsp http://music.columbia.edu/mailman/listinfo/music-dsp
Re: [music-dsp] oversampled Fourier Transform
If you just want higher frequency resolution, you don't need to do any oversampling in the time domain. Just zero-pad the time domain signal out to whatever long length you want, and then use an FFT of that size. E On Tue, Mar 31, 2015 at 3:12 PM, MF ukel...@gmail.com wrote: given a N-point window h(n), how do I create an oversampled Fourier transform H out of h(n)? I need to have very high resolution of the H so that I can find out the magnitude of arbitrary frequency with little frequency error. This is what I did: to oversample h(n) with factor 2, I first sampled the window h with N * 2 points, and then simply FFT the 'new h' with N*2-point. But the H with N*2 points looks very different from H with N points. What am I doing wrong? Here is detailed description with some sample output of the original and oversampled fft This is my windowing function: *function blackman(N)* a0 = 0.42659; a1 = 0.49656; a2 = 0.076849; w = zeros(N); for n = 1:N w[n] = a0 - a1 * cos(2*pi*n/(N-1)) + a2 * cos(4*pi*n/(N-1)); end return w; end I assume an oversampled fft of the window function by a factor of 2 would be *fft(blackman(N*2))* However, the fft values of the original and the oversampled value aren't increasing/decreasing in a similar way (I expect the oversampled version to takes twice as much samples to increase and decrease, in case the over-sample factor is 2). Here are some the fft sample values from 1 - 10. the original is FFT(blackman(512)), the oversampled is FFT(blackman(512 * 2)) k = 1 original: 217.9943833003818 + 0.0im oversampled: 873.2366098911153 + 0.0im k = 2 original: -127.11893711697411 - 2.340240508514901im oversampled: -508.4935935475258 - 2.3400747304536136im k = 3 original: 19.845020222085132 + 0.7309352010595286im oversampled: 78.87486492124337 + 0.7259756614936985im k = 4 original: 0.00044320678301414774 + 2.423663020980352e-5im oversampled: 0.0005516543659794568 + 7.61242746175689e-6im k = 5 original: 0.007485288714345478 + 0.0005517961643643162im oversampled: 0.00748836292444639 + 0.00013785416101184825im k = 6 original: 0.0060558634779024605 + 0.000558513387226844im oversampled: 0.00605535370050697 + 0.00013934970191932598im k = 7 original: 0.004581346143077173 + 0.0005075321354533686im oversampled: 0.004584081103264545 + 0.0001265980340830028im k = 8 original: 0.0035101418651299794 + 0.0004541979779443242im oversampled: 0.003516027383715841 + 0.0001132935089731223im k = 9 original: 0.0027517589243924367 + 0.00040747551324277084im oversampled: 0.0027600729094881945 + 0.00010164856218561597im k = 10 original: 0.002205670972923391 + 0.0003679932983093956im oversampled: 0.0022158009449307643 + 9.181309340287595e-5im What's wrong with my approach? Thanks in advance! -- dupswapdrop -- the music-dsp mailing list and website: subscription info, FAQ, source code archive, list archive, book reviews, dsp links http://music.columbia.edu/cmc/music-dsp http://music.columbia.edu/mailman/listinfo/music-dsp -- dupswapdrop -- the music-dsp mailing list and website: subscription info, FAQ, source code archive, list archive, book reviews, dsp links http://music.columbia.edu/cmc/music-dsp http://music.columbia.edu/mailman/listinfo/music-dsp
Re: [music-dsp] oversampled Fourier Transform
To expand on what Ethan wrote, it sounds like what you're trying to do is zero-pad the signal: http://www.dsprelated.com/dspbooks/mdft/Zero_Padding.html That said, whilst zero padding will give you an interpolated spectrum in the frequency domain, you may still miss the true location of your peaks, and it will also increase the computational cost. Another thing to look at is re-estimating the exact location of the peaks in your spectrum using parabolic interpolation: http://www.dsprelated.com/dspbooks/sasp/Quadratic_Interpolation_Spectral_Peaks.html Or, you could use the phase instead to compute the instantaneous frequency: http://en.wikipedia.org/wiki/Instantaneous_phase#Instantaneous_frequency Hope this helps, Best, Justin On Tue, Mar 31, 2015 at 6:31 PM, Ethan Duni ethan.d...@gmail.com wrote: If you just want higher frequency resolution, you don't need to do any oversampling in the time domain. Just zero-pad the time domain signal out to whatever long length you want, and then use an FFT of that size. E On Tue, Mar 31, 2015 at 3:12 PM, MF ukel...@gmail.com wrote: given a N-point window h(n), how do I create an oversampled Fourier transform H out of h(n)? I need to have very high resolution of the H so that I can find out the magnitude of arbitrary frequency with little frequency error. This is what I did: to oversample h(n) with factor 2, I first sampled the window h with N * 2 points, and then simply FFT the 'new h' with N*2-point. But the H with N*2 points looks very different from H with N points. What am I doing wrong? Here is detailed description with some sample output of the original and oversampled fft This is my windowing function: *function blackman(N)* a0 = 0.42659; a1 = 0.49656; a2 = 0.076849; w = zeros(N); for n = 1:N w[n] = a0 - a1 * cos(2*pi*n/(N-1)) + a2 * cos(4*pi*n/(N-1)); end return w; end I assume an oversampled fft of the window function by a factor of 2 would be *fft(blackman(N*2))* However, the fft values of the original and the oversampled value aren't increasing/decreasing in a similar way (I expect the oversampled version to takes twice as much samples to increase and decrease, in case the over-sample factor is 2). Here are some the fft sample values from 1 - 10. the original is FFT(blackman(512)), the oversampled is FFT(blackman(512 * 2)) k = 1 original: 217.9943833003818 + 0.0im oversampled: 873.2366098911153 + 0.0im k = 2 original: -127.11893711697411 - 2.340240508514901im oversampled: -508.4935935475258 - 2.3400747304536136im k = 3 original: 19.845020222085132 + 0.7309352010595286im oversampled: 78.87486492124337 + 0.7259756614936985im k = 4 original: 0.00044320678301414774 + 2.423663020980352e-5im oversampled: 0.0005516543659794568 + 7.61242746175689e-6im k = 5 original: 0.007485288714345478 + 0.0005517961643643162im oversampled: 0.00748836292444639 + 0.00013785416101184825im k = 6 original: 0.0060558634779024605 + 0.000558513387226844im oversampled: 0.00605535370050697 + 0.00013934970191932598im k = 7 original: 0.004581346143077173 + 0.0005075321354533686im oversampled: 0.004584081103264545 + 0.0001265980340830028im k = 8 original: 0.0035101418651299794 + 0.0004541979779443242im oversampled: 0.003516027383715841 + 0.0001132935089731223im k = 9 original: 0.0027517589243924367 + 0.00040747551324277084im oversampled: 0.0027600729094881945 + 0.00010164856218561597im k = 10 original: 0.002205670972923391 + 0.0003679932983093956im oversampled: 0.0022158009449307643 + 9.181309340287595e-5im What's wrong with my approach? Thanks in advance! -- dupswapdrop -- the music-dsp mailing list and website: subscription info, FAQ, source code archive, list archive, book reviews, dsp links http://music.columbia.edu/cmc/music-dsp http://music.columbia.edu/mailman/listinfo/music-dsp -- dupswapdrop -- the music-dsp mailing list and website: subscription info, FAQ, source code archive, list archive, book reviews, dsp links http://music.columbia.edu/cmc/music-dsp http://music.columbia.edu/mailman/listinfo/music-dsp -- Justin Salamon, PhD Postdoctoral researcher Music and Audio Research Laboratory (MARL) Center for Urban Science and Progress (CUSP) New York University, New York, NY www.justinsalamon.com -- dupswapdrop -- the music-dsp mailing list and website: subscription info, FAQ, source code archive, list archive, book reviews, dsp links http://music.columbia.edu/cmc/music-dsp http://music.columbia.edu/mailman/listinfo/music-dsp
Re: [music-dsp] oversampled Fourier Transform
On 3/31/15 6:53 PM, Justin Salamon wrote: To expand on what Ethan wrote, it sounds like what you're trying to do is zero-pad the signal: http://www.dsprelated.com/dspbooks/mdft/Zero_Padding.html That said, whilst zero padding will give you an interpolated spectrum in the frequency domain, you may still miss the true location of your peaks, how so? the only mechanism for missing the true location would be sidelobe interference from adjacent peaks, possibly including the peak(s) in negative frequency. and it will also increase the computational cost. Another thing to look at is re-estimating the exact location of the peaks in your spectrum using parabolic interpolation: http://www.dsprelated.com/dspbooks/sasp/Quadratic_Interpolation_Spectral_Peaks.html quadratic interpolation of peaks (given the three discrete samples around the peak) is a good ol' standby. i use it for autocorrelation to get the period to a fractional-sample precision. but i don't see why it would be more accurate than zero-padding before the FFT. Or, you could use the phase instead to compute the instantaneous frequency: http://en.wikipedia.org/wiki/Instantaneous_phase#Instantaneous_frequency needs a Hilbert transformer. and you would also need to isolate the sinusoidal partial from the other partials. i think so, anyway. -- r b-j r...@audioimagination.com Imagination is more important than knowledge. -- dupswapdrop -- the music-dsp mailing list and website: subscription info, FAQ, source code archive, list archive, book reviews, dsp links http://music.columbia.edu/cmc/music-dsp http://music.columbia.edu/mailman/listinfo/music-dsp
Re: [music-dsp] oversampled Fourier Transform
Hi MF, You might be misreading the frequency-domain versions of your windows. If you think about it, computing windows of different lengths (your oversampled windows) could be seen as obtaining windows of the same duration, but sampled /at different sample rates/. This will lead you to getting frequency-domain signals that are frequency-sampled at different frequency resolutions, and therefore your should be interpreting your frequency-domain samples differently. Let me clarify: Let's assume that you compute three different windows, w1(n), w2(n), and w3(n). All three windows last 1 second. w1(n) was computed by calling blackman(10). Equivalently, w2(n) was computed by calling blackman(20), and w3(n) by calling blackman(40). You can see this as having three windows at different sample rates of 10Hz, 20Hz, and 40Hz. Now if we take the FFT, we obtain W1[k], W2[k], and W3[k], with k being your frequency-domain sample (or bin). The number of frequency-domain samples (bins) of W1[k] is 10, while for W2[k] we have 20 bins, and for W3[k] we have 40 bins. From the sampling theory, we know that - the first 5 bins of W1[k] correspond to frequencies ranging from 0Hz to 5Hz (i.e., half your original sample rate); - the first 10 bins of W2[k] correspond to frequencies ranging from 0Hz to 10Hz (idem as above); - the first 20 bins of W3[k] correspond to frequencies ranging from 0Hz to 20Hz (idem as above). This means that the first 5 bins of all three W1[k], W2[k], and W3[k] will correspond to frequencies ranging from 0Hz to 5Hz. Bingo! Now you can see that all three time-domain windows, for which the only difference was the sample rate at which you obtained them, present /similar/ behavior//in the frequency domain. What we have learned here is that bin-to-bin magnitude (or phase) comparisons will only make sense if the bins correspond to the same frequency, and in this case they do : ) By oversampling the time-domain signal, the only thing you are doing is to add bins in the higher frequency region. In our example, W2[k] has more information than W1[k], but this new information only appears in the frequency region above 5Hz, which was /unknown/ to w1(n) in the first place. Now, going back to your objective. To obtain higher resolution in the lower frequency region, one could take W1[k] and directly resample or oversample it. This could be done by interpolating bins of W1[k]. Interpolation can be accomplished by convolving with a sinc function, which, in the time domain can be seen as multiplying by a rectangular function. Zero-padding is about changing the width of such rectangular function. As we all have done in the past : ), maybe you could read a bit on Windowing, Short-Time Fourier Transform, Zero-Padding, and Sinc Interpolation. It will all fall in place ; ) I totally recommend https://www.createspace.com/3177793 https://ccrma.stanford.edu/~jos/mdft/ , but there are many many other good materials out there! Cheers, Esteban On 3/31/2015 7:10 PM, robert bristow-johnson wrote: On 3/31/15 6:53 PM, Justin Salamon wrote: To expand on what Ethan wrote, it sounds like what you're trying to do is zero-pad the signal: http://www.dsprelated.com/dspbooks/mdft/Zero_Padding.html That said, whilst zero padding will give you an interpolated spectrum in the frequency domain, you may still miss the true location of your peaks, how so? the only mechanism for missing the true location would be sidelobe interference from adjacent peaks, possibly including the peak(s) in negative frequency. and it will also increase the computational cost. Another thing to look at is re-estimating the exact location of the peaks in your spectrum using parabolic interpolation: http://www.dsprelated.com/dspbooks/sasp/Quadratic_Interpolation_Spectral_Peaks.html quadratic interpolation of peaks (given the three discrete samples around the peak) is a good ol' standby. i use it for autocorrelation to get the period to a fractional-sample precision. but i don't see why it would be more accurate than zero-padding before the FFT. Or, you could use the phase instead to compute the instantaneous frequency: http://en.wikipedia.org/wiki/Instantaneous_phase#Instantaneous_frequency needs a Hilbert transformer. and you would also need to isolate the sinusoidal partial from the other partials. i think so, anyway. -- Esteban Maestre CIRMMT/CAML - McGill Univ MTG - Univ Pompeu Fabra http://ccrma.stanford.edu/~esteban -- dupswapdrop -- the music-dsp mailing list and website: subscription info, FAQ, source code archive, list archive, book reviews, dsp links http://music.columbia.edu/cmc/music-dsp http://music.columbia.edu/mailman/listinfo/music-dsp
