Suppose we have two real signals: x(t) and y(t) and we need to compute their respective spectra:

X(f)=FFT[x(t)]

and

Y(f)=FFT[y(t)]

This is possible by using just one (instead of two) FFT operations provided that the FFT algorithm used is capable to process a complex signal. There is a little extra cost in computation but it is possible. This “recipe” explains how:

a) build the complex signal: z(t) = x(t) + j * y(t)

b) computing the complex FFT of z(t) gives:

FFT[z(t)] = Z(f) = X(f) + j * Y(f) = Re[X(f)] + j * Im[X(f)] + j * [ Re[Y(f)] + j * Im[Y(f)] ]
FFT[z(t)] = Re[X(f)] - Im[Y(f)] + j * ( Im[X(f)] + Re[Y(f)] )

c) note the real and imaginary parts of FFT[z(t)] above and note also that for a real signal sampled at fs sampling frequency the following obvious relations hold:

Re[X(f)] = Re[X(fs - f)]
Im[X(f)] = - Im[X(fs - f)]

d) then it follows:

Re[Z(f)] = Re[X(f)] - Im[Y(f)]
Re[Z(fs - f)] = Re[X(fs - f)] - Im[Y(fs - f)] = Re[X(f)] + Im[Y(f)]
Im[Z(f)] = Im[X(f)] - Re[Y(f)]
Im[Z(fs - f)] = Im[X(fs - f)] - Re[Y(fs - f)] = -Im[X(f)] + Re[Y(f)]

e) and at the end it is possible to recover the two spectra as:

Re[X(f)] = (Re[Z(f)] + Re[Z(fs - f)]) / 2
Im[X(f)] = (Im[Z(f)] - Im[Z(fs - f)]) / 2
Re[Y(f)] = (Im[Z(f)] + Im[Z(fs - f)]) / 2
Im[Y(f)] = (Re[Z(fs - f)] - Re[Z(f)]) / 2

… done !

Curious to test the recipe?

USE THESE MATLAB LINES

% number of samples
N=ns;
% windowing
wstr=@blackmanharris;
nw=ww; % window width
wn = window(wstr,nw);
w=zeros(size(N));
w(fix((N-nw)/2)+(1:nw))=wn;
x=x.*w; % signal1 windowed
y=y.*w; % signal2 windowed
z=complex(x,y); % “complex” signal synthetized
Z = fft(z,N); % complex fft
% demultiplexing
a=zeros(size(N));
b=zeros(size(N));
a(2:N/2)=(real(Z(2:N/2))+real(fliplr(Z(N/2+2:N))))/2;
b(2:N/2)=(imag(Z(2:N/2))-imag(fliplr(Z(N/2+2:N))))/2;
X=complex(a,b);
a(2:N/2)=(imag(Z(2:N/2))+imag(fliplr(Z(N/2+2:N))))/2;
b(2:N/2)=(real(fliplr(Z(N/2+2:N)))-real(Z(2:N/2)))/2;
Y=complex(a,b);

COMMENTS:
According to MATLAB documentation: “For most N, real-input DFTs require roughly half the computation time of complex-input DFTs. However, when N has large prime factors, there is little or no speed difference. The execution time for fft depends on the length of the transform. It is fastest for powers of two. It is almost as fast for lengths that have only small prime factors. It is typically several times slower for lengths that are prime or which have large prime factors.”
This means that this recipe is of little or no value for using with MATLAB but in the case you need many FFTs of large prime factors length.
But if you need to do lots of FFT calculations with not optimized FFT algorithms this trick is very Good2Know. I first used this trick in 1981 and I was able to compute a 512 points FFT coded in Pascal on my old and faithfull Apple][ computer at the "breath-taking" speed of 0.022 FFTs/sec (or one single FFT in 45 seconds... it was just a 64 Kbytes Apple ][ !!!)

Have fun!
Stadero