Simultaneous ifft of two REAL signals

Hello,

I am playing a little bit with audio signals and Fourier analysis. Since 
audio is a real 2-channel thing, I used the twofft algorithm to 
simultaneously obtain both left and right channel Fourier transforms, see:

http://www.library.cornell.edu/nr/bookcpdf/c12-3.pdf

However, I cannot manage to simultaneously perform the inverse Fourier. It 
should be very easy, according to the refence link mentioned above:

quote: "use the fact that the FFT is linear and form the sum of the first 
transform plus i times the second."

Sounds easy enough, but how am I supposed to do this (read: what do they 
actually mean in C++ code)?

Any help is highly appreciated!

Regards,

Philip

Philip de Groot wrote:
> Hello,
>
> I am playing a little bit with audio signals and Fourier analysis. Since
> audio is a real 2-channel thing, I used the twofft algorithm to
> simultaneously obtain both left and right channel Fourier transforms, see:
>
> http://www.library.cornell.edu/nr/bookcpdf/c12-3.pdf
>
> However, I cannot manage to simultaneously perform the inverse Fourier. It
> should be very easy, according to the refence link mentioned above:
>
> quote: "use the fact that the FFT is linear and form the sum of the first
> transform plus i times the second."
>
> Sounds easy enough, but how am I supposed to do this (read: what do they
> actually mean in C++ code)?
>
> Any help is highly appreciated!

Well, I'm not sure this is "easy". First, we need to review some of the

basic properties of the Fouriert Transform (FT). The FT is linear,
so the quoted statement is partially correct. The problem, then, is to
find a way to compute the IFT so that channel 1 ends up in the real
part
of the complex-valued time series, while channel 2 ends up in the
imaginary part.

One could use such properties as the real part corresponding to the
even component of the spectrum and the imaginary part corresponding
to the odd part of the spectrum, but those are "artifical" properties
in the sense that your time signals are not even and odd.

We also note theat the spectrum of a real-valued sequence is conjugate
symmetric, F(-w) =  conj(F(w)), so there is 100% redundancy in the
spectrum. This is a property that ought to be exploited.

One can find the real-valued time series from the one-sided FT as

  x(t) = 2*real(IFT{X(w+)})

where X(w+) is the spectrum for positive w only. The trick, then,
is to find an expression in terms of X(w-) that gives a purely
imaginary signal x'(t), and then use the coefficents between N/2
and N-1 in the spectrum, to store these coefficients.

I am in too much of a vacation mode to do that, but I suspect
it would involve the spectrum with +j or -j, and possibly taking
a conjugate before or after that.

So, the procedure then becomes something like (in matlab code)

% X1: Spectrum of signal 1 (N points, N even)
% X2: Spectrum of signal 2 (N points, N even)

X= zeros(N,1);
X(1:N/2)= X1(1:N/2);
X(N/2+1:N)=conj(sqrt(-1)*X2(N/2+1:N)); % DON'T BELIEVE THE
                                       % DETAILS HERE!!!!

x=ifft(X);
x1=2*real(x);
x2=2*imag(x);

Again, my brain is not sharpened enough to help you out with all
the details, but I hope the general procedure outline is of some
help to you.

Rune

Philip wrote:

>However, I cannot manage to simultaneously perform the inverse Fourier. It
>should be very easy, according to the refence link mentioned above:
>
>quote: "use the fact that the FFT is linear and form the sum of the first
>transform plus i times the second."

You have two vectors, x1 and x2, with discrete Fourier Transforms X1
and X2 respectively. Then

x1 = Re[ IDFT[ X1 + i X2 ] ]

and

x2 = Im[ IDFT[ X1 + i X2 ] ].

The reason for this is the linearity of the DFT / IDFT and the fact
that a purely real vector gets transformed into a hermitian symmetric
vector, whereas a purely imaginary vector gets transfromed into an
anti-hermitian symmetric vector.

Regards,
Andor

Andor wrote:
> Philip wrote:
>
> >However, I cannot manage to simultaneously perform the inverse Fourier. It
> >should be very easy, according to the refence link mentioned above:
> >
> >quote: "use the fact that the FFT is linear and form the sum of the first
> >transform plus i times the second."
>
> You have two vectors, x1 and x2, with discrete Fourier Transforms X1
> and X2 respectively. Then
>
> x1 = Re[ IDFT[ X1 + i X2 ] ]
>
> and
>
> x2 = Im[ IDFT[ X1 + i X2 ] ].

Ouch! So ridiculously simple, obvious and straight-forward once
you see it. My brain have probably boiled, it's too hot for comfort
here. You people in more tempereate/tropical areas will probably
laugh in disguste when I end up in trouble from 25C for four days
straight...

I *did* spend an entire summer in Italy a few years back, but the
Italians are sensible enough to build houses that are cool inside
in summer time. Norwegian houses are built to stay warm inside
during cold winters. The design works, even in summertime...

> The reason for this is the linearity of the DFT / IDFT and the fact
> that a purely real vector gets transformed into a hermitian symmetric
> vector, whereas a purely imaginary vector gets transfromed into an
> anti-hermitian symmetric vector.

Perhaps you are right. I prefer to use the argument that the DFT is
linear and that a purely real and a purely imaginary vector can
not interact, regardless of any symmetries.

Rune

"Philip de Groot" <groot877@zonnet.nl> wrote in message 
news:Xns968C9F88E9EBDgroot877zonnetnl@137.224.11.5...
> Hello,
>
> I am playing a little bit with audio signals and Fourier analysis. Since
> audio is a real 2-channel thing, I used the twofft algorithm to
> simultaneously obtain both left and right channel Fourier transforms, see:
>
> http://www.library.cornell.edu/nr/bookcpdf/c12-3.pdf
>
> However, I cannot manage to simultaneously perform the inverse Fourier. It
> should be very easy, according to the refence link mentioned above:
>
> quote: "use the fact that the FFT is linear and form the sum of the first
> transform plus i times the second."
>
> Sounds easy enough, but how am I supposed to do this (read: what do they
> actually mean in C++ code)?
>
> Any help is highly appreciated!
>

Philip,

Andor's response is very nice.  However, I don't think it addresses what 
*seems* to be part of your question:

While you can pack purely real sequences by 2 in the forward transform, the 
result is two complex sequences or sequences with 4 samples per real input 
sample.  I don't see any way to pack two complex sequences to do the inverse 
transform.  So, you have to do *2* complex inverse transforms as Andor 
shows.

Does that help?

Fred 

"Fred Marshall" <fmarshallx@remove_the_x.acm.org> wrote in message 
news:wOudnT3I-Y_M8FDfRVn-rQ@centurytel.net...
>
> "Philip de Groot" <groot877@zonnet.nl> wrote in message 
> news:Xns968C9F88E9EBDgroot877zonnetnl@137.224.11.5...
>> Hello,
>>
>> I am playing a little bit with audio signals and Fourier analysis. Since
>> audio is a real 2-channel thing, I used the twofft algorithm to
>> simultaneously obtain both left and right channel Fourier transforms, 
>> see:
>>
>> http://www.library.cornell.edu/nr/bookcpdf/c12-3.pdf
>>
>> However, I cannot manage to simultaneously perform the inverse Fourier. 
>> It
>> should be very easy, according to the refence link mentioned above:
>>
>> quote: "use the fact that the FFT is linear and form the sum of the first
>> transform plus i times the second."
>>
>> Sounds easy enough, but how am I supposed to do this (read: what do they
>> actually mean in C++ code)?
>>
>> Any help is highly appreciated!
>>
>
> Philip,
>
> Andor's response is very nice.  However, I don't think it addresses what 
> *seems* to be part of your question:
>
> While you can pack purely real sequences by 2 in the forward transform, 
> the result is two complex sequences or sequences with 4 samples per real 
> input sample.

..... oops .. 2 samples per real input sample.  2 real N-length sequences 
in, 2 complex N-length sequences out.

Fred 

Fred, you're forgetting the Hermitian symmetry of the outputs; you can
indeed do the inverse in one complex transform because of this
redundancy.

All of these FFT algorithms are trivially invertible with the same
number of operations as the forward transform, since by unitarity the
inverse is just the conjugate transpose (i.e. do all the steps
backwards, in the sense of a linear network where you reverse the
edges, with conjugated constants).

Cordially,
Steven G. Johnson

Philip de Groot <groot877@zonnet.nl> wrote in
news:Xns968C9F88E9EBDgroot877zonnetnl@137.224.11.5: 

> Hello,
> 
> I am playing a little bit with audio signals and Fourier analysis.
> Since audio is a real 2-channel thing, I used the twofft algorithm to 
> simultaneously obtain both left and right channel Fourier transforms,
> see: 
> 
> http://www.library.cornell.edu/nr/bookcpdf/c12-3.pdf
> 
> However, I cannot manage to simultaneously perform the inverse
> Fourier. It should be very easy, according to the refence link
> mentioned above: 
> 
> quote: "use the fact that the FFT is linear and form the sum of the
> first transform plus i times the second."
> 
> Sounds easy enough, but how am I supposed to do this (read: what do
> they actually mean in C++ code)?
> 
> Any help is highly appreciated!
> 
> Regards,
> 
> Philip
> 

Hello,

Thank all of you for your help. The answer of Andor is in agreement with 
what is stated in the textbook, but it is also something I have tried to 
do (unsuccesfully):

>You have two vectors, x1 and x2, with discrete Fourier Transforms X1
>and X2 respectively. Then
>
>x1 = Re[ IDFT[ X1 + i X2 ] ]
>
>and
>
>x2 = Im[ IDFT[ X1 + i X2 ] ].

My problem is: how to do this in C++? The answer of Rune gave me the 
bright idea to test something out in Matlab: what happens if a complex 
number is multiplied with i. This way, I found out that:

      i(real+i*imag) = -imag + i*real.

When coding in C++, I did not realise this (stupid of course, but 
afterwards it's always easy)... I know what the problem is now and 
hopefully I will get it right soon (and enjoy a speed improvement of 
100%).

Note that I posted a second question to this mailing list regarding FIR 
filtering in the Fourier domain. If you want to help me out with that 
one too?

Regards,

Philip