Hi,

I am trying to compute the amplitude and phase of signal with specific
frequency. The frequency of the signal is known (measured with laser - the
rotation machine). I tried to simulate the computation in matlab first:

close all
clear all
% clc

f = 50;                 % Frequency of signal
fs = 64000;             % Sample frequency
% fs = 2560;
phase_deg = 10;         % Phase of signal [degre]
t = 0.5;                % Time of signal

x = (0:1/(fs):t-(1/fs));            % Sampled time
phase_rad = 2*pi*phase_deg/360;     % Phase in radians
y = cos(2*pi*f*x + phase_rad);      % Generate sampled signal
figure;plot(x,y);                   % Plot signal

% y = y + (0.1*randn(size(y)));       % Add some noise

% Apply windowing function on the signal
w = 2.*hanning(length(y));
y = y.*w';
figure;plot(x,y);

% Compute FFT and plot abs and angle
Y = fft(y);
figure;plot(abs(Y));
figure;plot(angle(Y));

% Find the signal in frequency domain
[max_num,max_pos] = max(abs(Y));

% Value of signal (Amplitude)
value_x = abs(Y(max_pos))/(length(x)/2)

% Angle of signal [degre]
angle_x = 360*(angle(Y(max_pos)))/(2*pi)

I used the FFT, found the main peak (that represent the signal of
interest), and compute the amplitude and phase of this peak.

Everything works great, even a small noise don't change the value of
amplitude and phase reasonable.

However, the problem occurs when the FFT is computed from Non-Integer
number of signal periods. What is in praxis all the time. Try to change the
frequency of the signal to e.g. 50.4, 51.7, 53.9 etc. The classical
approach should be to use the windowing function. I used Hanning (i tried
lot of them) in example. Windowing will help with amplitude. But the phase
of the signal is computed totally incorrectly, even when the windowing is
used. I found that the windowing function do not play any (or not
reasonable) role in phase computation.

Is there any way how to compute the phase in praxis, when the FFT is
computed from Non-Integer number of signal periods?

Thanks in advance.
Afinko

On Jul 13, 5:35=A0am, "Afinko" <afinko@n_o_s_p_a_m.gmail.com> wrote:
> Hi,
>
> I am trying to compute the amplitude and phase of signal with specific
> frequency. The frequency of the signal is known (measured with laser - th=
e
> rotation machine). I tried to simulate the computation in matlab first:
>
> close all
> clear all
> % clc
>
> f =3D 50; =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 % Frequency of signal
> fs =3D 64000; =A0 =A0 =A0 =A0 =A0 =A0 % Sample frequency
> % fs =3D 2560;
> phase_deg =3D 10; =A0 =A0 =A0 =A0 % Phase of signal [degre]
> t =3D 0.5; =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0% Time of signal
>
> x =3D (0:1/(fs):t-(1/fs)); =A0 =A0 =A0 =A0 =A0 =A0% Sampled time
> phase_rad =3D 2*pi*phase_deg/360; =A0 =A0 % Phase in radians
> y =3D cos(2*pi*f*x + phase_rad); =A0 =A0 =A0% Generate sampled signal
> figure;plot(x,y); =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 % Plot signal
>
> % y =3D y + (0.1*randn(size(y))); =A0 =A0 =A0 % Add some noise
>
> % Apply windowing function on the signal
> w =3D 2.*hanning(length(y));
> y =3D y.*w';
> figure;plot(x,y);
>
> % Compute FFT and plot abs and angle
> Y =3D fft(y);
> figure;plot(abs(Y));
> figure;plot(angle(Y));
>
> % Find the signal in frequency domain
> [max_num,max_pos] =3D max(abs(Y));
>
> % Value of signal (Amplitude)
> value_x =3D abs(Y(max_pos))/(length(x)/2)
>
> % Angle of signal [degre]
> angle_x =3D 360*(angle(Y(max_pos)))/(2*pi)
>
> I used the FFT, found the main peak (that represent the signal of
> interest), and compute the amplitude and phase of this peak.
>
> Everything works great, even a small noise don't change the value of
> amplitude and phase reasonable.
>
> However, the problem occurs when the FFT is computed from Non-Integer
> number of signal periods. What is in praxis all the time. Try to change t=
he
> frequency of the signal to e.g. 50.4, 51.7, 53.9 etc. The classical
> approach should be to use the windowing function. I used Hanning (i tried
> lot of them) in example. Windowing will help with amplitude. But the phas=
e
> of the signal is computed totally incorrectly, even when the windowing is
> used. I found that the windowing function do not play any (or not
> reasonable) role in phase computation.
>
> Is there any way how to compute the phase in praxis, when the FFT is
> computed from Non-Integer number of signal periods?
>
> Thanks in advance.
> Afinko

If you're only concerned with amplitude/phase of a specific frequency,
you don't need to go to the trouble of calculating a full FFT. Google
"Goertzel algorithm." One thing to remember is that phase is relative
to the time at which you started sampling the sinusoid; unless your
blocks of input data to the FFT are synchronized in some way, the
sinusoid's initial phase is going to be arbitrary.

Jason

On Jul 13, 8:21=A0am, Jason <cincy...@gmail.com> wrote:
> On Jul 13, 5:35=A0am, "Afinko" <afinko@n_o_s_p_a_m.gmail.com> wrote:
>
....
>
> > Is there any way how to compute the phase in praxis, when the FFT is
> > computed from Non-Integer number of signal periods?
>
>
> If you're only concerned with amplitude/phase of a specific frequency,
> you don't need to go to the trouble of calculating a full FFT. Google
> "Goertzel algorithm."

if the exact frequency, f0, is known in advance, i would just
correlate the sampled signal to cos(2*pi*f0*n) and sin(2*pi*f0*n),
after LPFing both separately, apply the 4 quadrant arg{} for phase
(amplitude is from square rooting the sums of squares).

r b-j

On Jul 13, 2:35=A0am, "Afinko" <afinko@n_o_s_p_a_m.gmail.com> wrote:
> Hi,
>
> I am trying to compute the amplitude and phase of signal with specific
> frequency. The frequency of the signal is known (measured with laser - th=
e
> rotation machine). I tried to simulate the computation in matlab first:
>
> close all
> clear all
> % clc
>
> f =3D 50; =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 % Frequency of signal
> fs =3D 64000; =A0 =A0 =A0 =A0 =A0 =A0 % Sample frequency
> % fs =3D 2560;
> phase_deg =3D 10; =A0 =A0 =A0 =A0 % Phase of signal [degre]
> t =3D 0.5; =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0% Time of signal
>
> x =3D (0:1/(fs):t-(1/fs)); =A0 =A0 =A0 =A0 =A0 =A0% Sampled time
> phase_rad =3D 2*pi*phase_deg/360; =A0 =A0 % Phase in radians
> y =3D cos(2*pi*f*x + phase_rad); =A0 =A0 =A0% Generate sampled signal
> figure;plot(x,y); =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 % Plot signal
>
> % y =3D y + (0.1*randn(size(y))); =A0 =A0 =A0 % Add some noise
>
> % Apply windowing function on the signal
> w =3D 2.*hanning(length(y));
> y =3D y.*w';
> figure;plot(x,y);
>
> % Compute FFT and plot abs and angle
> Y =3D fft(y);
> figure;plot(abs(Y));
> figure;plot(angle(Y));
>
> % Find the signal in frequency domain
> [max_num,max_pos] =3D max(abs(Y));
>
> % Value of signal (Amplitude)
> value_x =3D abs(Y(max_pos))/(length(x)/2)
>
> % Angle of signal [degre]
> angle_x =3D 360*(angle(Y(max_pos)))/(2*pi)
>
> I used the FFT, found the main peak (that represent the signal of
> interest), and compute the amplitude and phase of this peak.
>
> Everything works great, even a small noise don't change the value of
> amplitude and phase reasonable.
>
> However, the problem occurs when the FFT is computed from Non-Integer
> number of signal periods. What is in praxis all the time. Try to change t=
he
> frequency of the signal to e.g. 50.4, 51.7, 53.9 etc. The classical
> approach should be to use the windowing function. I used Hanning (i tried
> lot of them) in example. Windowing will help with amplitude. But the phas=
e
> of the signal is computed totally incorrectly, even when the windowing is
> used. I found that the windowing function do not play any (or not
> reasonable) role in phase computation.
>
> Is there any way how to compute the phase in praxis, when the FFT is
> computed from Non-Integer number of signal periods?

If you reference the phase to the center of your FFT aperture
(by means of an pre-fftshift, or post flipping of alternate
bins), the evenness to oddness ratio of sinusoid can be
interpolated, even if it is not an integer number of signal
periods in the FFT aperture.  And the evenness to oddness
ratio is cleanly related to the atan2() phase.

The idiocy where people say you can't interpolate phase
is because they are using the start of the FFT aperture;
and a non-integer period sinusoid in discontinuous at the
boundaries of an FFT aperture, so trying to reference the
phase to a discontinuity is certainly non-pictorial and
non-intuitive.

Also, most windows go to zero at the boundaries of the
FFT aperature.  So you end up looking for the phase of
nothing.  The signal, after windowing, is in the middle.

So just move your phase reference to the center.



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

>If you reference the phase to the center of your FFT aperture
>(by means of an pre-fftshift, or post flipping of alternate
>bins), the evenness to oddness ratio of sinusoid can be
>interpolated, even if it is not an integer number of signal
>periods in the FFT aperture.  And the evenness to oddness
>ratio is cleanly related to the atan2() phase.
>
>The idiocy where people say you can't interpolate phase
>is because they are using the start of the FFT aperture;
>and a non-integer period sinusoid in discontinuous at the
>boundaries of an FFT aperture, so trying to reference the
>phase to a discontinuity is certainly non-pictorial and
>non-intuitive.
>
>Also, most windows go to zero at the boundaries of the
>FFT aperature.  So you end up looking for the phase of
>nothing.  The signal, after windowing, is in the middle.
>
>So just move your phase reference to the center.

Hi Ron N.,
Thanks for your reply.
I read your post several times, but I do not understand what exactly you
mean by "move your phase reference to the center".
Can you pleas post some MATLAB example how to do this?
Or can you modify my MATLAB code and show how does it work?

Thanks in advance.
Afinko

>If you're only concerned with amplitude/phase of a specific frequency,
>you don't need to go to the trouble of calculating a full FFT. Google
>"Goertzel algorithm." One thing to remember is that phase is relative
>to the time at which you started sampling the sinusoid; unless your
>blocks of input data to the FFT are synchronized in some way, the
>sinusoid's initial phase is going to be arbitrary.

Hi Jason,
Thanks for your reply.
It is a great idea, because I can always use so many inputs to the
"Goertzel algorithm" that I will always take the integer number of signal
periods. I used "goertzel" function in MATLAB. It works great on simulated
data. But when I used the measured data, the phase information is still
very unstable :(

The only way how I obtained gut results on measured data, that I found till
now, is this approach:
I filtered the input data with sharp bandpass IIR filter with center
frequency that I measured with laser. Then I used synchronized averaging
with reference from laser. I obtained one period of examined signal (it
looks very similar to the one sinus period). Then I simply find the maximum
of this signal and obtain a phase.
The only problem is, that this kind of IIR filter has very nonlinear phase
and difference between delay of 50 Hz and 50.2Hz is more than one period of
50 Hz signal :(

I am going to try filter the measured signal, as well as reference signal
from laser with the same IIR. In this way, the shift of examined frequency
should by the same for both signals. Hopefully it will work.

Any other ideas?
Thanks in advance.
Afinko

>if the exact frequency, f0, is known in advance, i would just
>correlate the sampled signal to cos(2*pi*f0*n) and sin(2*pi*f0*n),
>after LPFing both separately, apply the 4 quadrant arg{} for phase
>(amplitude is from square rooting the sums of squares).

hi robert bristow-johnson,
Thanks for your reply.
The problem is, that in measured signal, there is not only this one
examined frequency component, but there is a lot of another frequency
components, with much higher amplitudes. The best solution would be to
filter this signal first, but such a sharp bandpass filter is a problem to
achieve with FIR, therefore I used IIR, as I posted earlier. The next
problem is, that the measured frequency from laser is changing for all
machine. Therefore I need to compute the filter coefficients realtime in
DSP before filter. To compute the one biquad of IIR such as with
"iirpeak.m" function in MATLAB looks quit easy.

Afinko

On Jul 14, 5:16=A0am, "Afinko" <afinko@n_o_s_p_a_m.gmail.com> wrote:
> >If you reference the phase to the center of your FFT aperture
> >(by means of an pre-fftshift, or post flipping of alternate
> >bins), the evenness to oddness ratio of sinusoid can be
> >interpolated, even if it is not an integer number of signal
> >periods in the FFT aperture. =A0And the evenness to oddness
> >ratio is cleanly related to the atan2() phase.
>
> >The idiocy where people say you can't interpolate phase
> >is because they are using the start of the FFT aperture;
> >and a non-integer period sinusoid in discontinuous at the
> >boundaries of an FFT aperture, so trying to reference the
> >phase to a discontinuity is certainly non-pictorial and
> >non-intuitive.
>
> >Also, most windows go to zero at the boundaries of the
> >FFT aperature. =A0So you end up looking for the phase of
> >nothing. =A0The signal, after windowing, is in the middle.
>
> >So just move your phase reference to the center.
>
> Hi Ron N.,
> Thanks for your reply.
> I read your post several times, but I do not understand what exactly you
> mean by "move your phase reference to the center".
> Can you pleas post some MATLAB example how to do this?
> Or can you modify my MATLAB code and show how does it work?

For a phase reference at t =3D t0:

f =3D (Fs/N)*[0:N-1]';
X =3D exp(i*2*pi*f*t0).*fft(x);

Hope this helps.

Greg

On Jul 13, 5:35=A0am, "Afinko" <afinko@n_o_s_p_a_m.gmail.com> wrote:
> Hi,
>
> I am trying to compute the amplitude and phase of signal with specific
> frequency. The frequency of the signal is known (measured with laser - th=
e
> rotation machine). I tried to simulate the computation in matlab first:
>
> close all
> clear all
> % clc
>
> f =3D 50; =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 % Frequency of signal
> fs =3D 64000; =A0 =A0 =A0 =A0 =A0 =A0 % Sample frequency
> % fs =3D 2560;
> phase_deg =3D 10; =A0 =A0 =A0 =A0 % Phase of signal [degre]
> t =3D 0.5; =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0% Time of signal
>
> x =3D (0:1/(fs):t-(1/fs)); =A0 =A0 =A0 =A0 =A0 =A0% Sampled time
> phase_rad =3D 2*pi*phase_deg/360; =A0 =A0 % Phase in radians
> y =3D cos(2*pi*f*x + phase_rad); =A0 =A0 =A0% Generate sampled signal
> figure;plot(x,y); =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 % Plot signal
>
> % y =3D y + (0.1*randn(size(y))); =A0 =A0 =A0 % Add some noise
>
> % Apply windowing function on the signal
> w =3D 2.*hanning(length(y));
> y =3D y.*w';
> figure;plot(x,y);
>
> % Compute FFT and plot abs and angle
> Y =3D fft(y);
> figure;plot(abs(Y));
> figure;plot(angle(Y));
>
> % Find the signal in frequency domain
> [max_num,max_pos] =3D max(abs(Y));
>
> % Value of signal (Amplitude)
> value_x =3D abs(Y(max_pos))/(length(x)/2)
>
> % Angle of signal [degre]
> angle_x =3D 360*(angle(Y(max_pos)))/(2*pi)
>
> I used the FFT, found the main peak (that represent the signal of
> interest), and compute the amplitude and phase of this peak.
>
> Everything works great, even a small noise don't change the value of
> amplitude and phase reasonable.
>
> However, the problem occurs when the FFT is computed from Non-Integer
> number of signal periods. What is in praxis all the time. Try to change t=
he
> frequency of the signal to e.g. 50.4, 51.7, 53.9 etc. The classical
> approach should be to use the windowing function. I used Hanning (i tried
> lot of them) in example. Windowing will help with amplitude. But the phas=
e
> of the signal is computed totally incorrectly, even when the windowing is
> used. I found that the windowing function do not play any (or not
> reasonable) role in phase computation.
>
> Is there any way how to compute the phase in praxis, when the FFT is
> computed from Non-Integer number of signal periods?

Replace "praxis" with "practice"?

I'm not absolutely sure of this; but maybe the
frequency interpolation obtained from time domain
zeropadding will work.

Hope this helps.

Greg

On Jul 13, 5:35=A0am, "Afinko" <afinko@n_o_s_p_a_m.gmail.com> wrote:
> Hi,
>
> I am trying to compute the amplitude and phase of signal with specific
> frequency. The frequency of the signal is known (measured with laser - th=
e
> rotation machine). I tried to simulate the computation in matlab first:
>
> close all
> clear all
> % clc
>
> f =3D 50; =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 % Frequency of signal
> fs =3D 64000; =A0 =A0 =A0 =A0 =A0 =A0 % Sample frequency
> % fs =3D 2560;
> phase_deg =3D 10; =A0 =A0 =A0 =A0 % Phase of signal [degre]
> t =3D 0.5; =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0% Time of signal
>
> x =3D (0:1/(fs):t-(1/fs)); =A0 =A0 =A0 =A0 =A0 =A0% Sampled time
> phase_rad =3D 2*pi*phase_deg/360; =A0 =A0 % Phase in radians
> y =3D cos(2*pi*f*x + phase_rad); =A0 =A0 =A0% Generate sampled signal
> figure;plot(x,y); =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 % Plot signal
>
> % y =3D y + (0.1*randn(size(y))); =A0 =A0 =A0 % Add some noise
>
> % Apply windowing function on the signal
> w =3D 2.*hanning(length(y));
> y =3D y.*w';
> figure;plot(x,y);
>
> % Compute FFT and plot abs and angle
> Y =3D fft(y);
> figure;plot(abs(Y));
> figure;plot(angle(Y));
>
> % Find the signal in frequency domain
> [max_num,max_pos] =3D max(abs(Y));
>
> % Value of signal (Amplitude)
> value_x =3D abs(Y(max_pos))/(length(x)/2)
>
> % Angle of signal [degre]
> angle_x =3D 360*(angle(Y(max_pos)))/(2*pi)
>
> I used the FFT, found the main peak (that represent the signal of
> interest), and compute the amplitude and phase of this peak.
>
> Everything works great, even a small noise don't change the value of
> amplitude and phase reasonable.
>
> However, the problem occurs when the FFT is computed from Non-Integer
> number of signal periods. What is in praxis all the time. Try to change t=
he
> frequency of the signal to e.g. 50.4, 51.7, 53.9 etc. The classical
> approach should be to use the windowing function. I used Hanning (i tried
> lot of them) in example. Windowing will help with amplitude. But the phas=
e
> of the signal is computed totally incorrectly, even when the windowing is
> used. I found that the windowing function do not play any (or not
> reasonable) role in phase computation.
>
> Is there any way how to compute the phase in praxis, when the FFT is
> computed from Non-Integer number of signal periods?
>
> Thanks in advance.
> Afinko

use a single point DFT that tracks with frequency, a tracking filter

On Jul 13, 2:35=A0am, "Afinko" <afinko@n_o_s_p_a_m.gmail.com> wrote:
....

> Is there any way how to compute the phase in praxis, when the FFT is
> computed from Non-Integer number of signal periods?
> ...

> Afinko

A group of methods called "interpolated fft" can be used to calculate
frequency, amplitude, and phase from the complex fft coefficients of
the peak magnitude bin and the adjacent bins on each side. The most
accurate results are achieved when the interpolation formula used is
designed to match the window applied to the data.

Examples:

http://cdsweb.cern.ch/record/738182/files/ab-2004-023.pdf

http://wwwtw.vub.ac.be/elec/Papers%20on%20web/Papers/JohanSchoukens/IM92Sch=
oukens-The%20Interpolated.pdf


Dale B. Dalrymple

On Jul 14, 5:39=A0am, "Afinko" <afinko@n_o_s_p_a_m.gmail.com> wrote:
> >if the exact frequency, f0, is known in advance, i would just
> >correlate the sampled signal to cos(2*pi*f0*n) and sin(2*pi*f0*n),
> >after LPFing both separately, apply the 4 quadrant arg{} for phase
> >(amplitude is from square rooting the sums of squares).
>
> hi robert bristow-johnson,
> Thanks for your reply.
> The problem is, that in measured signal, there is not only this one
> examined frequency component, but there is a lot of another frequency
> components, with much higher amplitudes. The best solution would be to
> filter this signal first,

no, the best solution is, *if* the frequency of interest is known in
advance, to heterodyne the signal down to DC first where ...

> but such a sharp bandpass filter is a problem to
> achieve with FIR,

   ... a sharp LPF is not a problem, either FIR or IIR.

r b-j

On Jul 14, 2:16=A0am, "Afinko" <afinko@n_o_s_p_a_m.gmail.com> wrote:
> >If you reference the phase to the center of your FFT aperture
> >(by means of an pre-fftshift, or post flipping of alternate
> >bins), the evenness to oddness ratio of sinusoid can be
> >interpolated, even if it is not an integer number of signal
> >periods in the FFT aperture. =A0And the evenness to oddness
> >ratio is cleanly related to the atan2() phase.
>
> >The idiocy where people say you can't interpolate phase
> >is because they are using the start of the FFT aperture;
> >and a non-integer period sinusoid in discontinuous at the
> >boundaries of an FFT aperture, so trying to reference the
> >phase to a discontinuity is certainly non-pictorial and
> >non-intuitive.
>
> >Also, most windows go to zero at the boundaries of the
> >FFT aperature. =A0So you end up looking for the phase of
> >nothing. =A0The signal, after windowing, is in the middle.
>
> >So just move your phase reference to the center.
>
> Hi Ron N.,
> Thanks for your reply.
> I read your post several times, but I do not understand what exactly you
> mean by "move your phase reference to the center".
> Can you pleas post some MATLAB example how to do this?
> Or can you modify my MATLAB code and show how does it work?

Read about and try using matlabs fftshift command before
the FFT.

It moves the phase reference (actually it leaves the phase
reference at the boundary, but rotates the data within the
aperture, giving the same result as moving the phase
reference.)


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

On Wed, 14 Jul 2010 09:05:11 -0700 (PDT), dbd <dbd@ieee.org> wrote:

>On Jul 13, 2:35�am, "Afinko" <afinko@n_o_s_p_a_m.gmail.com> wrote:
>...
>
>> Is there any way how to compute the phase in praxis, when the FFT is
>> computed from Non-Integer number of signal periods?
>> ...
>
>> Afinko
>
>A group of methods called "interpolated fft" can be used to calculate
>frequency, amplitude, and phase from the complex fft coefficients of
>the peak magnitude bin and the adjacent bins on each side. The most
>accurate results are achieved when the interpolation formula used is
>designed to match the window applied to the data.
>
>Examples:
>
>http://cdsweb.cern.ch/record/738182/files/ab-2004-023.pdf
>
>http://wwwtw.vub.ac.be/elec/Papers%20on%20web/Papers/JohanSchoukens/IM92Schoukens-The%20Interpolated.pdf
>
>
>Dale B. Dalrymple
 
    Hi Dale, thanks for those references.  I've only 
had a look at the Gasior and Gonzalez paper.  However, 
I'm suspicious of their Figure 1.  I wonder what is 
the meaning of the dashed curve in their Figure 1.

If that dashed curve is the *true* spectrum (the 
discrete-time Fourier transform) of a finite-duration 
input sinewave, then the mainlobe null-to-null width 
of that dashed curve should be only two bin-widths. 
However, their dashed curve has a mainlobe width of 
roughly five bin-widths.  That doesn't seem 
correct to me.

Also, if the input sinewave's freq is located exactly 
on a bin center, then the max spectral magnitude component 
will be large and its two neighboring spec magnitude 
components should be *VERY* small (ideally zero).
But that's not what they show in their Figure 1(a).

Am I missing something here?

See Ya',
[-Rick-]


   

On Sat, 17 Jul 2010 12:03:25 -0700, Rick Lyons
<R.Lyons@_BOGUS_ieee.org> wrote:

   [Snipped by Lyons]
>
>If that dashed curve is the *true* spectrum (the 
>discrete-time Fourier transform) of a finite-duration 
>input sinewave, then the mainlobe null-to-null width 
>of that dashed curve should be only two bin-widths. 
>However, their dashed curve has a mainlobe width of 
>roughly five bin-widths.  That doesn't seem 
>correct to me.

Oops, in the above where I wrote:

   "... mainlobe width of roughly five bin-widths ...",

I should have said:

   "... mainlobe width of roughly *four* bin-widths ..."

See Ya',
[-Rick-]

On Jul 17, 12:03=A0pm, Rick Lyons <R.Lyons@_BOGUS_ieee.org> wrote:
>...
> =A0 =A0 Hi Dale, thanks for those references. =A0I've only
> had a look at the Gasior and Gonzalez paper. =A0However,
> I'm suspicious of their Figure 1. =A0I wonder what is
> the meaning of the dashed curve in their Figure 1.
>
> If that dashed curve is the *true* spectrum (the
> discrete-time Fourier transform) of a finite-duration
> input sinewave, then the mainlobe null-to-null width
> of that dashed curve should be only two bin-widths.
> However, their dashed curve has a mainlobe width of
> roughly five bin-widths. =A0That doesn't seem
> correct to me.
>
> Also, if the input sinewave's freq is located exactly
> on a bin center, then the max spectral magnitude component
> will be large and its two neighboring spec magnitude
> components should be *VERY* small (ideally zero).
> But that's not what they show in their Figure 1(a).
>
> Am I missing something here?
>
> See Ya',
> [-Rick-]

Hi Rick

The red-dashed spectrum in figure 1 if the G&G paper is the shape of
the frequency response of the window function. You may note that they
do not include the rectangular window in their collection of windows.
In their first paragraph of conclusions they suggest that parabolic
interpolation does not start to have "significant" accuracy until the
window frequency mainlobe width is at least 3 bins.

We could have a long discussion of whether the assumption that the
choice of fft size equal to the number of samples and a rectangular
window is "true" or just lazy and ignorant, but comp.dsp does not
serve beer, so we should discuss that issue elsewhere.

The parabolic interpolation is a quick and dirty "know nothing but
that there are 3 points about a peak" approach. To use it with a
rectangular window, the transform size is usually increased by zero-
filling to improve accuracy. The bandwidth of the frequency mainlobe
of the window doesn't change, but the samples in the frequency domain
become closer together: more samples per mainlobe width. The parabolic
interpolater is an exact match for the log of the power spectrum with
a Gaussian window, but some error sneaks back in when we sample and
truncate the window to finite extent. The errors so introduced are
smaller for Gaussian windows with high sidelobe rejection (and thus
small tails to truncate).The bottom entries in Table 2 illustrate
this.

"Interpolated FFT" results are common in the Transactions on
Instrumentation  and Measurement. Other extensions include estimators
of greater than 3 points and estimators designed for specific windows.

Dale B. Dalrymple

On Jul 17, 3:10=A0pm, dbd <d...@ieee.org> wrote:
> We could have a long discussion of whether the assumption that the
> choice of fft size equal to the number of samples and a rectangular
> window is "true" or just lazy and ignorant, but comp.dsp does not
> serve beer, so we should discuss that issue elsewhere.
>
> The parabolic interpolation is a quick and dirty "know nothing but
> that there are 3 points about a peak" approach. To use it with a
> rectangular window, the transform size is usually increased by zero-
> filling to improve accuracy.

A couple variations:  Since zero-padding is roughly equivalent
to Sinc interpolation, you could just locally interpolate a
few samples before using parabolic interpolator.  The other
trick it to pick a window whose transform is very similar to
a parabola in shape, which, of course, should reduce the error
of a parabolic interpolator.


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

On Jul 18, 12:51=A0am, "Ron N." <rhnlo...@yahoo.com> wrote:

> ...
> A couple variations: =A0Since zero-padding is roughly equivalent
> to Sinc interpolation, you could just locally interpolate a
> few samples before using parabolic interpolator. =A0

When the mainlobe shape is simple and well enough known for easy
accurate local interpolation, the same characteristics can be applied
to estimating the position of the peak directly. That is what the
window specific interpolators do.

>The other
> trick it to pick a window whose transform is very similar to
> a parabola in shape, which, of course, should reduce the error
> of a parabolic interpolator.
>

This is done with the log of the Gaussian windowed power spectrum, as
shown in the G&G paper.

Some authors have taken this a different direction by designing a
window that has a triangular mainlobe and using a peak estimator
designed for the triangle shape.

Dale B. Dalrymple

On Sat, 17 Jul 2010 15:10:56 -0700 (PDT), dbd <dbd@ieee.org> wrote:

   [Snipped by Lyons]
>
>Hi Rick
>
>The red-dashed spectrum in figure 1 if the G&G paper is the shape of
>the frequency response of the window function. You may note that they
>do not include the rectangular window in their collection of windows.
>In their first paragraph of conclusions they suggest that parabolic
>interpolation does not start to have "significant" accuracy until the
>window frequency mainlobe width is at least 3 bins.
>
>We could have a long discussion of whether the assumption that the
>choice of fft size equal to the number of samples and a rectangular
>window is "true" or just lazy and ignorant, but comp.dsp does not
>serve beer, so we should discuss that issue elsewhere.
>
>The parabolic interpolation is a quick and dirty "know nothing but
>that there are 3 points about a peak" approach. To use it with a
>rectangular window, the transform size is usually increased by zero-
>filling to improve accuracy. The bandwidth of the frequency mainlobe
>of the window doesn't change, but the samples in the frequency domain
>become closer together: more samples per mainlobe width. The parabolic
>interpolater is an exact match for the log of the power spectrum with
>a Gaussian window, but some error sneaks back in when we sample and
>truncate the window to finite extent. The errors so introduced are
>smaller for Gaussian windows with high sidelobe rejection (and thus
>small tails to truncate).The bottom entries in Table 2 illustrate
>this.
>
>"Interpolated FFT" results are common in the Transactions on
>Instrumentation  and Measurement. Other extensions include estimators
>of greater than 3 points and estimators designed for specific windows.
>
>Dale B. Dalrymple

Hi Dale,
  Thanks a lot for your detailed explanation.  
I saw those words: "...for efficient interpolation the 
minimum width of the spectral peak is 3 bins.", but it 
didn't sink into my head that what they really meant was:

   "For accurate (not efficient) interpolation, we need 
    at least three samples within the mainlobe of the
    windowed signal's "true" spectrum."

Ya' know, this morning (before I read this post of yours) 
I thought, "Maybe those guys were zero-padding their original 
input time sequence to have more samples within the 
signal's mainlobe."  But now, thanks to you, I'm reminded 
that the original input sinusoid can also be windowed in 
order to widen its mainlobe.

See Ya',
[-Rick-]