Hello, I am seeing a strange problem with the mechanism of zero-padding and performing FFT. The problem is bin center is shifting, instead of getting the bin center at 520 Hz. I am getting it at 448 Hz. I am attaching the code. It would be useful if someone can run the code, analyse the results and point out the mistake. The x-axis has been translated to frequency axis. Regards Bharat Pathak clear; close all; N1 = 16; % generate a sequence with 16 points N5 = 16*N1; % this will be the zero-padded seq len fin = 520; Fs = 8192; n1 = 0 : N1-1; n5 = 0 : N5-1; xin = sin(2*pi*fin*n1/Fs); x1 = xin; x5 = [xin zeros(1, N5-N1)]; h1 = fft(x1)/N1; h5 = fft(x5)/N5; figure; clf; stem(n1*(Fs/N1), abs(h1)); figure; clf; stem(n5*(Fs/N5), abs(h5));
zero padding of FFT
Started by ●February 17, 2008
Reply by ●February 17, 20082008-02-17
"bharat pathak" <bharat@arithos.com> wrote in message news:ppmdnR77HsOr8yXanZ2dnUVZ_vamnZ2d@giganews.com...> Hello, > > I am seeing a strange problem with the mechanism > of zero-padding and performing FFT. > > The problem is bin center is shifting, instead of > getting the bin center at 520 Hz. I am getting it > at 448 Hz. > > I am attaching the code. It would be useful if > someone can run the code, analyse the results and > point out the mistake. The x-axis has been translated > to frequency axis. > > Regards > Bharat Pathak > > clear; > close all; > > N1 = 16; % generate a sequence with 16 points > N5 = 16*N1; % this will be the zero-padded seq len > > fin = 520; > Fs = 8192; > n1 = 0 : N1-1; > n5 = 0 : N5-1; > > xin = sin(2*pi*fin*n1/Fs); > > x1 = xin; > x5 = [xin zeros(1, N5-N1)]; > > h1 = fft(x1)/N1; > h5 = fft(x5)/N5; > > figure; clf; stem(n1*(Fs/N1), abs(h1)); > figure; clf; stem(n5*(Fs/N5), abs(h5)); >At least it would be good to have the comments in the code make sense. "N1 = 16; % generate a sequence with 16 points" does no such thing. It sets N1 to the value 16 and that's all. No sequence. It may sound picky but that's the essence of coding effectively - being very picky. I can't make any sense out of your question because you haven't defined "bin center" and 520Hz makes no sense if one uses conjecture to guess what it means. I don't see an obvious problem here. The frequency samples of h1 are at intervals of 512Hz and the frequency intervals of h5 are 32Hz. That seems right if you have Fs=8,192. Thus the bin centers would generally be considered to be integer multiples of these numbers. A lot depends on what you're wanting to do. Fred
Reply by ●February 17, 20082008-02-17
On Feb 17, 9:26 am, "bharat pathak" <bha...@arithos.com> wrote:> Hello, > > I am seeing a strange problem with the mechanism > of zero-padding and performing FFT. > > The problem is bin center is shifting, instead of > getting the bin center at 520 Hz. I am getting it > at 448 Hz. > > I am attaching the code. It would be useful if > someone can run the code, analyse the results and > point out the mistake. The x-axis has been translated > to frequency axis. > > Regards > Bharat Pathak > > clear; > close all; > > N1 = 16; % generate a sequence with 16 points > N5 = 16*N1; % this will be the zero-padded seq len > > fin = 520; > Fs = 8192; > n1 = 0 : N1-1; > n5 = 0 : N5-1; > > xin = sin(2*pi*fin*n1/Fs); > > x1 = xin; > x5 = [xin zeros(1, N5-N1)]; > > h1 = fft(x1)/N1; > h5 = fft(x5)/N5; > > figure; clf; stem(n1*(Fs/N1), abs(h1)); > figure; clf; stem(n5*(Fs/N5), abs(h5));What you are seeing is the result of a combination of things, including sampling the spectrum at bin centers defined by your DFT length (multiples of 32 Hz in your example), so-called "spectral leakage" of non-bin-centered frequency content (e.g. 520 is not a multiple of 32) in finite length or windowed sample sets, plus leakage from the negative frequency image of the spectrum near DC, which tends to make low frequency magnitude peaks appear in even lower bins. The magnitude spectrum of a non-bin-centered sinusoid is not a single spike, but a hump whose shape depends depends on the transfrom of the window. The closer two non-bin-centered spectral frequencies are, the more the "humps" can interfere with each other in a spectral plot. The standard DFT procedure produces a two-sided spectrum from real input, splitting the magnitude between positive and negative frequencies (plot from -n5/2 to n5/2 to see this). Thus with spectral content near DC, the positive and negative frequency peaks may be close enough to interfere with each other, which can "pull" the apparent magnitude peaks lower in frequency. The only effect from zero padding will be to increase the frequency resolution of the DFT bins (e.g. the bin spacing will get closer in terms of frequency difference, which will sample the spectrum more finely). This will expose the shallow slope near any spectral peaks. IMHO. YMMV. -- rhn A.T nicholson d.0.t C-o-M http://www.nicholson.com/rhn/dsp.html
Reply by ●February 17, 20082008-02-17
On Feb 17, 12:41 pm, "Ron N." <rhnlo...@yahoo.com> wrote:> On Feb 17, 9:26 am, "bharat pathak" <bha...@arithos.com> wrote: > > > > > Hello, > > > I am seeing a strange problem with the mechanism > > of zero-padding and performing FFT. > > > The problem is bin center is shifting, instead of > > getting the bin center at 520 Hz. I am getting it > > at 448 Hz. > > > I am attaching the code. It would be useful if > > someone can run the code, analyse the results and > > point out the mistake. The x-axis has been translated > > to frequency axis. > > > Regards > > Bharat Pathak > > > clear; > > close all; > > > N1 = 16; % generate a sequence with 16 points > > N5 = 16*N1; % this will be the zero-padded seq len > > > fin = 520; > > Fs = 8192; > > n1 = 0 : N1-1; > > n5 = 0 : N5-1; > > > xin = sin(2*pi*fin*n1/Fs); > > > x1 = xin; > > x5 = [xin zeros(1, N5-N1)]; > > > h1 = fft(x1)/N1; > > h5 = fft(x5)/N5; > > > figure; clf; stem(n1*(Fs/N1), abs(h1)); > > figure; clf; stem(n5*(Fs/N5), abs(h5)); > > What you are seeing is the result of a combination of things, > including sampling the spectrum at bin centers defined by your > DFT length (multiples of 32 Hz in your example), so-called > "spectral leakage" of non-bin-centered frequency content (e.g. > 520 is not a multiple of 32) in finite length or windowed > sample sets, plus leakage from the negative frequency image > of the spectrum near DC, which tends to make low frequency > magnitude peaks appear in even lower bins. > > The magnitude spectrum of a non-bin-centered sinusoid is not > a single spike, but a hump whose shape depends depends on the > transfrom of the window. The closer two non-bin-centered > spectral frequencies are, the more the "humps" can interfere > with each other in a spectral plot. The standard DFT > procedure produces a two-sided spectrum from real input, > splitting the magnitude between positive and negative > frequencies (plot from -n5/2 to n5/2 to see this). Thus > with spectral content near DC, the positive and negative > frequency peaks may be close enough to interfere with each > other, which can "pull" the apparent magnitude peaks lower > in frequency.I'll add that a similar phenomena can occur for spectral content just below the Nyquist rate, where the frequency peaks can be "pulled upward". But this is seldom seen in practice because realistic anti-alias filters will remove most of this spectral content prior to sampling and applying the DFT/FFT to the data set.> The only effect from zero padding will be to increase the > frequency resolution of the DFT bins (e.g. the bin spacing > will get closer in terms of frequency difference, which will > sample the spectrum more finely). This will expose the > shallow slope near any spectral peaks. > > IMHO. YMMV. > -- > rhn A.T nicholson d.0.t C-o-M > http://www.nicholson.com/rhn/dsp.html
Reply by ●February 17, 20082008-02-17
Thanks Ron,
Your explaination seems to be most appropriate. By the way
any suggestions on how to model this and prove that the
interference phenomena, you explained, is the exact cause
of the frequency center drift?
Also is there any bare minimum number of samples required
to identify a discrete signal's frequency correctly for a
given Fs value?? Or we have to live with the frequency drifts.
Interestingly this phenomena is missed in Rick's book. Or
I couldn't locate it. Any comments Rick?
Regards
Bharat Pathak
Arithos Designs
www.arithos.com
P.S: Sorry for my code commenting style. When someone points to the moon,
stop looking at his finger tip. Look at the moon instead.
Reply by ●February 17, 20082008-02-17
bharat pathak wrote:> Thanks Ron, > > Your explaination seems to be most appropriate. By the way > any suggestions on how to model this and prove that the > interference phenomena, you explained, is the exact cause > of the frequency center drift? > > Also is there any bare minimum number of samples required > to identify a discrete signal's frequency correctly for a > given Fs value?? Or we have to live with the frequency drifts. > > Interestingly this phenomena is missed in Rick's book. Or > I couldn't locate it. Any comments Rick? > > Regards > Bharat Pathak > > Arithos Designs > www.arithos.com > > P.S: Sorry for my code commenting style. When someone points to the moon, > stop looking at his finger tip. Look at the moon instead.When someone says one thing known to be inaccurate, it's hard to credit the rest of his remarks. If you want a screw, don't ask for a nail. A better comment would be "Length of generated sequence". You got the next comment right. Jerry -- Engineering is the art of making what you want from things you can get. �����������������������������������������������������������������������
Reply by ●February 17, 20082008-02-17
On Feb 17, 6:16 pm, "bharat pathak" <bha...@arithos.com> wrote:> Thanks Ron, > > Your explaination seems to be most appropriate. By the way > any suggestions on how to model this and prove that the > interference phenomena, you explained, is the exact cause > of the frequency center drift?Just plot the spectra on *both* sides of zero, for the FT of rectangularly windowed low frequency sinusoids. Your own example shows it fairly clearly if you plot from -N5/2 to N5/2. As the frequency (fin) gets lower, you will see the two humps interfere with each other's shape, and thus shift their peak locations away from the known frequency of the sinusoidal input.> Also is there any bare minimum number of samples required > to identify a discrete signal's frequency correctly for a > given Fs value?? Or we have to live with the frequency drifts.Accuracy of frequency estimation depends more on the signal-to-noise ratio than the number of samples (as has been said before: in the zero noise case, only three non-aliased samples should be needed to solve for the three unknowns). In a high signal-to-noise situation, you should be able to compensate for low-frequency spectral self-interference. I don't know of a closed-form solution for this kind of estimation, but it's fairly easy to get decent sub-bin accuracy using a pre-computed table look-up for low-numbered DFT bins if the S/N ratio is good enough. I don't know if these minor details of ad-hoc frequency estimation are worth a DSP tips-and-tricks article... IMHO. YMMV. -- rhn A.T nicholson d.0.t C-o-M http://www.nicholson.com/rhn/dsp.html
Reply by ●February 18, 20082008-02-18
Fred Marshall wrote: "bharat pathak" <bharat@arithos.com> wrote (snip)>>N1 = 16; % generate a sequence with 16 points >>N5 = 16*N1; % this will be the zero-padded seq len(snip)> At least it would be good to have the comments in the code make sense.> "N1 = 16; % generate a sequence with 16 points" > does no such thing. It sets N1 to the value 16 and that's all.> No sequence.> It may sound picky but that's the essence of coding > effectively - being very picky.I have been in discussions on comment style before: Consider: n1=16; % set n1 to 16. The comment exactly describes the statement, but doesn't tell the reader what is important: Why is nl set to 16? Consider; nl=32; % generate a sequence with 16 points Now the reader knows it is related to the sequence length, but why does the value disagree with the comment? nl=16; % specify the sequence length. Now changes won't make the comment wrong, it doesn't claim something that it doesn't do, and does explain the reason for the statement. -- glen
Reply by ●February 18, 20082008-02-18
On 17 Feb, 18:26, "bharat pathak" <bha...@arithos.com> wrote:> Hello, > > � �I am seeing a strange problem with the mechanism > � �of zero-padding and performing FFT. > > � �The problem is bin center is shifting, instead of > � �getting the bin center at 520 Hz. I am getting it > � �at 448 Hz. > > � �I am attaching the code. It would be useful if > � �someone can run the code, analyse the results and > � �point out the mistake. The x-axis has been translated > � �to frequency axis.I ran a slightly modified version of your code (below) where I took care of some scaling factors and plotted the spectra in the same plot for easier comparisions. It seems you have dethroned a piece of 'common DSP knowledge.' Lots of people take for granted that the maximum of a zero-padded spectrum will coincide with the 'true' single-tone sinusoidal in the noiseless case. Your demo clearly shows that this is not ruse. Even with the fin = 512 case the zero-padded spectrum misses the 'true' frequency. The only place I can remember having seen zero-padding 'officially' mentioned as being able to help with frequency estimation is Kay: "Modern Spectrum Estimation" Prentice Hhall, 1988, section 13.3.1. On page 410, just below eqn 13.8 he mentions that data should be zero padded prior to a frequency estimation step such that 'a large number of frequency samples are computed.' Kay in turn refers to two previous papers, one of which is Palmer: "Coarse Frequency Estimation Using the DFT" IEEE Trans, Inf. Theory, 1974. I haven't looked at that paper but the appearance of the word 'coarse' is a bit interesting. All in all, I can't remember ever having seen a *proof* that zero-padding should produce a maximum at the 'true' frequency; I have only seen heuristic arguments. Your demo does what is required by zero padding an N-length sequence by M*N, which is that the DFT of the zero-padded sequence contains the original DFT points. Thanks for an eye-opening demo. Rune %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear; close all; N1 = 16; % generate a sequence with 16 points N5 = 16*N1; % this will be the zero-padded seq len %fin = 520; fin = 512; Fs = 8192; n1 = 0 : N1-1; n5 = 0 : N5-1; xin = sin(2*pi*fin*n1/Fs); x1 = xin; x5 = [xin zeros(1, N5-N1)]; h1 = fft(x1)/sqrt(N1); h5 = fft(x5)/sqrt(N1); clf plot(n1*(Fs/N1), abs(h1),'rx',n5*(Fs/N5), abs(h5),'ob');
Reply by ●February 18, 20082008-02-18
