Simultaneous ifft of two REAL signals

Started by Philip de Groot July 7, 2005
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"  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"  wrote in message 
news:wOudnT3I-Y_M8FDfRVn-rQ@centurytel.net...
> > "Philip de Groot" 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  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
Philip de Groot  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
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

"Fred Marshall"  wrote in message 
news:wOudnT3I-Y_M8FDfRVn-rQ@centurytel.net...
> > "Philip de Groot" 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
"Philip de Groot"  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

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 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

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
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