FPGA Central

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));

"bharat pathak" <[email protected]> wrote in message
news[email protected] ...
> 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

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

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

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.

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.
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

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

Fred Marshall wrote:
"bharat pathak" <[email protected]> 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

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');

On 18 Feb, 10:16, Rune Allnor <all...@tele.ntnu.no> wrote:

Typo correction:

> Your demo clearly shows that this is not ruse.

Your demo clearly shows that this is not true.

Rune

Thanks Rune,

You made my day. One beer is sure for you and Ron
when we meet.

I was not sure if I was making any mistake,
in understanding, and hence posted it to the
forum.

Ron's explanation was good. I tried his idea's
in my demo and saw the drift even when I move closer
to Fs/2, and no drift when fin was closer to Fs/4.

Not sure if it is possible to derive a closed form
equation to tell the drift, which will be function
of N and window shape. Maybe this will make a nice
project by itself.

I feel this topic can be studied further, and may
be a nice article can come in DSP Tips and Tricks
column, or even forthcoming books of Rick and Steve.

Regards
Bharat Pathak

Arithos Designs
www.Arithos.com

On Mon, 18 Feb 2008 01:16:16 -0800 (PST), Rune Allnor
<[email protected]> wrote:

(snipped by Lyons)
>
>Thanks for an eye-opening demo.
>
>Rune

Hi Rune,
Yes, this is a very good example of how
leakage crosses over the 'zero-Hz'
boundary, and "pulls" the zero-padded spectral
peak toward zero Hz. (Ron. N mentioned this
behavior in his post.)

It's instructional to replace the "xin"
with a positive-frequency only (analytic) sinusoid
as:

xin = exp(j*2*pi*fin*n1/Fs);

In this case no negative-freq spectral component
leakage exists to distort the positive-freq spectral
component. And the maximum of a zero-padded spectrum
DOES coincide with the 'true' single-tone sinusoidal
in the noiseless case.

See Ya',
[-Rick-]

On Sun, 17 Feb 2008 20:16:55 -0600, "bharat pathak"
<[email protected]> 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

Hi Bharat,
Well, I did point out, on pages 72-73, how spectral
leakage will "cross over" the 'zero-hz' and 'Fs/2'
boundaries to make spectral 'humps' look asymmetrical.
I also mentioned that this asymmetry will be minimized
when the input sinusoid's frequency is close to Fs/4.

See Ya',
[-Rick-]

On Feb 18, 1:16 am, Rune Allnor <all...@tele.ntnu.no> wrote:

> ...
> 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.
>
> ...
>
> Rune
> ...

In the paper, Palmer explains that:
"In practice, this coarse approximation could form the first step in a
two-step acquisition procedure where the ultimate goal is to estimate
both phase and frequency."

The paper is an illustration of the procedure for the coarse component
of the procedure and an analysis of the errors produced by performing
only the coarse search.

For the performance of the two part search:
Single tone parameter estimation from discrete-time observations
Rife, D. Boorstyn, R.
Information Theory, IEEE Transactions on
Sep 1974 page(s): 591- 598 Volume: 20, Issue: 5
gives a formulation of a maximum likelihood estimator based on a
"coarse search" like Palmer's followed by a "fine search" and
discusses strategies for and errors from fine search implementations.

Zero-padding is a possible approach to the "fine search" component of
the parameter estimation from discrete observations.

These authors recognized the need to deal with finite sets of discrete
observations and the fact that this necessitated windowing. Palmer
compares 3 windows.

Dale B. Dalrymple
http://dbdimages.com

"bharat pathak" <[email protected]> wrote in message
news[email protected] ...

Bharat,

I'm sorry if my response was too terse. I surely missed the point.
I guess because "bin" to me always means the analysis cell and would have
had to do with sample spacing in frequency.

However, it now seems clear that you meant "the frequency of maximum
magnitude in the FFT" when you referred to "bin center". And now I see the
issue clearly.

I believe the "pulling" happens because we are dealing with a summation of
sincs or Dirichlets in frequency. As their peaks get closer and closer to
zero, their tails add differently until as they get very close to zero, they
overlap in the main lobes. 0

If you look at a sum of two Dirichlets that are equispaced and located at Fo
and fs-Fo, the "tails" will generally add diffferently from zero to Fo
(accompanied with fs-Fo to fs) than from from Fo to fs-Fo.

This is the same as saying "it's the leakage" but goes on to describe why.
I haven't tried it but suspect it will be both fractional frequency and
window length dependent as that's what affects how the tails add.

Fred

On Feb 18, 8:29 am, dbd <d...@ieee.org> wrote:
> On Feb 18, 1:16 am, Rune Allnor <all...@tele.ntnu.no> wrote:
>
>
>
> > ...
> > 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.
>
> > ...
>
> > Rune
> > ...
>
> In the paper, Palmer explains that:
> "In practice, this coarse approximation could form the first step in a
> two-step acquisition procedure where the ultimate goal is to estimate
> both phase and frequency."
>
> The paper is an illustration of the procedure for the coarse component
> of the procedure and an analysis of the errors produced by performing
> only the coarse search.
>
> For the performance of the two part search:
> Single tone parameter estimation from discrete-time observations
> Rife, D. Boorstyn, R.
> Information Theory, IEEE Transactions on
> Sep 1974 page(s): 591- 598 Volume: 20, Issue: 5
> gives a formulation of a maximum likelihood estimator based on a
> "coarse search" like Palmer's followed by a "fine search" and
> discusses strategies for and errors from fine search implementations.
>
> Zero-padding is a possible approach to the "fine search" component of
> the parameter estimation from discrete observations.
>
> These authors recognized the need to deal with finite sets of discrete
> observations and the fact that this necessitated windowing. Palmer
> compares 3 windows.

Interesting... Windowing helps diminish interfering spectral
interference and noise, but actually seems counter-productive
for spectral estimation after zero-padding. Zero-padding
brings out the shape of the main Sinc lobe, which seems to
make a "fine search" (by autocorrelation methods, or some
such) more accurate (in the case of isolated spectral
content). Whereas windowing "smears out" the main Sinc lobe,
as well as reduces the total information content available
from the zero-padded DFT (in exchange for lowering the level
of adjacent spectral interference of course).

> Dale B. Dalrymplehttp://dbdimages.com


IMHO. YMMV.
--
rhn A.T nicholson d.0.t C-o-M
http://www.nicholson.com/rhn/dsp.html

On Feb 18, 9:16 am, "Fred Marshall" <fmarshallx@remove_the_x.acm.org>
wrote:
> "bharat pathak" <bha...@arithos.com> wrote in message
>
> news[email protected] ...
>
> Bharat,
>
> I'm sorry if my response was too terse. I surely missed the point.
> I guess because "bin" to me always means the analysis cell and would have
> had to do with sample spacing in frequency.
>
> However, it now seems clear that you meant "the frequency of maximum
> magnitude in the FFT" when you referred to "bin center". And now I see the
> issue clearly.
>
> I believe the "pulling" happens because we are dealing with a summation of
> sincs or Dirichlets in frequency. As their peaks get closer and closer to
> zero, their tails add differently until as they get very close to zero, they
> overlap in the main lobes. 0
>
> If you look at a sum of two Dirichlets that are equispaced and located at Fo
> and fs-Fo, the "tails" will generally add diffferently from zero to Fo
> (accompanied with fs-Fo to fs) than from from Fo to fs-Fo.
>
> This is the same as saying "it's the leakage" but goes on to describe why.
> I haven't tried it but suspect it will be both fractional frequency and
> window length dependent as that's what affects how the tails add.

It's dependent on both of these (the reciprocal of the window
length vs. the frequency offset from DC, or from Fs/2), and on
the shape of the window (rectangular vs. something else), *and*
on whether the sinusoid being estimated is more odd or even
in the DFT/FFT aperture (since that also affects how the two
complex-conjugate waveforms interfere).

I have a formula which includes some of these factors somewhere
on my dsp blog/web page:

http://www.nicholson.com/rhn/dsp.html

One promising method for assisting in the frequency estimation
of zero-padded FFTs near the low bins is to pre-compute (for
the given window) look-up tables with compensation factors for
the amount of "pull" vs. frequency.


IMHO. YMMV.
--
rhn A.T nicholson d.0.t C-o-M

On Feb 18, 2:01 pm, "Ron N." <rhnlo...@yahoo.com> wrote:
> On Feb 18, 8:29 am, dbd <d...@ieee.org> wrote:
>
>
>
> > On Feb 18, 1:16 am, Rune Allnor <all...@tele.ntnu.no> wrote:
>
> > > ...
> > > 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.
>
> > > ...
>
> > > Rune
> > > ...
>
> > In the paper, Palmer explains that:
> > "In practice, this coarse approximation could form the first step in a
> > two-step acquisition procedure where the ultimate goal is to estimate
> > both phase and frequency."
>
> > The paper is an illustration of the procedure for the coarse component
> > of the procedure and an analysis of the errors produced by performing
> > only the coarse search.
>
> > For the performance of the two part search:
> > Single tone parameter estimation from discrete-time observations
> > Rife, D. Boorstyn, R.
> > Information Theory, IEEE Transactions on
> > Sep 1974 page(s): 591- 598 Volume: 20, Issue: 5
> > gives a formulation of a maximum likelihood estimator based on a
> > "coarse search" like Palmer's followed by a "fine search" and
> > discusses strategies for and errors from fine search implementations.
>
> > Zero-padding is a possible approach to the "fine search" component of
> > the parameter estimation from discrete observations.
>
> > These authors recognized the need to deal with finite sets of discrete
> > observations and the fact that this necessitated windowing. Palmer
> > compares 3 windows.
>
> Interesting... Windowing helps diminish interfering spectral
> interference and noise, but actually seems counter-productive
> for spectral estimation after zero-padding. Zero-padding
> brings out the shape of the main Sinc lobe, which seems to
> make a "fine search" (by autocorrelation methods, or some
> such) more accurate (in the case of isolated spectral
> content). Whereas windowing "smears out" the main Sinc lobe,
> as well as reduces the total information content available
> from the zero-padded DFT (in exchange for lowering the level
> of adjacent spectral interference of course).

Sorry, I meant "non-rectangular windowing". Perhaps the
paper compared rectangular windows vs. others.

> > Dale B. Dalrymplehttp://dbdimages.com
>
> IMHO. YMMV.
> --
> rhn A.T nicholson d.0.t C-o-M
> http://www.nicholson.com/rhn/dsp.html

On Feb 18, 2:01 pm, "Ron N." <rhnlo...@yahoo.com> wrote:
> On Feb 18, 8:29 am, dbd <d...@ieee.org> wrote:
>
....
> > These authors recognized the need to deal with finite sets of discrete
> > observations and the fact that this necessitated windowing. Palmer
> > compares 3 windows.
>
> Interesting... Windowing helps diminish interfering spectral
> interference and noise, but actually seems counter-productive
> for spectral estimation after zero-padding. Zero-padding
> brings out the shape of the main Sinc lobe, which seems to
> make a "fine search" (by autocorrelation methods, or some
> such) more accurate (in the case of isolated spectral
> content). Whereas windowing "smears out" the main Sinc lobe,
> as well as reduces the total information content available
> from the zero-padded DFT (in exchange for lowering the level
> of adjacent spectral interference of course).
>
> > Dale B. Dalrymplehttp://dbdimages.com
>
> IMHO. YMMV.
> --
> rhn A.T nicholson d.0.t C-o-M
> http://www.nicholson.com/rhn/dsp.html

Of course our mileage may vary, we may have interfering tones (near
tones of similar strength or distant tones of greater strength),
colored noise or minimum bias concerns that can dictate desirable
windows. Zero padding works to decrease scalloping loss error as does
windowing. The history of frequency determination from sets of
discrete values contains a number of windows designed specifically to
determine frequency. We can't avoid windowing, it has uses and if we
learn how to use it, it can vary our mileage, it is a useful tool for
the kit.

Dale B. Dalrymple

On Feb 18, 1:16 am, Rune Allnor <all...@tele.ntnu.no> wrote:
> On 17 Feb, 18:26, "bharat pathak" <bha...@arithos.com> wrote:
....
> 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;

Since zero-padding gives the same result as Sinc
interpolation, in the zero noise case it should still
converge on the input frequency peak near Fs/4, where the
effects from the image spectra balance out. Away from Fs/4
it should still converge of the sum on the two frequency
spectra. Since this sum is deterministic, the input
frequency can still be reverse engineered.




IMHO. YMMV.
--
rhn A.T nicholson d.0.t C-o-M

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));

If you want to see another interesting phenomena
with the zero-padded FFT, try replacing

> xin = sin(2*pi*fin*n1/Fs);

with:

xin = cos(2*pi*fin*n1/Fs);

keeping the same frequency, and see what happens
to the frequency peak in the upsampled DFT result.




YMMV.
--
rhn A.T nicholson d.0.t C-o-M

On Mon, 18 Feb 2008 08:16:43 -0600, "bharat pathak"
<[email protected]> wrote:

>Thanks Rune,
>
> You made my day. One beer is sure for you and Ron
> when we meet.
>
> I was not sure if I was making any mistake,
> in understanding, and hence posted it to the
> forum.
>
> Ron's explanation was good. I tried his idea's
> in my demo and saw the drift even when I move closer
> to Fs/2, and no drift when fin was closer to Fs/4.
>
> Not sure if it is possible to derive a closed form
> equation to tell the drift, which will be function
> of N and window shape. Maybe this will make a nice
> project by itself.

(snipped by Lyons)

Hi Bharat,
I imagine that it wouldn't be too difficult to
come up with an equation for the "drifted" spectrum, and
I'm thinking it would be a function of the DFT input
sinusoid's frequency, DFT-size (N), and the DFT of any
window that's used in the time domain. But I bet
such an equation would be moderately complicated
(fairly messy).

As for an equation that predicts the "drift in frequency"
caused by leakage wrap-around (wrapping around zero Hz &
Fs/2 Hz), now that sounds pretty difficult to me.
Wouldn't we have to take the derivative of the "drifted"
spectrum equation, and set that derivative equal to zero
in order to find the peak of the "drifted" spectrum?

Humm, ...sounds like a homework problem from the Oppenheim
and Schafer DSP book. :-)

[-Rick-]

"Ron N." <[email protected]> wrote in message
news:[email protected]...
> On Feb 18, 1:16 am, Rune Allnor <all...@tele.ntnu.no> wrote:
>> On 17 Feb, 18:26, "bharat pathak" <bha...@arithos.com> wrote:
> ...
>> 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;
>
> Since zero-padding gives the same result as Sinc
> interpolation, in the zero noise case it should still
> converge on the input frequency peak near Fs/4, where the
> effects from the image spectra balance out. Away from Fs/4
> it should still converge of the sum on the two frequency
> spectra. Since this sum is deterministic, the input
> frequency can still be reverse engineered.

That's a good observation Ron. So sinusoids near zero and fs/2 are affected
similarly but in the opposite direction due to the Dirichlet overlaps being
asymmetrical.

Fred

On Feb 19, 4:54 am, Rick Lyons <R.Lyons@_BOGUS_ieee.org> wrote:
> On Mon, 18 Feb 2008 08:16:43 -0600, "bharat pathak"
> <bha...@arithos.com> wrote:
> >Thanks Rune,
>
> > You made my day. One beer is sure for you and Ron
> > when we meet.
>
> > I was not sure if I was making any mistake,
> > in understanding, and hence posted it to the
> > forum.
>
> > Ron's explanation was good. I tried his idea's
> > in my demo and saw the drift even when I move closer
> > to Fs/2, and no drift when fin was closer to Fs/4.
>
> > Not sure if it is possible to derive a closed form
> > equation to tell the drift, which will be function
> > of N and window shape. Maybe this will make a nice
> > project by itself.
>
> (snipped by Lyons)
>
> Hi Bharat,
> I imagine that it wouldn't be too difficult to
> come up with an equation for the "drifted" spectrum, and
> I'm thinking it would be a function of the DFT input
> sinusoid's frequency, DFT-size (N), and the DFT of any
> window that's used in the time domain.

It also depends on the even/oddness of the sinusoid
being estimated in the DFT window. To see this yourself,
try the OP's example with a cosine instead of a sine test
waveform at the same frequency.

> As for an equation that predicts the "drift in frequency"
> caused by leakage wrap-around (wrapping around zero Hz &
> Fs/2 Hz), now that sounds pretty difficult to me.

I have an attempt at an describing the shape of this spectrum
somewhere on by dsp web page/blog:
http://www.nicholson.com/rhn/dsp.html

IMHO. YMMV.
--
rhn A.T nicholson d.0.t C-o-M

"Ron N." <[email protected]> wrote in message
news:[email protected]...
> 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));
>
> If you want to see another interesting phenomena
> with the zero-padded FFT, try replacing
>
>> xin = sin(2*pi*fin*n1/Fs);
>
> with:
>
> xin = cos(2*pi*fin*n1/Fs);
>
> keeping the same frequency, and see what happens
> to the frequency peak in the upsampled DFT result.
>

A gross generalization on this would be:
The sin is cut off at its peaks, causing "steps".
The cos is cut off at its zeros causing "steps" only in in the first
derivative.
Anything in between (i.e. different phases) would have a bit of both and, I
conjecture, that the leakage energy would fall in between as well.
The frequency content of the sin thus chopped is greater - so the tails must
have some degree of constructive interference.
The frequency content of the cos thus chopped is minimum (?) - so the tails
must have some degree of destructive interference.

An interesting student's exercise.

Then expand it to non-integer numbers of cycles in the window - i.e.
different frequencies for the same window length.

Related to this and Rick's comment about formulae:

Taking the thing in a continuous frequency context (sampled time -
periodic/continuous frequency) one can compute the Dirichlet pair sum and
find the peaks and plot their positions. It's not a "formula" but actually
I'll bet one can be found. Someone like Robert Israel at ubc.ca I'll bet
could do it. He's helped in the past with such "math" things.

Then, with this, no doubt one could find the peak sample in a discrete
frequency context. The caution would be that there's time domain aliasing
to deal with.

Of course, in the original case there's frequency domain aliasing that we
usually call leakage. It goes like this:
There are two cases:

Case I: Pick a sinusoid at some frequency. Decide on a reasonable sample
rate. Sample it. There is no aliasing. Now window it ... that's where the
aliasing comes in. Not obvious? Look at Case II.

Case II: Pick a sinusoid at some frequency - same as Case I. Decide on a
reasonable sample rate - same as Case I. Now window it - same as Case I.
Now sample it. Oops! The originally selected sample rate is now too low
because of the frequency content of the edges created by the windowing (i.e.
truncation)! So aliasing is more evident in this case.

We stipulate that the resulting samples are the same for Case I and Case II.
So, if there's aliasing with one, there's aliasing with the other. One
cannot tell the difference between the two as they have the same starting
point and the same ending result.

Fred