FFT & averaging

Mark Borgerding wrote:
> Andor wrote:
> >>>In simpler words, since the DFT is unitary, the euclidian length
of
> >>>vector is invariant under the DFT.
...
> Ah, yes.  Much simpler.  *POP* *(struggles to extract tongue from
cheek)*

Well, if you are familiar with the definitions, then my point was that
one can state Parseval's relation in simple words (see above).

> Sorry Andor, I just couldn't resist.

No problem. Sometimes, I get carried away by my affinity for linear
algebra ... :-)


> I view Parseval's equality in terms of what it means to the
> practitioner, not the mathematical properties that make it so.
>
> I'm sure both viewpoints are needed at various times.
>

There is another nice property that follows from Parseval. Let's say,
you have some linear transform U, and you have verified Parseval's
relation (ie.(x, x) = (U x, U x), for all x). Then this already
suffices to show that U is unitary:

We need to show (U x, U y) = (x, y) for all vectors x, y. So let x and
y be any non-zero vectors (otherwise the statement is trivial). Write y
= y_1 + y_2, where y_1 is parallel to x, and y_2 is perpendicular to x
(in any Hilbert space, this decomposition is unique), ie. y_1 = x
(x,y)/(x,x), and y_2 = y - y_1. Then you have

(U x, U y)
= (U x, U (y_1 + y_2))
= (U x, U y_1)
= (x, y)/(x, x) (U x, U x)
= (x, y)

where the last equation follows because of Parseval's relation. So
there you have it, U is unitary. This implies again, for example, that
U maps sets of orthonormal vectors to orthonormal vectors (ie. the
angle between the vectors is also invariant under U) - something not
immediately obvious from Parseval's relation! Neat, no? This of course
is the reason why fast correlation works (as the angle between two
vectors is the same in both domains, thanks to unitarity of the DFT).

>  From a practitioner's standpoint, it allows one to make a sliding
power
> detecter by squaring each sample and boxcar filtering that sequence
--
> both operations have minimal cpu requirements independent of the
length
> of the window.  Pretty nifty.

Yup. The whole of DSP is basically about being nifty :-).

> 
> -- Mark

Andor

"Andor" <an2or@mailcircuit.com> wrote in message
news:1115973864.786042.142780@g43g2000cwa.googlegroups.com...
> Bhaskar Thiagarajan wrote:
>
> > > In simpler words, since the DFT is unitary, the euclidian length of
> a
> > > vector is invariant under the DFT.
> >
> > Andor,
> > Perhaps the words are shorter (or the sentence is) but by no means
> are they
> > 'simpler' :-)
> > It took me a while to digest what you said and I'm still not very
> sure if I
> > truly understand it's implication.
>
> Allow me in this case to elaborate by giving the definitions: Let V be
> some n-dimensional (wher n < infinity) complex vector space with an
> inner product. Write the inner product for vectors x and y as (x,y),
> meaning
>
> (x,y) := x^H y =  sum_{k=1}^n x_k* y_k,
>
> where x^H is the adjoint of x (write x as a column vector and x^H as
> the complex conjugate row vector). Let F: V->V be a linear map, and F^H
> be its adjoint (complex conjugate transpose). By definition, if
>
> F^H F = 1        (where "1" denotes the identity map on V)
>
> then F is called unitary. The simple consequence of this definition is
> that
>
> (F x, F y) = (F x)^H (F y) = x^H F^H F y = x^H y = (x, y).
>
> Now if you let F be the DFT, and X = F x, then Parseval's relation is a
> simple consequence of the unitarity of F:
>
> (x, x) = (X, X).
>
> Note that for the DFT to be unitary, you need to scale the matrices F
> and F^H by a factor 1/sqrt(n). This is usual in physics but not so in
> ee, but I find it makes notation a lot more compact.

Given that I made a C in Linear Algebra when I was supposed to study it,
this was primarily a translation of short phrases of difficult words to long
phrases (of difficult words) :-))

Seriously though - I actually tried to understand this (despite my poor
grasp of this space of math...well, who am I kidding...math in general). I
should say, it actually made *some* sense to me. Enough to appreciate this
alternate view of explaining Parseval's theorem - certainly 'nifty'!!

I only wish I'd paid more attention in my linear algebra class instead of
trying to peak at bras (forgive my play on words)...

Cheers
Bhaskar

> Regards,
> Andor

Andor wrote:
> Mark Borgerding wrote:
> > Andor wrote:
> > >>>In simpler words, since the DFT is unitary, the euclidian length
> of
> > >>>vector is invariant under the DFT.
> ...
> > Ah, yes.  Much simpler.  *POP* *(struggles to extract tongue from
> cheek)*
>
> Well, if you are familiar with the definitions, then my point was
that
> one can state Parseval's relation in simple words (see above).
>
> > Sorry Andor, I just couldn't resist.
>
> No problem. Sometimes, I get carried away by my affinity for linear
> algebra ... :-)

I have to side with Andor here. DSP is the art of making
use of linear algebra for solving practical problems.
Lots of people find that to be a contradiciction in terms,
but that's how I view the world.

> > I view Parseval's equality in terms of what it means to the
> > practitioner, not the mathematical properties that make it so.
> >
> > I'm sure both viewpoints are needed at various times.

The mathemathical properties of unitary transforms are what
make Parseval's theorem useful to the practitioner.

> There is another nice property that follows from Parseval. Let's say,
> you have some linear transform U, and you have verified Parseval's
> relation (ie.(x, x) = (U x, U x), for all x). Then this already
> suffices to show that U is unitary:
>
> We need to show (U x, U y) = (x, y) for all vectors x, y. So let x
and
> y be any non-zero vectors (otherwise the statement is trivial). Write
y
> = y_1 + y_2, where y_1 is parallel to x, and y_2 is perpendicular to
x
> (in any Hilbert space, this decomposition is unique), ie. y_1 = x
> (x,y)/(x,x), and y_2 = y - y_1. Then you have
>
> (U x, U y)
> = (U x, U (y_1 + y_2))
> = (U x, U y_1)
> = (x, y)/(x, x) (U x, U x)
> = (x, y)
>
> where the last equation follows because of Parseval's relation. So
> there you have it, U is unitary. This implies again, for example,
that
> U maps sets of orthonormal vectors to orthonormal vectors (ie. the
> angle between the vectors is also invariant under U) - something not
> immediately obvious from Parseval's relation! Neat, no? This of
course
> is the reason why fast correlation works (as the angle between two
> vectors is the same in both domains, thanks to unitarity of the DFT).

This is, incidentially, the types of arguments that form the
basis (no pun intended) for subspace methods. Such DSP techniques
are basically (eh...) a study of linear transforms.

> >  From a practitioner's standpoint, it allows one to make a sliding
> power
> > detecter by squaring each sample and boxcar filtering that sequence
> --
> > both operations have minimal cpu requirements independent of the
> length
> > of the window.  Pretty nifty.
>
> Yup. The whole of DSP is basically about being nifty :-).

Agreed.

Rune

"Rune Allnor" <allnor@tele.ntnu.no> wrote in message
news:1116009121.392786.87770@o13g2000cwo.googlegroups.com...

> Andor wrote:

> >                                  This implies again, for example,
> that
> > U maps sets of orthonormal vectors to orthonormal vectors (ie. the
> > angle between the vectors is also invariant under U) - something not
> > immediately obvious from Parseval's relation! Neat, no? This of
> course
> > is the reason why fast correlation works (as the angle between two
> > vectors is the same in both domains, thanks to unitarity of the DFT).

> This is, incidentially, the types of arguments that form the
> basis (no pun intended) for subspace methods. Such DSP techniques
> are basically (eh...) a study of linear transforms.

Since the two of you seem to be able to relate unitarity to
correlation, maybe one of you can make the connection for us
mere mortals.  I mean, unitarity seems to relate to how a
linear transformation transforms inner products, but the product
you do after transformation in the case of correlation is an
elemental product that yields another vector rather than an
inner product that yields a scalar.  And the unitarity of the
DFT is not really necessary for fast correlation to work --
the methods that lead to minimum arithmetic operation counts
don't directly use transformations that are strictly unitary.
Also the successful proof should demonstrate the well-known
fact that autocorrelation requires only about 3/4 the arithmetic
operations that correlation with a known constant vector does
whereas the opcounts for convolution and the convolution-square
are identical.

-- 
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end

ya know, James, i think i agree with you.  this whole thing is getting very
thick.  i think i know what the DFT is, what Parseval's Theorem is, what
unitary transforms are, what an inner product is, what fast convolution is,
but as i was watching this thread i was beginning to wonder similarly what
all this was about.  at least what this is becoming.  (so fill me in, too,
guys.)

i have no advice for the OP, so i'll just continue to watch.

r b-j


in article 6M6dnbPah-8chhjfRVn-gQ@comcast.com, James Van Buskirk at
not_valid@comcast.net wrote on 05/13/2005 17:26:

> "Rune Allnor" <allnor@tele.ntnu.no> wrote in message
> news:1116009121.392786.87770@o13g2000cwo.googlegroups.com...
> 
>> Andor wrote:
> 
>>> This implies again, for example, that
>>> U maps sets of orthonormal vectors to orthonormal vectors (ie. the
>>> angle between the vectors is also invariant under U) - something not
>>> immediately obvious from Parseval's relation! Neat, no? This of course
>>> is the reason why fast correlation works (as the angle between two
>>> vectors is the same in both domains, thanks to unitarity of the DFT).
> 
>> This is, incidentially, the types of arguments that form the
>> basis (no pun intended) for subspace methods. Such DSP techniques
>> are basically (eh...) a study of linear transforms.
> 
> Since the two of you seem to be able to relate unitarity to
> correlation, maybe one of you can make the connection for us
> mere mortals.  I mean, unitarity seems to relate to how a
> linear transformation transforms inner products, but the product
> you do after transformation in the case of correlation is an
> elemental product that yields another vector rather than an
> inner product that yields a scalar.  And the unitarity of the
> DFT is not really necessary for fast correlation to work --
> the methods that lead to minimum arithmetic operation counts
> don't directly use transformations that are strictly unitary.
> Also the successful proof should demonstrate the well-known
> fact that autocorrelation requires only about 3/4 the arithmetic
> operations that correlation with a known constant vector does
> whereas the opcounts for convolution and the convolution-square
> are identical.

-- 

r b-j                  rbj@audioimagination.com

"Imagination is more important than knowledge."

James Van Buskirk wrote:

> Since the two of you seem to be able to relate unitarity to
> correlation, maybe one of you can make the connection for us
> mere mortals.  I mean, unitarity seems to relate to how a
> linear transformation transforms inner products, but the product
> you do after transformation in the case of correlation is an
> elemental product that yields another vector rather than an
> inner product that yields a scalar.

(Circular) Convolution is just correlation with one of the sequences
inverted in time. The result is a vector with each coordinate the inner
product between two vectors, where one of them gets cyclically shifted.
Fast correlation is basically the same as fast convolution (modulo some
negative indexing and complex conjugate operation). Numerical Recipes
Chpt. 13.2 can fill in the details:

http://www.library.cornell.edu/nr/bookcpdf/c13-2.pdf


> And the unitarity of the DFT is not really necessary for fast
> correlation to work --

Unitarity is necessary. Correlation is the inner product between two
vectors. Unitarity implies that the inner product is invariant under
the transform -- that is why you can compute the inner product between
two vectors in both time and frequency domain and get the same result.

> the methods that lead to minimum arithmetic operation counts
> don't directly use transformations that are strictly unitary.

The FFT is the algorithm that speeds up fast correlation (I know you
know that). It computes a unitary transform.

Regards,
Andor

James Van Buskirk wrote:
> "Rune Allnor" <allnor@tele.ntnu.no> wrote in message
> news:1116009121.392786.87770@o13g2000cwo.googlegroups.com...
>
> > Andor wrote:
>
> > >                                  This implies again, for example,
> > that
> > > U maps sets of orthonormal vectors to orthonormal vectors (ie.
the
> > > angle between the vectors is also invariant under U) - something
not
> > > immediately obvious from Parseval's relation! Neat, no? This of
> > course
> > > is the reason why fast correlation works (as the angle between
two
> > > vectors is the same in both domains, thanks to unitarity of the
DFT).
>
> > This is, incidentially, the types of arguments that form the
> > basis (no pun intended) for subspace methods. Such DSP techniques
> > are basically (eh...) a study of linear transforms.
>
> Since the two of you seem to be able to relate unitarity to
> correlation, maybe one of you can make the connection for us
> mere mortals. I mean, unitarity seems to relate to how a
> linear transformation transforms inner products,

A "unitary tranform" is a "linear transform" that "preserves
the norm". That might sound like gobbledigook at first, but
let's define the terms and see if things might make some sense.

First, a "linear transform" is a general matrix-vector product
on the form

   Ax = b                                                   [1]

where a is a matrix and x and b are column vectors.

Now, there are two ways of viewing an expression like [1]. The
first is to see it as some system of equations, the other is to
view it as a geometrical operation on vectors.

If x is a 2D vector [x1, x2]^T, it can be plotted in a 2D
coordinate system. If A is of dimension 2x2, b will also
be a 2x1 vector that can be plotted in 2D.

Let's see an example: If

x = [1 1]^T                                                  [2]

A = | cos(t)  -sin(t) |                                      [3]
    | sin(t)   cos(t) |

then

b = [cos(t)-sin(t) sin(t)+cos(t)]^T                          [4]

and multiplication with A corresponds to rotating the
vector x an angle t around origo in the anticlockwise
direction.

Det(A) = cos^(t)+sin^2(t) = 1                                [5]

and so |x| = |b|, and multiplication with A does not change the
length of the vector. The matrix A "preserves the norm".

Undoing the effect of A (multiplying with inv(A)) amounts
to rotating an angle -t, and

inv(A) = | cos(t)  sin(t) | = A^T.                           [6]
         |-sin(t)  cos(t) |

A mathrix that preserves the norm and where

    inv(A) = conj(A^T)                                       [7]

is "unitary".

> but the product
> you do after transformation in the case of correlation is an
> elemental product that yields another vector rather than an
> inner product that yields a scalar.

Correlation is a sequence of inner products. To see this,
define a data matrix as

X = | x(0)  0    0  ...  0  |                                [8]
    | x(1) x(0)  0  ...  0  |
    | x(2) x(1) x(0)        |
    |                       |
    | x(N)     ....     x(0)|
    | 0   x(N) ...      x(1)|
    |                       |
    | 0    0            x(N)|

and the data vector

x = [x(0) x(1) ... x(N) 0 ... 0]^T                           [9]

The autocorrelation sequence becomes

r = [r(0) r(1) ... r(N)]^T = conj(X^T)x.                    [10]

So each coefficient in the correlation sequence r(k)
corresponds to the inner product between column k in the
matrix X and the vector x. Convolution can be defined
in a similar way.

> And the unitarity of the
> DFT is not really necessary for fast correlation to work --
> the methods that lead to minimum arithmetic operation counts
> don't directly use transformations that are strictly unitary.

This is a different question. The arithemtic algorithm that
yeilds the FFT has more to do with symmetries in the coefficients
than the DFT being unitary. In fact, in lots of implementations
one saves a number of arithmetic operations by skipping an 1/sqrt(N)
normalization factor. This is the reason for quite a few
"why doesn't Parseval's identity hold" questions here on
comp.dsp and comp.soft-sys.matlab.

> Also the successful proof should demonstrate the well-known
> fact that autocorrelation requires only about 3/4 the arithmetic
> operations that correlation with a known constant vector does
> whereas the opcounts for convolution and the convolution-square
> are identical.

Again, this has nothing to do with correlations being sequences
of inner products, but rather with properties of the vectors
in question.

Rune

"Andor" <an2or@mailcircuit.com> wrote in message
news:1116059723.682422.139790@o13g2000cwo.googlegroups.com...

> James Van Buskirk wrote:

> > And the unitarity of the DFT is not really necessary for fast
> > correlation to work --

> Unitarity is necessary. Correlation is the inner product between two
> vectors. Unitarity implies that the inner product is invariant under
> the transform -- that is why you can compute the inner product between
> two vectors in both time and frequency domain and get the same result.

Well, in at least the lowest operation count fast correlation
algorithms I know of, the linear transformation applied to
the inputs isn't even square, much less unitary.

> > the methods that lead to minimum arithmetic operation counts
> > don't directly use transformations that are strictly unitary.

> The FFT is the algorithm that speeds up fast correlation (I know you
> know that). It computes a unitary transform.

Well, the FFT is the algorithm most responsible for reducing
the operation count for large input vectors, but the sequence
of operations as described in NR is not going to reduce the
operation count to a minimum.  Fast correlation is the same
as fast convolution (except for autocorrelation, which is
fundamentally different) and the algorithms I have right at
the tips of my fingers compute convolution rather than
correlation, so let me present an example of fast cyclic
convolution on vectors of length 8.  For this length the
algorithm generated is the same as Nussbaumer, but hopefully
you can see that the NR sequence isn't going to get the
operation count this low:

! File: const8.i90
! Generated by conv2b.f90 05/14/2005 08:58:04.123
subroutine const8(h,a)
   implicit real(wp) (x,y)
   real(wp), intent(in) :: h(0:7)
   real(wp), intent(out) :: a(0:9)
   real(wp), parameter :: d1 = 0.2500000000000000000000000000000000_wp
   real(wp), parameter :: d2 = 0.5000000000000000000000000000000000_wp
   real(wp), parameter :: d3 = 0.7071067811865475244008443621048491_wp
   xr1_0 = h(0)+h(4)
   xr1_1 = h(1)+h(5)
   xr1_2 = h(2)+h(6)
   xr1_3 = h(3)+h(7)
   a(0) = d1*(xr1_0+xr1_2)
   a(1) = d1*(xr1_1+xr1_3)
   a(2) = d1*(xr1_0-xr1_2)
   a(3) = d1*(xr1_1-xr1_3)
   xr1_4 = h(0)-h(4)
   xr1_5 = h(1)-h(5)
   xr1_6 = h(2)-h(6)
   xr1_7 = h(3)-h(7)
   xr2_4 = d2*xr1_4
   xr2_5 = d2*xr1_5
   xr2_6 = d2*xr1_6
   xr2_7 = d2*xr1_7
   a(4) = xr2_4
   a(5) = xr2_5-xr2_4
   a(6) = xr2_7+xr2_4
   a(7) = xr2_6
   a(8) = xr2_7-xr2_6
   a(9) = xr2_5-xr2_6
end subroutine const8

! File: eval8.i90
! Generated by conv2b.f90 05/14/2005 08:58:04.139
subroutine eval8(h,a)
! 40 adds, 20 muls
   implicit real(wp) (x,y)
   real(wp), intent(inout) :: h(0:7)
   real(wp), intent(in) :: a(0:9)
   real(wp), parameter :: d1 = 0.7071067811865475244008443621048491_wp
   xr1_0 = h(0)+h(4)
   xr1_1 = h(1)+h(5)
   xr1_2 = h(2)+h(6)
   xr1_3 = h(3)+h(7)
   xr2_0 = xr1_0+xr1_2
   xr2_1 = xr1_1+xr1_3
   xr2_2 = xr1_0-xr1_2
   xr2_3 = xr1_1-xr1_3
   xr3_0 = a(0)*xr2_0+a(1)*xr2_1
   xr3_1 = a(1)*xr2_0+a(0)*xr2_1
   xr3_2 = a(2)*xr2_2-a(3)*xr2_3
   xr3_3 = a(3)*xr2_2+a(2)*xr2_3
   xr4_0 = xr3_0+xr3_2
   xr4_1 = xr3_1+xr3_3
   xr4_2 = xr3_0-xr3_2
   xr4_3 = xr3_1-xr3_3
   xr1_4 = h(0)-h(4)
   xr1_5 = h(1)-h(5)
   xr1_6 = h(2)-h(6)
   xr1_7 = h(3)-h(7)
   xr2_4 = xr1_4+xr1_5
   xr2_5 = xr1_6+xr1_7
   xr3_6 = a(4)*xr2_4-a(7)*xr2_5
   xr3_7 = a(7)*xr2_4+a(4)*xr2_5
   xr4_4 = xr3_6-a(6)*xr1_5-a(9)*xr1_7
   xr4_5 = xr3_6+a(5)*xr1_4-a(8)*xr1_6
   xr4_6 = xr3_7+a(9)*xr1_5-a(6)*xr1_7
   xr4_7 = xr3_7+a(8)*xr1_4+a(5)*xr1_6
   h(0) = xr4_0+xr4_4
   h(1) = xr4_1+xr4_5
   h(2) = xr4_2+xr4_6
   h(3) = xr4_3+xr4_7
   h(4) = xr4_0-xr4_4
   h(5) = xr4_1-xr4_5
   h(6) = xr4_2-xr4_6
   h(7) = xr4_3-xr4_7
end subroutine eval8

! File: eval8.i90
! Generated by conv2b.f90 05/14/2005 08:58:04.139
module mykinds
   implicit none
   integer, parameter :: dp = selected_real_kind(15,300)
end module mykinds

module test
   use mykinds, only: wp => dp
   implicit none
   private wp
   contains
include 'const8.i90'
include 'eval8.i90'
end module test

program conv_test
   use mykinds
   use test
   implicit none
   integer, parameter :: k = 3
   integer, parameter :: n = 2**k
   integer, parameter :: nta = n+(1+sign(1,n-4))/2*(n/2-2)
   real(dp) ta(0:nta-1)
   integer i
   integer j
   real(dp) h0(0:n-1)
   real(dp) h1(0:n-1)
   real(dp) error
   real(dp) max_error

   max_error = 0
   do i = 0, n-1
      h0 = 0
      h0(i) = 1
      call const8(h0,ta)
      error = 0
      do j = 0, n-1
         h0 = 0
         h0(j) = 1
         call eval8(h0,ta)
         h1 = 0
         h1(modulo(i+j,n)) = 1
         error = max(error,maxval(abs(h0-h1)))
      end do
      write(*,*) i,error
      max_error = max(error,max_error)
   end do
   write(*,*)
   write(*,*) 'max_error = ',max_error
end program conv_test

I was kind of surprised when I saw the generated code because
my generator

http://home.comcast.net/~kmbtib/conv2b.f90

was generating constants for an order-n DFT when it only
needs the constants for an order-n/2 DFT.  Well, if you
pay no attention to the little man in that particular booth,
hopefully you will be able to see some of the points I am
trying to make in this thread.

-- 
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end

James Van Buskirk schrieb:
> "Andor" <an2or@mailcircuit.com> wrote in message
> news:1116059723.682422.139790@o13g2000cwo.googlegroups.com...
>
> > James Van Buskirk wrote:
>
> > > And the unitarity of the DFT is not really necessary for fast
> > > correlation to work --
>
> > Unitarity is necessary. Correlation is the inner product between
two
> > vectors. Unitarity implies that the inner product is invariant
under
> > the transform -- that is why you can compute the inner product
between
> > two vectors in both time and frequency domain and get the same
result.
>
> Well, in at least the lowest operation count fast correlation
> algorithms I know of, the linear transformation applied to
> the inputs isn't even square, much less unitary.
>
> > > the methods that lead to minimum arithmetic operation counts
> > > don't directly use transformations that are strictly unitary.
>
> > The FFT is the algorithm that speeds up fast correlation (I know
you
> > know that). It computes a unitary transform.
>
> Well, the FFT is the algorithm most responsible for reducing
> the operation count for large input vectors, but the sequence
> of operations as described in NR is not going to reduce the
> operation count to a minimum. Fast correlation is the same
> as fast convolution (except for autocorrelation, which is
> fundamentally different) and the algorithms I have right at
> the tips of my fingers compute convolution rather than
> correlation, so let me present an example of fast cyclic
> convolution on vectors of length 8.  For this length the
> algorithm generated is the same as Nussbaumer, but hopefully
> you can see that the NR sequence isn't going to get the
> operation count this low:

James, if I am to believe Steven G. Johnson (and I have no reason no
to), then you are responsible for the current minimum FLOP count FFT
algorithm record, as described in http://www.fftw.org/newsplit.p=ADdf.
Congratulations!

I just wished you stop posting code and start talking a language I
understand (=3D maths). I have to invest a large amount of time just to
decipher you indices, never mind the rest of the code. Can't you
describe the general idea which makes your convolution (or FFT) faster
than everybody else's? Then we can discuss the (non-) unitarity of
these transforms.

Furthermore, while perhaps theoretically interesting, certain
reductions in FLOP have absolutley no impact on instruction counts on
specific hardwares (ie DSPs). An easy example is the symmetric FIR.
Theoretically, a N tap symmetric FIR can be computed with N/2
multiplies and N adds, making a total of 3/2 N FLOP (let's forget about
fast convolution for a moment). However, DSPs have a built in MAC
instruction that executes one multiply, one add and the necessary
memory transfers in one instruction, thus any N tap FIR can be computed
in N instructions. Therefore, this theoretical improvement gains you
exactly nothing on most DSP hardwares (in fact, a naive implementation
actually takes more instruction cycles).

Now I don't know if you or Steven have implemented your FFT codes on
machines where you explicitly take advantage of parallel computation
units. It would be really interesting to see the improvement in number
of execution cycles (which can be quite independent on the FLOPs count
of an algorithm) with respect to the standard DIT FFT (for example).

> ! File: const8.i90

What is that, Fortran?

Regards,
Andor

Andor wrote:
> Unitarity is necessary. Correlation is the inner product between two
> vectors. Unitarity implies that the inner product is invariant under
> the transform -- that is why you can compute the inner product
> between two vectors in both time and frequency domain and get
> the same result.

Andor, I think you've made a fundamental mistake here (besides the fact
that correlation is more than just a single inner product).  In order
for a fast convolution (or correlation) algorithm to work, you just
have to get the correct answer in time domain.  You do *not* have to
get the same answer in frequency domain (or whatever "domain" your
intermediate results correspond to, if anything), and thus unitarity
(or similar) is unnecessary.

Trivial counter-example: the unitary FFT would have a 1/sqrt(N) scale
factor, but this is actually the *wrong* normalization for a fast
convolution/correlation algorithm...you will end up with an extra
factor of 1/sqrt(N) at the end. You need to use a non-unitary
normalization where you have 1/N only once (typically, for the inverse
FFT, as in Matlab).

(There are a very large number of fast convolution algorithms out
there...I haven't checked this explicitly, but I suspect that many of
them, not just James', correspond to non-unitary transforms.)

Cordially,
Steven G. Johnson