Circular convolution of real signal N by N/2 sized complex IFFT

Question

There are two related questions with good answers; for calculating the forward FFT in N/2 values, and calculating inverse FFT from N real values, now I'm in search for the missing link of packing the N/2+1 element complex signal X back to the N/2 element complex signal that it was derived in the first place for memory and arithmetic efficient inverse FFT.

Having calculated N/2+1 valued complex FFT from N real values (code from the first link):

x = [ 1 2 3 4 5 6 7 8 ]; n = 8; n2 = 4;
Z = fft(x(1:2:end) + i * x(2:2:end);
Ze = .5*( Z + conj([Z(1),Z(n2:-1:2)]));
Zo = -.5*j*( Z - conj([Z(1),Z(n2:-1:2)]));
X = [Ze,Ze(1)] + exp(-j*2*pi/n*(0:n2)).*[Zo,Zo(1)];

I have actually packed the real valued element X(n2+1) to imaginary part of X(1), saving some memory. Then I've point wise multiplied X with G representing the FFT of a same sized real signal g.

In order to finish the (circular) convolution of conv(x, g), I'd need to calculate the packed N/2 sized complex signal Z from (X.*G) for IFFT.

The first coefficient can be reverted as Z(1) = real(X(1)+X(n2+1))/2 + i*(imag(X(1)-X(n2+1)))/2;

As the other values X(2:n2-1) are originally produced from a bufferfly operation of x(k) +- exp(j*k/N) * x(n2-n), I would hope that all the elements would be recoverable. The question is how.

EDIT

I need to recheck the formulas, but since for a pair A=Z(k+1), B=Z(n2+1-k):

                                         [ 1-s 1+s   c     c ] [Real(A)]
 [Real(A') Imag(A') Real(B') Imag(B')] = [ 1+s 1-s  -c    -c ]*[Imag(A)]
                                         [ -c   c    1-s -1-s] [Real(B)]
                                         [ -c   c   -1-s  1-s] [Imag(B)]

with invertible matrix, there really is a solution.

Your code would be lot easier you to read if you add a few comments. — Hilmar, Apr 08 '20 at 13:56
Let me fully understand your question: Are you asking how to compute the circular convolution of two real sequences efficiently using the relationship $xcorr = ifft(fft(a)*fft(b))$ where a and b are real sequences and specifically fft(a) and fft(b) are computed using N/2 complex values as given by the first link you gave? Also in your case is N = 64 (It's confusing how the index 32 and 33 show up as it's no where else in the code you gave or your description)? — Dan Boschen, Apr 08 '20 at 14:55
@DanBoschen: Yes, that's what I'm after. And additionally I want to use the N/2 IFFT, which is just one element shorter than the X in the Matlab script. (My N might indeed be 64, but it's now edited, thanks for pointing out...) — Aki Suihkonen, Apr 08 '20 at 16:57
please learn to use $\LaTeX$ in your mathematical expression. — robert bristow-johnson, Apr 08 '20 at 19:02
just to clarify in my comment above it should be $xcorr = ifft(fft(a)*conj(fft(b))$ I was thinking they were real so showing the conjugate was not needed but the fft's are indeed complex! — Dan Boschen, Apr 09 '20 at 03:08

Dan Boschen · Accepted Answer · 2020-04-08T22:59:16.463

This is the general solution to recover N real samples using an N/2 Length IFFT. N/2+1 samples of the FFT are required for an even sequence N.

First generating a sample sequence of the N/2+1 samples of an FFT of a real signal:

x = [1 2 3 4 5 6 7 8];
N = length(x);
N2 = N/2;
fx = fft(x);
X = fx(1:N2+1);             % starting result: positive spectrum of real sequence

So above is the same result as can be alternatively computed with an N/2 FFT as already demonstrated.

Given such an result, below shows how to get the original sequence x using an N/2 IFFT. The code is self-explantory:

W = @(k,N) exp(-j*2*pi/N*[0:k-1]);                   % FFT Twiddle factors
Xe =  (X(1:N2) + conj(X(N2+1:-1:2)))/2;              % FFT of even samples
Xo = (X(1:N2) - conj(X(N2+1:-1:2)))./(2*W(N2,N));    % FFT of odd samples

newX = Xe + j*Xo;                      

xr = ifft(newX);                      % N/2 IFFT 
newxo =  real(xr)                     % recovered odd x samples
newxe =  imag(xr)                     % recovered even x samples

Since each FFT in your formula for cross correlation is the FFT of a real signal, the product will also be the FFT of a real signal, hence the technique outlined above can be used.

Circular convolution of real signal N by N/2 sized complex IFFT

1 Answers1