AAC | ADPCM | Convolution | DAFx | FFT | IIR | Mixer | MP3 | MPEG | MPEG-4

Greetings,

I have been experimenting with several FFT routines and none of them seem to work like they are
supposed to.

From what I understand of the theory, an impulse in the time domain is supposed to correspond
to a sinusoid in the frequency domain with a frequency the same as the offset of the impulse.

I have been testing using 16 point signals (imaginary parts set to 0) and there are 2 problems
I have observed which are common to all the FFT routines I have tried :-

1. A very large number sometimes appears where the zero crossings should be.

2. The resulting sinusoids do not have the correct frequency (as far as I can tell this only
happens for numbers which are not a power of 2)

Can anyone explain this?
	

(You need to be a member of audiodsp -- send a blank email to audiodsp-subscribe@yahoogroups.com )

Shadow-

> I have been experimenting with several FFT routines and none of them seem to work like
they are supposed to.
> 
> From what I understand of the theory, an impulse in the time domain is supposed to
correspond to a sinusoid in the frequency domain with a frequency the same as the offset of the
impulse.

No -- impulse in time domain is a flat line in freq domain.  The ideal impulse is
infinitely narrow and infinitely tall -- it contains all possible frequencies in
equal amounts.

You need a good DSP book stuff, here is a one:

  http://www.prenhall.com/~dspfirst/

Straightforward and easy to understand, but written by famous experts nevertheless.

-Jeff

(You need to be a member of audiodsp -- send a blank email to audiodsp-subscribe@yahoogroups.com )

Oops - forgot to hit reply all instead of reply :-

I should also have mentioned that the results of inverse FFT are nowhere
 near the same as the original signal - this applies to all routines I 
have tested.  The errors are far too big to count as round off noise.  
Something is clearly wrong here.
	--- Jeff Brower <jbrower@jbro...> wrote:

From: Jeff Brower <jbrower@jbro...>
Date: Sun, 27 Feb 2005 18:55:15 -0600
To: shadowsedge@shad...
Cc: audiodsp@audi...
Subject: Re: [audiodsp] Impulse Testing the FFT

Shadow-

> I have been experimenting with several FFT routines and none of them seem to work like
they are supposed to.
> 
> From what I understand of the theory, an impulse in the time domain is supposed to
correspond to a sinusoid in the frequency domain with a frequency the same as the offset of the
impulse.

No -- impulse in time domain is a flat line in freq domain.  The ideal impulse is
infinitely narrow and infinitely tall -- it contains all possible frequencies in
equal amounts.

You need a good DSP book stuff, here is a one:

  http://www.prenhall.com/~dspfirst/

Straightforward and easy to understand, but written by famous experts nevertheless.

-Jeff
	_____________________________________________________________
Get free email that shows your true colors...
http://www.YellowSnow.com
	

(You need to be a member of audiodsp -- send a blank email to audiodsp-subscribe@yahoogroups.com )

Hi Shadow,
 
      Too add little more to Jeffs comments:
      Time Domain unit step function ---FFT--> unit impluse. This means that input signal
contains only DC component. 
 
     Time Domain unit impluse function ---FFT--> unit step.  This means that unit impluse
function contains all the frequency components............!. Why it contains all the frequency
components already jeff has given hints.
 
Amaresh

Jeff Brower <jbrower@jbro...> wrote:
Shadow-

> I have been experimenting with several FFT routines and none of them seem to work like
they are supposed to.
> 
> From what I understand of the theory, an impulse in the time domain is supposed to
correspond to a sinusoid in the frequency domain with a frequency the same as the offset of the
impulse.

No -- impulse in time domain is a flat line in freq domain.  The ideal impulse is
infinitely narrow and infinitely tall -- it contains all possible frequencies in
equal amounts.

You need a good DSP book stuff, here is a one:

  http://www.prenhall.com/~dspfirst/

Straightforward and easy to understand, but written by famous experts nevertheless.

-Jeff
	   To
	---------------------------------
Do you Yahoo!?
 Yahoo! Sports -  Sign up for Fantasy Baseball.
	

(You need to be a member of audiodsp -- send a blank email to audiodsp-subscribe@yahoogroups.com )

Martin-

> From: "Jeff Brower" <jbrower@jbro...>
> 
> > No -- impulse in time domain is a flat line in freq domain.
> > The ideal impulse is infinitely narrow and infinitely tall
> 
> Alas, both of these statements are inaccurate. That strange
> thing, the Dirac delta is the ideal *continuous-time* impulse. In
> discrete time, the Kronecker delta d[n] = {1, n=0; 0, else} is as
> ideal as it gets. But it is true of both that they contain
> > all possible frequencies in equal amounts
> under the Fourier transform or D(T)FT, respectively -- that's the
> definition of ideality here and leads to the other point: While
> the *magnitude* spectrum of a displaced impulse is indeed flat,
> its transform is a complex exponential in the frequency variable
> as the OP has found. It is a constant exactly when the impulse is
> at the time origin.

Thanks for pointing out my inaccuracies.  I hope Shadow has gained full confidence in
the FFT algorithm :-)

-Jeff

(You need to be a member of audiodsp -- send a blank email to audiodsp-subscribe@yahoogroups.com )

From: "Jeff Brower" <jbrower@jbro...>

> No -- impulse in time domain is a flat line in freq domain.
> The ideal impulse is infinitely narrow and infinitely tall

Alas, both of these statements are inaccurate. That strange
thing, the Dirac delta is the ideal *continuous-time* impulse. In
discrete time, the Kronecker delta d[n] = {1, n=0; 0, else} is as
ideal as it gets. But it is true of both that they contain
> all possible frequencies in equal amounts
under the Fourier transform or D(T)FT, respectively -- that's the
definition of ideality here and leads to the other point: While
the *magnitude* spectrum of a displaced impulse is indeed flat,
its transform is a complex exponential in the frequency variable
as the OP has found. It is a constant exactly when the impulse is
at the time origin.
	Martin
	

(You need to be a member of audiodsp -- send a blank email to audiodsp-subscribe@yahoogroups.com )

> Date: Mon, 28 Feb 2005 07:05:13 -0800 (PST)
> From: Amaresh patil <amaresh_p@amar...>
> Subject: Re: Impulse Testing the FFT
>
> Time Domain unit step function ---FFT--> unit impluse.

Yet another myth. It is the derivative of the Heavyside function
(unit step fct) that is the Dirac delta function. No Fourier transform
is involved.

> Time Domain unit impluse function ---FFT--> unit step.

Wrong point. It is the integral of the Dirac function that is the
unit step function.

> This means that unit impluse function contains all the frequency components

What is an FT coefficient then?

Originally, the Fourier theorem was about decomposing a function in a sum
of sines and cosines of different frequencies. Hence, the Fourier coefficients
are magnitudes of the corresponding base functions.

Best regards,

Andrew

P.S. Speaking about small non-zero coefficients shown in an earlier post
to this thread, there is something wrong with them. They are unrealistically
small for single precision and unrealistically large for double precision,
if IEEE-754 numbers were used. The double precision epsilon is about 1.0E-22
if my memory serves me correct. Perhaps the Fourier matrix (the twiddle
factors) were calculated with lesser accuracy.
	

(You need to be a member of audiodsp -- send a blank email to audiodsp-subscribe@yahoogroups.com )

Andrew-

> > Time Domain unit step function ---FFT--> unit impluse.
> 
> Yet another myth. It is the derivative of the Heavyside function
> (unit step fct) that is the Dirac delta function. No Fourier transform
> is involved.
> 
> > Time Domain unit impluse function ---FFT--> unit step.
> 
> Wrong point. It is the integral of the Dirac function that is the
> unit step function.
> 
> > This means that unit impluse function contains all the frequency components
> 
> What is an FT coefficient then?
> 
> Originally, the Fourier theorem was about decomposing a function in a sum
> of sines and cosines of different frequencies. Hence, the Fourier coefficients
> are magnitudes of the corresponding base functions.

These are great comments.  The finer points are appreciated.

I would however point out that the OP appeared to be entry-level (at least in DSP). 
I continue to be surprised by the number of BSEEs and MSEEs who I interview, and when
I ask them to **sketch** (sketch, not elaborate with theory) the magnitude response
of a delta function, they do not draw anything like a flat line.

For entry-level engineers, I feel like if they know in their heads what to expect
from Fourier transform of a) pure sinusoid, b) delta function, and c) a constant,
then they have some basic idea of when and why an FFT could be helpful in their work,
and they came out of school in good shape.

-Jeff

(You need to be a member of audiodsp -- send a blank email to audiodsp-subscribe@yahoogroups.com )

Hi Andrew,
 
  Thanks for correcting my statements. your absolutely right, Heavy side   
These are practical observation what we come across put it without mapping them to   classical
integral mathematics. If the Pulse is wide in time domain, it  is narrow in the frequency
domain and vice versa. 
 
Once i agree to your statements then i am left with two questions...!
 
 Time Domain unit step function ---FFT--> unit impluse.

Yet another myth. It is the derivative of the Heavyside function
(unit step fct) that is the Dirac delta function. No Fourier transform is involved. 
 
Then what would be the Fourier tranform of (unit) step function. 
 
> Time Domain unit impluse function ---FFT--> unit step.

Wrong point. It is the integral of the Dirac function that is the
unit step function.
 
What will be the fourier transform of impluse function. 
 
What we will get when we compute on the Digital Machines.
 
regards
Amaresh
	

(You need to be a member of audiodsp -- send a blank email to audiodsp-subscribe@yahoogroups.com )

I have to add something here. I believe Jeff was trying to help a novice,
and his comments were very much appropriate for a novice audience.
I don't think the finer, theoretical points made by others, although
technically accurate and welcome by myself and, I'm sure, also by Jeff,
would be much help to a novice. If fact, the novice would only get more
confused by them.

As an analogy, V=RI is universally known as Ohm's Law, but strictly
speaking, this is a simplified version of the original Ohm's Law, in which
current, resistance, etc. are actually functions, and not constants. 
However,
since,in practical terms, V=RI is perfectly accurate 99.999% of the time , I 
don't think
I'd be wrong to tell a novice that ohm's law is simply V=RI.

Tony
----- Original Message ----- 
From: "Jeff Brower" <jbrower@jbro...>
To: "Andrew Nesterov" <andrew.nesterov@andr...>
Cc: <audiodsp@audi...>
Sent: Thursday, March 03, 2005 8:41 AM
Subject: Re: [audiodsp] Re: Impulse Testing the FFT
	>
> Andrew-
>
>> > Time Domain unit step function ---FFT--> unit impluse.
>>
>> Yet another myth. It is the derivative of the Heavyside function
>> (unit step fct) that is the Dirac delta function. No Fourier transform
>> is involved.
>>
>> > Time Domain unit impluse function ---FFT--> unit step.
>>
>> Wrong point. It is the integral of the Dirac function that is the
>> unit step function.
>>
>> > This means that unit impluse function contains all the frequency 
>> > components
>>
>> What is an FT coefficient then?
>>
>> Originally, the Fourier theorem was about decomposing a function in a sum
>> of sines and cosines of different frequencies. Hence, the Fourier 
>> coefficients
>> are magnitudes of the corresponding base functions.
>
> These are great comments.  The finer points are appreciated.
>
> I would however point out that the OP appeared to be entry-level (at least 
> in DSP).
> I continue to be surprised by the number of BSEEs and MSEEs who I 
> interview, and when
> I ask them to **sketch** (sketch, not elaborate with theory) the magnitude 
> response
> of a delta function, they do not draw anything like a flat line.
>
> For entry-level engineers, I feel like if they know in their heads what to 
> expect
> from Fourier transform of a) pure sinusoid, b) delta function, and c) a 
> constant,
> then they have some basic idea of when and why an FFT could be helpful in 
> their work,
> and they came out of school in good shape.
>
> -Jeff
>
>
> 
>
>
>
>
>
>
>
	

(You need to be a member of audiodsp -- send a blank email to audiodsp-subscribe@yahoogroups.com )

I am taking the opportunity to correct the errors
I made, without getting into too much math though :)

First, I forgot that the Fourier image of a derivative of
a function is the Fourier image of the function multiplied
by a i*2*PI*f, thus if d(step(t))/dt = delta(t) and
F(step(t)) = 0.5*delta(f)+1/(i*2*PI*f)
then
F(d(step(t))/dt) = i*2*PI*f*(0.5*delta(f) + 1/(i*2*PI*f)) = 1
note that f # delta(f) = 0

> Date: Thu, 3 Mar 2005 06:46:05 -0800 (PST)
> From: Amaresh patil <amaresh_p@amar...>
> 
> Then what would be the Fourier tranform of (unit) step function.

To answer that question I just looked at my copy of Applied Time
Series Analysis, vol.1 Basic Techniques by R.Otnes and L.Enochson,
Wiley, 1978, where it says that for

x(t)= {1 if t >= 0, and 0 if t < 0} the Fourier image would be
X(f)= 0.5*delta(f) + 1/(i*2*PI*f)

The Fourier transform of x(t)=1 (a constant) is delta(f), or
simply speaking, an infinite impulse at 0 Hz, therefore

> What will be the fourier transform of impluse function.

The Fourier image of delta(t) is 1 for all f.

For the shifted delta function, delta(t-t0) the image is a complex
exponent exp(-i*2*PI*f*t0) or cos(2*PI*f*t0)-i*sin(2*PI*f*t0)
It can easily be seen that the magnitude is still equal to 1.

However, if we consider the image of two delta functions,
(delta(t-t0)+delta(t+t0))/2, it would be cos(2*PI*f*t0)
and vice versa, two deltas in the Fourier domain define
a real harmonic in the function domain.

We can think of it the other way: applying the Dirac function
delta(x) to another function f(x) results in the value of the
function at 0. I will denote it as delta(x)#f(x)=f(0), where
the # sign does NOT mean multiplication. Using this fact
we easily can obtain that the Fourier transform of delta(t) is

delta(t) # exp(-i*2*PI*f*t) = exp(0) = 1
delta(t-t0) # exp(-i*2*PI*f*t) = exp(-i*2*PI*f*t0)
delta(f-f0) # exp(-i*2*PI*f*t) = exp(-i*2*PI*f0*t)

> What we will get when we compute on the Digital Machines.

This question (as well as the above) involves a lot of background
theory. To mention a small portion of it, let's say that the Discrete
Fourier Transform is a zeroth order approximation of the Fourier integral
over a finite interval. The 0th order means that an integral is approximated
as a sum of rectangles f(x)*(x2-x1), where x is located somewhere in the
interval [x1,x2], perhaps in x1 or x2 or (x1+x2)/2.

Eventually this leads to approximation of exp(-i*2*PI*f*t) as a
square N x N matrix with elements F(l,k) = exp(-i*2*PI*l*k/N),
and the Discrete Fourier Transform becomes simply multiplication
of the Fourier matrix and the discrete vector that approximates
the transformed function on the finite interval:

        N-1                           N-1
X(k) = SUM exp(-i*2*PI*l*k/N)*x(l) = SUM F(l,k)*x(k)
        l=0                           l=0

Of course, more accurate approximation of the integral can be used,
e.g. 1st order or trapezoid formula, Simpson's method, but the main
advantage of the 0th order approximation lies is in the high symmetry
of the Fourier matrix. It can be decomposed into a series of matrices
with very small number of non-zero elements, significantly decreasing
the number of operations to obtain the final result. This is the base
for the Fast Fourier Transform methods. The common multiplication of
a vector by an N x N matrix takes about N^2 operations, the FFT methods
reduce this number to something proportional to N*log(N)

The narrowing of the integration interval from the whole real axis to
a finite interval results in convolving the continuous spectrum with
the sync(x) function, sin(x)/x, thus making the infinitely thin
spectral lines to have a finite width and the spectrum becames
perodical.

Approximating the Fourier integral by an infinite summation makes continuous
spectrum to be dicrete and a part of the spectrum in the frequencies
larger than half of discretization (sampling) frequency would be folded,
which is called aliasing, e.g. if Fs is the sampling ferquency (also
called sampling rate), then the folding frequency is F = Fs/2, and
frequencies in the [2F,F],[2F,3F],[4F,3F],[4F,5F],... etc intervals all get
folded into [0,F] interval.

Also, since we cannot perform infinite summation on a computer, reducing the
summation to be finite makes the discrete spectrum periodical.

In total, DFT makes the original spectrum to be discrete, periodical, the
original discrete spectrum lines would be spread over the whole discrete
spectrum and higher frequencies will be aliased.

The other side of the calculations is approximation or real numbers by
machine numbers, which sometimes is referred to as "quantization".
This results in that the spectral values would be restricted to the
range of the computer number system: overflow (when a number is outside
the range) or underflow (when a number is less than the smallest representable 
number) might occur, and the numbers would be exact to the precision of
that particular system, this introduces roundoff errors.

Sometimes a computer's arithmetic would roundoff the overflowed values
to the maximum magninute value possible, which is called saturation. These
saturated values would be known with undetermined error, however the sign
of the error is known.

The wraparound arithmetic in the case of overflow gives completely
undetermined error.
	--
Andrew
	

(You need to be a member of audiodsp -- send a blank email to audiodsp-subscribe@yahoogroups.com )

Hey Tony,

I am absolutely glad to hear that and of course you are absolutely
correct, but remember there were two ways of teaching introduced
in the 1st half of 20th century.

I am referring to what is supposedly known as Bohr's method and Landau's
method. Roughly speaking, the methods differ by the definition of who is
more (I am sorry for the language here) "stupid", the teacher or the pupil.
Both methods gave good results :)

That leaves two options to a novice: either to get not enough inspired
and quit or to be scared to death and quit :)

P.S. Don't remember exactly if Dr.Ohm used in fact the inverse of R,
conductivity sigma=R^-1...

On Thu, 3 Mar 2005, Tony Zampini wrote:

> I have to add something here. I believe Jeff was trying to help a novice,
> and his comments were very much appropriate for a novice audience.
> I don't think the finer, theoretical points made by others, although
> technically accurate and welcome by myself and, I'm sure, also by Jeff,
> would be much help to a novice. If fact, the novice would only get more
> confused by them.
>
> As an analogy, V=RI is universally known as Ohm's Law, but strictly
> speaking, this is a simplified version of the original Ohm's Law, in which
> current, resistance, etc. are actually functions, and not constants. However,
> since,in practical terms, V=RI is perfectly accurate 99.999% of the time , I 
> don't think
> I'd be wrong to tell a novice that ohm's law is simply V=RI.
>
> Tony
	

(You need to be a member of audiodsp -- send a blank email to audiodsp-subscribe@yahoogroups.com )

Andrew--

>without getting into too much math though :)

You were teasing us there werent you? Wow! You some math wizard or 
something?

"I have medical certificate from my doctor certifying am allergic to 
derivatives"---One of the cartoons characters the worlds most over 
simplified DSP Book--the North American Classic DSP and the Microcontroller 
Vy Grover & Deller (The publisher tells me only 10 copies were sold in 
India!)

--Bhooshan.N.Iyer

----Original Message Follows----
From: Andrew Nesterov <andrew.nesterov@andr...>
To: audiodsp@audi..., Amaresh patil <amaresh_p@amar...>,  Jeff 
Brower <jbrower@jbro...>
Subject: Re: [audiodsp] Impulse Testing the FFT
Date: Thu, 3 Mar 2005 12:02:48 -0800 (Pacific Standard Time)
	I am taking the opportunity to correct the errors
I made, without getting into too much math though :)

First, I forgot that the Fourier image of a derivative of
a function is the Fourier image of the function multiplied
by a i*2*PI*f, thus if d(step(t))/dt = delta(t) and
F(step(t)) = 0.5*delta(f)+1/(i*2*PI*f)
then
F(d(step(t))/dt) = i*2*PI*f*(0.5*delta(f) + 1/(i*2*PI*f)) = 1
note that f # delta(f) = 0

 > Date: Thu, 3 Mar 2005 06:46:05 -0800 (PST)
 > From: Amaresh patil <amaresh_p@amar...>
 >
 > Then what would be the Fourier tranform of (unit) step function.

To answer that question I just looked at my copy of Applied Time
Series Analysis, vol.1 Basic Techniques by R.Otnes and L.Enochson,
Wiley, 1978, where it says that for

x(t)= {1 if t >= 0, and 0 if t < 0} the Fourier image would be
X(f)= 0.5*delta(f) + 1/(i*2*PI*f)

The Fourier transform of x(t)=1 (a constant) is delta(f), or
simply speaking, an infinite impulse at 0 Hz, therefore

 > What will be the fourier transform of impluse function.

The Fourier image of delta(t) is 1 for all f.

For the shifted delta function, delta(t-t0) the image is a complex
exponent exp(-i*2*PI*f*t0) or cos(2*PI*f*t0)-i*sin(2*PI*f*t0)
It can easily be seen that the magnitude is still equal to 1.

However, if we consider the image of two delta functions,
(delta(t-t0)+delta(t+t0))/2, it would be cos(2*PI*f*t0)
and vice versa, two deltas in the Fourier domain define
a real harmonic in the function domain.

We can think of it the other way: applying the Dirac function
delta(x) to another function f(x) results in the value of the
function at 0. I will denote it as delta(x)#f(x)=f(0), where
the # sign does NOT mean multiplication. Using this fact
we easily can obtain that the Fourier transform of delta(t) is

delta(t) # exp(-i*2*PI*f*t) = exp(0) = 1
delta(t-t0) # exp(-i*2*PI*f*t) = exp(-i*2*PI*f*t0)
delta(f-f0) # exp(-i*2*PI*f*t) = exp(-i*2*PI*f0*t)

 > What we will get when we compute on the Digital Machines.

This question (as well as the above) involves a lot of background
theory. To mention a small portion of it, let's say that the Discrete
Fourier Transform is a zeroth order approximation of the Fourier integral
over a finite interval. The 0th order means that an integral is approximated
as a sum of rectangles f(x)*(x2-x1), where x is located somewhere in the
interval [x1,x2], perhaps in x1 or x2 or (x1+x2)/2.

Eventually this leads to approximation of exp(-i*2*PI*f*t) as a
square N x N matrix with elements F(l,k) = exp(-i*2*PI*l*k/N),
and the Discrete Fourier Transform becomes simply multiplication
of the Fourier matrix and the discrete vector that approximates
the transformed function on the finite interval:

         N-1                           N-1
X(k) = SUM exp(-i*2*PI*l*k/N)*x(l) = SUM F(l,k)*x(k)
         l=0                           l=0

Of course, more accurate approximation of the integral can be used,
e.g. 1st order or trapezoid formula, Simpson's method, but the main
advantage of the 0th order approximation lies is in the high symmetry
of the Fourier matrix. It can be decomposed into a series of matrices
with very small number of non-zero elements, significantly decreasing
the number of operations to obtain the final result. This is the base
for the Fast Fourier Transform methods. The common multiplication of
a vector by an N x N matrix takes about N^2 operations, the FFT methods
reduce this number to something proportional to N*log(N)

The narrowing of the integration interval from the whole real axis to
a finite interval results in convolving the continuous spectrum with
the sync(x) function, sin(x)/x, thus making the infinitely thin
spectral lines to have a finite width and the spectrum becames
perodical.

Approximating the Fourier integral by an infinite summation makes continuous
spectrum to be dicrete and a part of the spectrum in the frequencies
larger than half of discretization (sampling) frequency would be folded,
which is called aliasing, e.g. if Fs is the sampling ferquency (also
called sampling rate), then the folding frequency is F = Fs/2, and
frequencies in the [2F,F],[2F,3F],[4F,3F],[4F,5F],... etc intervals all get
folded into [0,F] interval.

Also, since we cannot perform infinite summation on a computer, reducing the
summation to be finite makes the discrete spectrum periodical.

In total, DFT makes the original spectrum to be discrete, periodical, the
original discrete spectrum lines would be spread over the whole discrete
spectrum and higher frequencies will be aliased.

The other side of the calculations is approximation or real numbers by
machine numbers, which sometimes is referred to as "quantization".
This results in that the spectral values would be restricted to the
range of the computer number system: overflow (when a number is outside
the range) or underflow (when a number is less than the smallest 
representable
number) might occur, and the numbers would be exact to the precision of
that particular system, this introduces roundoff errors.

Sometimes a computer's arithmetic would roundoff the overflowed values
to the maximum magninute value possible, which is called saturation. These
saturated values would be known with undetermined error, however the sign
of the error is known.

The wraparound arithmetic in the case of overflow gives completely
undetermined error.
	--
Andrew

_________________________________________________________________
It's easy to send photos with Hotmail. Click here to find out how:  
http://www.imagine-msn.com/Hotmail/Post/Communicate/SendPhotos.aspx
	

(You need to be a member of audiodsp -- send a blank email to audiodsp-subscribe@yahoogroups.com )