practical FFT

Dear Steve, Cedron, Eric,

I am wondering if I have a requirement to transform N sample inputs and
the FFT I have is a N/2 pt. FFT. 

In this case, can I pass first N/2 samples to instance 1 and instance 2 in
parallel (in HW world) or call FFT twice to get the transform. I guess
this will not work as the basis functions in the N/2 pt. FFT versus N pt.
FFT are different.

---------------------------------------
Posted through http://www.DSPRelated.com

>Dear Steve, Cedron, Eric,
>
>I am wondering if I have a requirement to transform N sample inputs and
>the FFT I have is a N/2 pt. FFT. 
>
>In this case, can I pass first N/2 samples to instance 1 and instance 2
in
>parallel (in HW world) or call FFT twice to get the transform. I guess
>this will not work as the basis functions in the N/2 pt. FFT versus N
pt.
>FFT are different.
>
>---------------------------------------
>Posted through http://www.DSPRelated.com

Depends on what you are trying to do, and the frequency range that is of
interest to you.

You could also sample every other point.  If your signal contains only
lower frequencies ( < N/4 ), then the first N/4 bins will be very close in
your N/2 FFT as they would be in a full N FFT.

If you want to try to combine the two FFTs as you described, you won't get
one that is similar to a full N FFT, but it can be done, sort of.  Because
of phase issues, you can't simple add the two smaller FFTs, but you can
get an RMS value by averaging the squares of the corresponding bin
magnitudes and then take the square roots.

Ced
---------------------------------------
Posted through http://www.DSPRelated.com

On Wed, 06 Jan 2016 09:03:34 -0600, "Sharan123" <99077@DSPRelated>
wrote:

>Dear Steve, Cedron, Eric,
>
>I am wondering if I have a requirement to transform N sample inputs and
>the FFT I have is a N/2 pt. FFT. 
>
>In this case, can I pass first N/2 samples to instance 1 and instance 2 in
>parallel (in HW world) or call FFT twice to get the transform. I guess
>this will not work as the basis functions in the N/2 pt. FFT versus N pt.
>FFT are different.

You just need to do one more stage using the outputs of the N/2 FFTs
to compute the length-N FFT.    The most common decimation-in-time FFT
algorithms break down each stage as groups of smaller transforms of
previous stages, so the last stage computes the length-N output from
two length N/2 inputs.

The computation of the last stage is not tremendously difficult, and
can be found in many DSP texts or books specifically about the FFT.


Eric Jacobsen
Anchor Hill Communications
http://www.anchorhill.com

On Friday, January 1, 2016 at 10:39:29 AM UTC-6, Sharan123 wrote:
> I am having some questions related to FFT while using inbuilt Matlab or
> Octave FFT functions.
> 
> These functions seem to take N samples of the signal that needs to be
> analyzed. 
> 
> Are there any requirement on the length of the input signal that needs to
> be fed into.

The most easily-implemented transforms are for those values of N possessing no large prime factors (the case where the only prime factor of N is 2 being the most commonly-cited and easiest-to-implement one). The matrix for DFT(N) that leads to an FFT factors as N factors. So, for instance DFT(144) = DFT(3) DFT(3) DFT(4) DFT(4) (actually more efficient to do DFT(4) instead of DFT(2) DFT(2)). The one complication is that additional factors are needed if you use the most common way of indexing the matrix components in the intermediate states.

I posed an implementation of the DFT for arbitrary N here a year or two ago, you might want to do a search on it.

As for large primes (which I don't handle in my version), the usual approach is to re-express DFT(p) for a large prime p as a convolution problem, and express this IN TURN, in terms of a DFT (since convolution is easily expressible in terms of the Fourier components), where the size of the 2nd stage DFT is made into a convenient composite number large enough to handle the convolution.

The best understanding of the DFT and what values are associated with its components comes by restoring the units and expressing everything as a sum patterned after the integral relations. For convenience (as I always do) I will write exp(2 pi i x) as 1^x (since exp(2 pi i) = 1). The integral relations:
   f^(nu) = integral f(t) 1^{-nu t} dt
   f(t) = integral f^(nu) 1^{nu t} d(nu)
where nu denotes the CYCLIC frequency (cycles per second), not ANGULAR frequency, become in finite form the following relations:
   f^(nu) = sum_{t in T} f(t) 1^{-nu t} Dt
   f(t) = sum_{nu in F} f^{nu} 1^{nu t} D(nu)
where the finite increments Dt and D(nu) are subject to the identity
   Dt D(nu) = 1/N.
(The "quantization" condition for a DFT(N)).

The domains T and F are what your question pertains to. For the above formulas to work, they have to both be regularly sampled domains and are therefore of the form:
   T = { t_0 + (k + 1/2) Dt: k in N }
   F = { nu_0 + (n + 1/2) D(nu): n in N } 
where I use the statement "n in N" to mean "n = 0, 1, 2, ..., N - 1".

The corresponding matrices would be
   f[k] = f(t_0 + (k + 1/2) Dt),
   corresponding to the half-open interval
   [t_0 + k Dt, t_0 + (k + 1) Dt)
   centered on t_0 + (k + 1/2) Dt.

   F[n] = f^(nu_0 + (n + 1/2) D(nu)).
   corresponding to the half-open interval
   [nu_0 + n D(nu), nu_0 + (n + 1) D(nu))
   centered on nu_0 + (n + 1/2) D(nu).

Now ... as to your question, the time bandwidth covered by this transform is the union of the intervals for k in N, which comes out to
   [t_0, t_0 + N Dt) = [t_0, t_0 + 1/(D(nu))).
This interval starts at time t_0 (which is usually taken to be 0), and has a bandwidth 1/(D(nu)). That's the duration of the signal to be transformed.

As far as the usual treatment of the transform makes it specific, I think the usual convention amounts to taking t_0 = -Dt/2, so that the k interval would be
   [(k - 1/2) Dt, (k + 1/2) Dt),
   centered on k Dt.

But it's actually better to take t_0 = 0, for the sake of symmetry, or to leave it as an adjustable parameter. The choice of t_0 (and nu_0) affects the transform -- as described below -- by the inclusion of additional phase factors.

The frequency bandwidth is likewise determined as the union of its respective intervals:
   [nu_0, nu_0 + N D(nu)) = [nu_0, nu_0 + 1/Dt)
so that it's bandwidth is 1/Dt. The prevailing conventions for choosing nu_0 and interpreting the frequencies are the following:
(1) nu_0 = 0, bandwidth = [0, 1/Dt).
   The interval [1/(2Dt), 1/Dt)
   is understood to refer to the NEGATIVE frequencies
   residing on the interval [-1/(2Dt), 0)
   obtained by shifting the above interval back by -1/Dt
(2) nu_0 = -1/(2Dt), bandwidth = [-1/(2Dt), +1/(2Dt)).
   No reinterpretation is required.

Again, I think the more usual version of (1) actually amounts to an offset by -D(nu)/2, so that the n interval would be
   [(n - 1/2) D(nu), (n + 1/2) D(nu))
   centered on n D(nu).
But this makes a mess of things for n = 0! So I prefer either (1) ... or better yet (2).

In all cases, if f(t) is real-valued then f^(nu) and f^(-nu) are conjugates of one another, so that the negative frequency part of the spectrum is redundant.

As to Dt and D(nu), the usual conventions unfortunately obfuscate the physical units attached to the quantities in question and complicate the picture. They are:
(1) Dt = 1, D(nu) = 1/N
(2) Dt = root(1/N) = D(nu), or
(3) Dt = 1/N, D(nu) = 1
with (1) being the most common one used.

The FFT will only directly implement a special case of the above DFT, which is slightly more general. If you write the factors in terms of k and n, you get
   F{n] = sum_{k in N} f[k] 1^{-t nu} Dt
with
   t nu = (t_0 + (k+1/2) Dt)(nu_0 + (n+1/2) D(nu))
   = t_0 nu_0 + t_0 D(nu) (n+1/2)
   + nu_0 Dt (k+1/2) + (k+1/2)(n+1/2)/N

This splits into
  t nu = A + B k + C nu + kn/N
for constant A, B, C given by
   A = (t_0 + Dt/2) (nu_0 + D(nu)/2)
   B = (t_0 + Dt/2) D(nu)
   C = (nu_0 + D(nu)/2) Dt
so
   1^{-t nu} = 1^{-A} 1^{-B k} 1^{-C n} 1^{-kn/N}

Only the nk/N part is directly implemented through the FFT. The other phase factors would have to be included with f[k] or F[n], e.g.
   F[n] 1^{Cn} = sum_{k in N} (f[k] 1^{-A-Bk}) 1^{-kn/N}

If one uses t_0 = -Dt/2, and nu_0 = -D(nu)/2 then the above decomposition works out to
   t nu = (k Dt)(n D(nu)) = kn/N
so that the DFT formula would read
   F[n] = sum_{k in N} f[k] 1^{-kn/N}.

But you're better off keeping the convention t_0 = 0 and nu_0 = 0 (or nu_0 = -1/(2 Dt)) in order to have a more consistent interpretation, even at the cost of adding the extra phase factors.

In the program WtoB.c (which I mentioned several articles back, and which I'm making freely available -- send e-mail note if you want a copy) I use the above convention in general form, except for making t_0 = 0. So the resulting spectrograph is actually a ARBITRARY frequency band [nu_0, nu_0 + 1/Dt), where Dt is the sound source file's sampling rate; the most common values being 1/44100 second 1/48000 second.

>No worries.  Here's an explanation that should help:
>
>cycles per time unit = Hz (usually)
>
>samples per time unit = Sampling Rate
>
>samples per frame = N
>
>cycles per frame = bin number
>
>frames per time unit = frame rate = Sampling Rate / N
>
>Hz = bin number * ( Sampling Rate / N )
>
>Hz (cycles per second) stands for your actual frequency.  It could just
as
>well been RPM (revolutions per minute), depends on your time unit.

Dear Cedron,

I am working out a sample calculation below:

cycles per time unit = 100 (100 Hz signal)
samples per time unit = 200 (twice highest frequency)
samples per frame = N = 200 (one frame contains one complete cycle of 100
Hz signal)
frames per time unit = frame rate = sampling rate/N = 200/200
Hz = bin no * (sampling rate/N) = Max bin number * (200/200) = 200 *
(200/20) = 200

Above, I have assumed that Max bin number for a frame size of 200 = 200

I don't know if my calculations are correct, as I expected higher
frequency to be 100 Hz instead of 200 Hz.
---------------------------------------
Posted through http://www.DSPRelated.com

>Dear Cedron,
>
>I am working out a sample calculation below:
>
>cycles per time unit = 100 (100 Hz signal)
>samples per time unit = 200 (twice highest frequency)
>samples per frame = N = 200 (one frame contains one complete cycle of
100
>Hz signal)
>frames per time unit = frame rate = sampling rate/N = 200/200
>Hz = bin no * (sampling rate/N) = Max bin number * (200/200) = 200 *
>(200/20) = 200
>
>Above, I have assumed that Max bin number for a frame size of 200 = 200
>
>I don't know if my calculations are correct, as I expected higher
>frequency to be 100 Hz instead of 200 Hz.


Dear Cedron,

One more thing I just thought of adding to the above steps -

Hz per bin = Sampling rate / N 
             In the above case, it is 200/200 = 1 Hz per bin
---------------------------------------
Posted through http://www.DSPRelated.com

>
>Dear Cedron,
>
>I am working out a sample calculation below:
>
>cycles per time unit = 100 (100 Hz signal)
>samples per time unit = 200 (twice highest frequency)
>samples per frame = N = 200 (one frame contains one complete cycle of
100
>Hz signal)
>frames per time unit = frame rate = sampling rate/N = 200/200
>Hz = bin no * (sampling rate/N) = Max bin number * (200/200) = 200 *
>(200/20) = 200
>
>Above, I have assumed that Max bin number for a frame size of 200 = 200
>
>I don't know if my calculations are correct, as I expected higher
>frequency to be 100 Hz instead of 200 Hz.
>---------------------------------------
>Posted through http://www.DSPRelated.com

Ok, what you are trying to calculate is a bit confusing.  If you are
saying that you have a 100Hz signal, then I'm going to assume you are
talking about a sinusoidal tone and not your random number generator.  But
your next line uses the term "highest frequency".

For a sampling rate of 200 samples per second, the upper limit of
frequencies, the Nyquist limit, is 100Hz.  The frame size doesn't matter,
you are still going to be at two samples per cycle.  You can not get a
meaningful answer at the limit.  For the best frequency resolution using
my formula you want the frequency to be in the neighborhood of half the
Nyquist limit.

So let's work with the numbers as you have them, ignoring your 100Hz
signal specification.

Samples per second = 200

For audio, that is really slow, but this may not be audio.

Samples per frame = 200

Now you know your frame rate.

Frames per second = 1

Suppose you take the DFT.  The bin values range from 0 to 199, where bin
100 is the Nyquist limit.  So the "meaningful" bins are 0 through 99.

Suppose you have a spike at bin 38, and all the other bin values are
basically zero.  This means there is a 38 ( = 38 * (200/200) ) Hz signal
present in your frame.

Suppose instead that there is a peak at bin 54, with nearly as big a value
at bin 55, and the magnitudes fall off on either side.  You also note that
the value of bin 54 is nearly in the opposite direction of the value of
bin 55 on the complex plane.  This means if you have a sinusoidal signal,
then it has a frequency between 54 and 55 Hz, a little closer to 54.

On the other hand, if you are talking about analysing your random noise
generator, unless there is a significant spike in any bin value, any small
peaks are not going to be indicative of a pattern in your number stream. 
If there is a spike, or a peak like the one described previously,
determining an accurate frequency and phase is no going to be possible.

This set of equations was my answer to your question of how do you
correlate a bin number to a frequency per unit time.

If you are trying to determine if your random number generator is actually
random, rather than somehow measuring how random it is, I present the
following code snippet in BASICA (or gwbasic as it later became known).

randomize timer  ' reset the PRG

screen 1  ' 320 x 480 four colors

for i& = 1 to 1000000
 x% = int( rnd * 320 )
 y% = int( rnd * 480 )
 c% = 1 + int( 3 * rnd )
 pset (x%,y%), c%
next

Did it from memory, so please forgive any typos.  You'll want to implement
on a more modern platform anyway.

The point is your eyes are very good at picking out patterns.  If there is
any pattern in your number stream, it will show in the display.  I chose
the basica example, because at the time you would get a set of clearly
discernible diagonal lines on the screen.  The PRG wasn't very good at
all.

PRGs (Pseudo Random Generators) are heavily studied in the field of
cryptography because the security of many ciphers depends on the
randomness of the number stream.  They have ways to measure that are
significantly superior to using a DFT.

Ced
---------------------------------------
Posted through http://www.DSPRelated.com

>Ok, what you are trying to calculate is a bit confusing.  If you are
>saying that you have a 100Hz signal, then I'm going to assume you are
>talking about a sinusoidal tone and not your random number generator. 
But
>your next line uses the term "highest frequency".
>
>For a sampling rate of 200 samples per second, the upper limit of
>frequencies, the Nyquist limit, is 100Hz.  The frame size doesn't
matter,
>you are still going to be at two samples per cycle.  You can not get a
>meaningful answer at the limit.  For the best frequency resolution using
>my formula you want the frequency to be in the neighborhood of half the
>Nyquist limit.
>
>So let's work with the numbers as you have them, ignoring your 100Hz
>signal specification.
>
>Samples per second = 200
>
>For audio, that is really slow, but this may not be audio.
>
>Samples per frame = 200
>
>Now you know your frame rate.
>
>Frames per second = 1
>
>Suppose you take the DFT.  The bin values range from 0 to 199, where bin
>100 is the Nyquist limit.  So the "meaningful" bins are 0 through 99.
>
>Suppose you have a spike at bin 38, and all the other bin values are
>basically zero.  This means there is a 38 ( = 38 * (200/200) ) Hz signal
>present in your frame.
>
>Suppose instead that there is a peak at bin 54, with nearly as big a
value
>at bin 55, and the magnitudes fall off on either side.  You also note
that
>the value of bin 54 is nearly in the opposite direction of the value of
>bin 55 on the complex plane.  This means if you have a sinusoidal
signal,
>then it has a frequency between 54 and 55 Hz, a little closer to 54.

Dear Cedron,

My sincere apologies. My question above has nothing to do with PRNG I am
discussing in the other thread.

So, to understand better. Are half the bins at the FFT out always sort of
throw away?. Because this is the fundamental question I got after trying
out the above example. In the example above, I assumed a single tone
signal at 100 Hz and sampled at twice its frequency, which is 200 samples
(roughly 2 samples per cycle).

Also, the Hz calculation I am showing is for the highest bin number after
FFT. 

One more question - for a given highest frequency component in a signal or
for a single tone, what is the implication of increasing the sample rate
beyond Nyquist rate? (say 4*Nyquist/6*Nyquist etc. rate). 

I see that except the resolution things remain consistent?

---------------------------------------
Posted through http://www.DSPRelated.com

>
>Dear Cedron,
>
>My sincere apologies. My question above has nothing to do with PRNG I am
>discussing in the other thread.
>
>So, to understand better. Are half the bins at the FFT out always sort
of
>throw away?. 

Yes, for real valued signals.  For complex valued signals interpreting the
upper half of the DFT as "negative frequencies", is more meaningful.

>Because this is the fundamental question I got after trying
>out the above example. In the example above, I assumed a single tone
>signal at 100 Hz and sampled at twice its frequency, which is 200
samples
>(roughly 2 samples per cycle).
>

Right at the Nyquist frequency, which is useless.  Suppose your sampling
hit at every zero crossing, it would appear that there is no signal at
all.

Signals with a higher frequency than the Nyquist frequency will appear to
be a signal of a lower frequency.  This phenomena is called "aliasing".

This also means when you detect a signal within your frame, you can't be
sure it isn't an alias unless you know something about the signal you are
analysing that rules it out.

This is reflected in my frequency formula because the result on an inverse
cosine function can have many values.

>Also, the Hz calculation I am showing is for the highest bin number
after
>FFT. 
>
>One more question - for a given highest frequency component in a signal
or
>for a single tone, what is the implication of increasing the sample rate
>beyond Nyquist rate? (say 4*Nyquist/6*Nyquist etc. rate). 

This is sort of an apples and oranges question, which is why I avoided the
term "rate" in my units labels.  The DFT itself does not care what the
sampling rate is.  The Nyquist frequency will always be located at bin
N/2.  (Which corresponds to two samples per signal cycle.)

>
>I see that except the resolution things remain consistent?
>
>---------------------------------------
>Posted through http://www.DSPRelated.com

"Resolution" is another slippery term.  In order to "resolve" two signals
of nearby frequencies you need to increase the sampling interval (include
more cycles), the sampling density has nothing to do with it.

For example, suppose your interval of interest is one second long.

If you have 200 sample points you can correctly measure all pure tones
with a frequency up to 100Hz.  If your signal tone is 3Hz, then bin 3 will
be non-zero and all the rest of the bins will be zero.

If you have 44100 sample points (CD quality), you can correctly find all
pure tones up to 22050Hz.  If your signal tone is 3Hz, then bin 3 will be
non-zero and all the rest will be zero.

A frequency of 110 Hz will appear to be 90Hz in the first example and
correctly 100 Hz in the second.

Hope this helps.  It can be confusing, until you "get it", of course.

Ced
---------------------------------------
Posted through http://www.DSPRelated.com

>
>A frequency of 110 Hz will appear to be 90Hz in the first example and
>correctly 100 Hz in the second.
>

Oops, typo.  It should have read "correctly 110 Hz in the second."

Ced
---------------------------------------
Posted through http://www.DSPRelated.com