dft : property of symmetry for real & even seq

Hi every one,
% As per theory, if there are real even seq x[n] then it will have Xi
(k) img spectrum to become zero
% for periodic even xp[n] = xp [- n] ; if L=4=x[n] seq length
%  x[-1] =x[1]; x[-2]=x[2] and have to include x[0] additionally
% so that the length of the seq becomes odd . i.e., L=5 with x[0] as
center point
% Instead with 6 given i/p seq x[n] = [1 1 1 1 1 1]
% To check for symmetry property: given x[n] to be real & even seq
% it is enough to compute X(k) for  0 : (N / 2)  if N even. ( DFT have
symmetry)
% then 1+ N/2 points are unique points and the rest is symmetrically
redundant
% if N is odd then X(k) can be computed for 0 : (N-1)/2 . then (N+1)/2
are unique pts and the rest are symmetrically redundant.

code:
L=6;    % i/p signal seq length
x=[1 1 1 1 1 1];    % i/p seq
N=2.^(nextpow2(length(x)));  %N-pt DFT
fftx=fft(x,N);   %calculaing DFT
Numofuniquepoints=ceil((N+1)/2); % Number of Unique Points
fftx=fftx(1:Numofuniquepoints); % Number of Unique Points-DFT
mx=abs(fftx)/length(x);  % scaling fft magnitude so that is not the
function of x
real(fftx);    % real DFT values
imag(fftx);   % imaginary DFT values
o/p:
fftx =
  Columns 1 through 5
   6.0000            -0.7071 - 1.7071i   1.0000 - 1.0000i   0.7071 +
0.2929i        0
  Columns 6 through 8
   0.7071 - 0.2929i   1.0000 + 1.0000i  -0.7071 + 1.7071i

Numofuniquepoints =      5
fftx =    6.0000            -0.7071 - 1.7071i   1.0000 - 1.0000i
0.7071 + 0.2929i        0
mx =     1.0000    0.3080    0.2357    0.1276         0
real(fftx) =     6.0000   -0.7071    1.0000    0.7071         0
imag(fftx) =      0   -1.7071   -1.0000    0.2929         0

Now my doubt is
1.i have X(k) o/p values for 5 seq length instead of length(x[n])
=6;
(whether 50% computation is saved because of removing redundancies,
i.e, DFT symmetry property)
2.Imag part are not exactly zero as x[n] is of real & even seq demands
imaginary spectrum to be zero.

Whether this procedure is wrong
I tried to calculate fft for Numofuniquepoints instead of N points:

code:
L=6;    % i/p signal seq length
x=[1 1 1 1 1 1];    % i/p seq
N=2.^(nextpow2(length(x)));  %N-pt DFT
Numofuniquepoints=ceil((N+1)/2); % Number of Unique Points
fftx=fft(x,Numofuniquepoints); %calculating DFT for saving 50%
computation and storage.
mx=abs(fftx)/length(x);  % scaling fft magnitude so that is not the
function of x
real(fftx);    % real DFT values
imag(fftx);   % imaginary DFT values
o/p:
fftx =  5     0     0     0     0
mx=  5     0     0     0     0
real(fftx)=  5     0     0     0     0
imag(x) = 0     0     0     0     0


Which one is correct ???

Thank  You.

On Mar 8, 10:30 am, rramya <ramyarav...@yahoo.co.in> wrote:
> L=6;    % i/p signal seq length
> x=[1 1 1 1 1 1];    % i/p seq
> N=2.^(nextpow2(length(x)));  %N-pt DFT

By padding up to the next power of 2, you are breaking the symmetry.
There is no reason to pad; an FFT is perfectly fine on 6 points, and
by padding you are changing the transform you are computing.

That is, by padding, you are changing the input array to [1 1 1 1 1 1
0 0]  (length 8), which is not even.

If you really wanted to pad, you woud have to pad in an "even"
fashion.  In this case, you would pad to
                     [1 1 1 0.5 0 0.5 1  1]
which is real-even.   The 0.5 is because one has to do something with
the Nyquist element 1 in your original array, and the symmetric thing
to do is to spit this up equally between positive and negative
frequencies.

It's a little easier to see if your data aren't all 1's.  For example,
the array [1 2 3 8 3 2] is a real-even array, where "1" is the zero-
frequency (DC) element, "8" is the Nyquist element, and "2 3" and "3
2" are the symmetric positive and negative frequency amplitudes,
respectively.  To pad it symmetrically to length 8 (although, as I
said above, you really don't need to), you would do [1 2 3 4 0 4 3 2],
where the Nyquist frequency "8" has been split equally between
positive and negative frequencies.

But really, the simplest advice is: don't pad.  The FFT function in
Matlab is O(N log N) for all N, even for prime N (although prime N are
a few times slower than highly composite N).

Understanding aliasing is the key to understanding symmetry in the
DFT.

Steven

On 8 Mar, 17:53, "Steven G. Johnson" <stev...@alum.mit.edu> wrote:

> &#4294967295;The FFT function in
> Matlab is O(N log N) for all N, even for prime N (although prime N are
> a few times slower than highly composite N).

Just out of curiosity, how does this order reduction work for prime N?
I would assume off the top of my head that the FFT only speeds things
up for composite N?

Rune

On Mar 8, 12:59&#4294967295;pm, Rune Allnor <all...@tele.ntnu.no> wrote:
> On 8 Mar, 17:53, "Steven G. Johnson" <stev...@alum.mit.edu> wrote:
> > &#4294967295;The FFT function in
> > Matlab is O(N log N) for all N, even for prime N (although prime N are
> > a few times slower than highly composite N).
>
> Just out of curiosity, how does this order reduction work for prime N?
> I would assume off the top of my head that the FFT only speeds things
> up for composite N?

That hasn't been true since 1968, when Bluestein's algorithm was
published; another algorithm for prime N (due to Rader) was published
in 1969.  They both work by expressing the DFT in terms of a
convolution, which can be computed exactly (in exact arithmetic) in
terms of DFTs of a different length that is not prime.  I've put
articles on Wikipedia with the details and references (google
Bluestein FFT or Rader FFT).

The root of your error is expressed in the phrase "the FFT".  There
are many distinct FFT algorithms, in addition to the most well known
algorithms (radix-2, mixed-radix, split radix, etc.) which are
variations on Cooley-Tukey.

(I've been repeatedly correcting the popular misconception that N log
N FFT algorithms only exist for composite sizes for over 10 years now
on comp.dsp and elsewhere.  I'm hoping that someday there will be a
critical mass of people who understand that this is not correct, and
the myth will finally die off.)

Steven

On 8 Mar, 18:30, "Steven G. Johnson" <stev...@alum.mit.edu> wrote:
> On Mar 8, 12:59&#4294967295;pm, Rune Allnor <all...@tele.ntnu.no> wrote:

> (I've been repeatedly correcting the popular misconception that N log
> N FFT algorithms only exist for composite sizes for over 10 years now
> on comp.dsp and elsewhere. &#4294967295;I'm hoping that someday there will be a
> critical mass of people who understand that this is not correct, and
> the myth will finally die off.)

Maybe you should talk to the people who write textbooks? That's
where most people are exposed to FFT algorithms. Very few people
roll their own these days, so no one check out the journals.

Rune

On Mar 8, 1:37&#4294967295;pm, Rune Allnor <all...@tele.ntnu.no> wrote:
> Maybe you should talk to the people who write textbooks? That's
> where most people are exposed to FFT algorithms. Very few people
> roll their own these days, so no one check out the journals.

Oh, I do, when I have the opportunity.

The worst offender (due to its popularity) is Numerical Recipes, whose
discussion of FFTs (which apparently hasn't been updated since 1986)
contains several blatant errors in its discussion of algorithms other
than radix-2 Cooley-Tukey, but they've never responded to my emails,
and I've never had a chance to meet the authors in person.  Maybe I
should write them snail mail, but I'm not sure whether they will read
that either (as they probably are deluged with questions/comments and
people urging them to update some section or other).

Ultimately, I'm hoping that the pressure will come from practitioners,
many of whom use Matlab and/or FFTW and/or Intel MKL and should
eventually notice that they are explicitly documented as having N log
N algorithms for all N, including prime N.  As you say, very few
people roll their own FFTs these days (with good reason).  The
difficulty with Matlab is that there are a lot of legacy programs and
tutorials out there that instinctively use the old zero-padding option
in the "fft" function (which didn't support fast prime-length
transforms until Matlab 6 in 2001-ish incorporated FFTW).

Steven

On Mar 8, 10:56&#4294967295;pm, "Steven G. Johnson" <stev...@alum.mit.edu> wrote:
> On Mar 8, 1:37&#4294967295;pm, Rune Allnor <all...@tele.ntnu.no> wrote:
>
> > Maybe you should talk to the people who write textbooks? That's
> > where most people are exposed to FFT algorithms. Very few people
> > roll their own these days, so no one check out the journals.
>
> Oh, I do, when I have the opportunity.
>
> The worst offender (due to its popularity) is Numerical Recipes, whose
> discussion of FFTs (which apparently hasn't been updated since 1986)
> contains several blatant errors in its discussion of algorithms other
> than radix-2 Cooley-Tukey, but they've never responded to my emails,
> and I've never had a chance to meet the authors in person. &#4294967295;Maybe I
> should write them snail mail, but I'm not sure whether they will read
> that either (as they probably are deluged with questions/comments and
> people urging them to update some section or other).
>
> Ultimately, I'm hoping that the pressure will come from practitioners,
> many of whom use Matlab and/or FFTW and/or Intel MKL and should
> eventually notice that they are explicitly documented as having N log
> N algorithms for all N, including prime N. &#4294967295;As you say, very few
> people roll their own FFTs these days (with good reason). &#4294967295;The
> difficulty with Matlab is that there are a lot of legacy programs and
> tutorials out there that instinctively use the old zero-padding option
> in the "fft" function (which didn't support fast prime-length
> transforms until Matlab 6 in 2001-ish incorporated FFTW).
>
> Steven

thank you i understood your point.

1.In the padded symmetric array [1 2 3 4 0 4 3 2],
  this is symmetric with respect to what point? this doubt arises if
one want to visualize
  the symmetric plot.
 ( Actually for circularly symmetric we consider with zero point in
the circle.)

2. In case of real-odd array [ 1  2  3  -3  -2]
(based on the DFT symmetry for odd seq x[N-n]= - x[n] )
In this case there is no nyquist freq so it is made zero.To pad it
symmetrically to length 8,
whether is padded symmetric array is correct? [ 1 2 3 0 0 0 -3 -2 ]
where middle value=nyquist freq=0.
this is symmetric with respect to what point?

3. Since the above odd array N=5, is it neccesary to pad three zeros
symmateically to become N=8
    Instead can one make odd array to become even array by padding one
zero to become N=6?
If it is allowed. then for N=8 case, dense freq components will be
seen in the spectrum (better display)
whereas in case of N=6 , spectral freq resolution will be good.
Am i correct in this regard.

4. You said fft can be done even for prime N. consider N=5; real odd
array [1 2 3 -3 -2]
is there no need to pad zeros symmetrically to N as power of 2 .? then
what about spectrum symmetry?

can you recommend any book or article to study fft ?

Thank you.

On Mar 8, 1:56&#4294967295;pm, "Steven G. Johnson" <stev...@alum.mit.edu> wrote:
> On Mar 8, 1:37&#4294967295;pm, Rune Allnor <all...@tele.ntnu.no> wrote:
>
> > Maybe you should talk to the people who write textbooks? That's
> > where most people are exposed to FFT algorithms. Very few people
> > roll their own these days, so no one check out the journals.
>
> Oh, I do, when I have the opportunity.
>
> The worst offender (due to its popularity) is Numerical Recipes, whose
> discussion of FFTs (which apparently hasn't been updated since 1986)
> contains several blatant errors in its discussion of algorithms other
> than radix-2 Cooley-Tukey, but they've never responded to my emails,
> and I've never had a chance to meet the authors in person. &#4294967295;Maybe I
> should write them snail mail, but I'm not sure whether they will read
> that either (as they probably are deluged with questions/comments and
> people urging them to update some section or other).

Steven, i dunno all of the mistakes in NR, of course, but one would
think that if it's as specific as "There <is>/<is no> FFT alg for
prime N", it would have an objective factuality (hell, it's in my O&S)
and they would have enough sense of the current state of the art to
know about FFTW.  they republish newer versions of NR with changes and
corrections.  do they *intend* to continue to say it wrong?

just curious.

r b-j

On 9 Mar, 04:30, robert bristow-johnson <r...@audioimagination.com>
wrote:
> On Mar 8, 1:56&#4294967295;pm, "Steven G. Johnson" <stev...@alum.mit.edu> wrote:
>
>
>
>
>
> > On Mar 8, 1:37&#4294967295;pm, Rune Allnor <all...@tele.ntnu.no> wrote:
>
> > > Maybe you should talk to the people who write textbooks? That's
> > > where most people are exposed to FFT algorithms. Very few people
> > > roll their own these days, so no one check out the journals.
>
> > Oh, I do, when I have the opportunity.
>
> > The worst offender (due to its popularity) is Numerical Recipes, whose
> > discussion of FFTs (which apparently hasn't been updated since 1986)
> > contains several blatant errors in its discussion of algorithms other
> > than radix-2 Cooley-Tukey, but they've never responded to my emails,
> > and I've never had a chance to meet the authors in person. &#4294967295;Maybe I
> > should write them snail mail, but I'm not sure whether they will read
> > that either (as they probably are deluged with questions/comments and
> > people urging them to update some section or other).
>
> Steven, i dunno all of the mistakes in NR, of course, but one would
> think that if it's as specific as "There <is>/<is no> FFT alg for
> prime N", it would have an objective factuality (hell, it's in my O&S)
> and they would have enough sense of the current state of the art to
> know about FFTW. &#4294967295;they republish newer versions of NR with changes and
> corrections. &#4294967295;do they *intend* to continue to say it wrong?

I've looked up the 3rd edition (2007) of NR (ISBN 978-0-521-88068-8).
In chapter 12.2 they introduce the 'clearest derivation' of the FFT
(page 609). On page 610:

"Although there are ways of treating other cases, by far the easiest
case is the one in which the original N is an integer power of 2"

I can agree so far, if pedagogic clarity is a major factor.
And the above is general enough to include general N, including
prime N. But then the authors take it one step too far:

"In fact, we categorically recommend that you /only/ [sic!]
use FFTs with N a power of 2."

Immediately followed by what seems to be a rebutal of the
previous sentence, rendering the paragraph all but meaningless:

"(We will give more sophisticated suggestions in subsequent
sections below.)"

Then, in section 12.2.3, 'Other FFT algorithms', some more
general FFTs are briefly mentioned, including composite N.
However, on page 616:

"If N is prime, then no subdivision is possible, and the
user (whether he knows it or not) is taking a /slow/ [sic!]
Fourier transform, of order N^2 instead of order N log_2 N.
Our advice is to stay clear of such FFT implementations,
with perhaps the one class of exceptions, the /Winograd
Fourier transform algorithms/."

This is pretty weird stuff, considering that algorithms
for prime N apparently have been out there for 40 years.
I'd consider sending a mail to the publishers of NR if
the authors don't respond.

I always used to consider NR as a good intro text to
get the big picture when one starts out in a new field.
This certainly makes me reconsider!

Rune

Steven G. Johnson <stevenj@alum.mit.edu> wrote:
> On Mar 8, 1:37?pm, Rune Allnor <all...@tele.ntnu.no> wrote:
>> Maybe you should talk to the people who write textbooks? That's
>> where most people are exposed to FFT algorithms. 
 
> Oh, I do, when I have the opportunity.
 
> The worst offender (due to its popularity) is Numerical Recipes, whose
> discussion of FFTs (which apparently hasn't been updated since 1986)
> contains several blatant errors in its discussion of algorithms other
> than radix-2 Cooley-Tukey, but they've never responded to my emails,
> and I've never had a chance to meet the authors in person.  

I suppose as I understand it that is about right.

NR gives you a good start on a variety of numerical topics.
Not just the programs, but some explanation on what went into them.

It is then up to the reader to continue on to other sources when
more depth is needed.  I don't remember exactly how they explain
it, so that is my paraphrase of it.  If you understand the radix-2
algorithm you will be able to read and understand the others from
other sources.  

-- glen