Understanding the DFT of a complex-valued time signal

209 views (last 30 days)
I am trying to understand how the fft function in MATLAB deals with the case of a complex signal used as input.
For a real input signal, I understand that the fft of a given signal results in a vector of Fourier coefficients that contains both positive and negative frequencies including the DC offset and the Nyquist frequencies. An example is shown below for a signal consisting of 2 cosine waves with different amplitudes, frequencies and phases:
%% Sampling Characteristics
fs = 1000; % Sample rate, Hz
N = 4000; % No. of samples taken in T
dt = 1/fs; % Time increment between samples,s
t = (0:(N-1))*dt; % Create time increment array from 0 to T
df = fs/N; % Frequency lines delF
%% Example of a REAL input signal
yy = 4*cos(2*pi*5*t + pi/6) + 2*cos(2*pi*60*t - pi/2); % Create REAL time signal
YY = fft(yy); % Calculates FFT with N "frequency bins"
fax_bins = 0:N-1; % Array of "frequency bins"
fax_Hz = fax_bins*df; % Corresponding frequencies, Hz
YY_norm = abs(YY)/N; % Normalise
YY_norm_spec = YY_norm(:,1:N/2); % Take positive freqs only
YY_norm_spec(:,2:end) = 2*YY_norm_spec(:,2:end); % Calculate 1-sided amplitude spectrum, don't multiply DC
figure (1)
subplot(3,1,1); plot(t,yy,'k') % Plot original time series
title('Original Time Series')
xlabel('Time (s)')
ylabel('Amplitude')
grid on
hold on
subplot(3,1,2); plot(fax_Hz,YY_norm,'.-') % Plot raw output from FFT
title('Raw FFT Output')
xlabel('"Frequency bins" or FFT points')
ylabel('Amplitude')
hold on
grid on
subplot(3,1,3); plot(fax_Hz(1:N/2),YY_norm_spec(1:N/2),'.-') % Plot amplitude spectrum
title('1-sided Amplitude Spectrum')
xlabel('Frequency (Hz)')
ylabel('Amplitude')
grid on
Resulting in:
Which makes sense as the raw FFT output is reflected about the Nyquist frequency with half the amplitude of the two original signals coming from the positive frequencies and half the amplitude coming from the negative frequencies. This can be transformed into a 1-sided amplitude spectrum by neglecting the negative frequencies and doubling the amplitude of the positive frequencies.
However, I am struggling to understand what happens when you start with a complex input signal. For example:
%% Example of a COMPLEX input signal
zz = 4*exp(1i*(2*pi*5*t + pi/6)) + 2*exp(1i*(2*pi*60*t - pi/2); % Complex rep. of two signals
p_spectrum = fft(zz); % Calculates fft in N "frequency bins"
P2 = abs(p_spectrum/N); % Normalise amplitudes by N "frequency bins"
P1 = P2(:,1:N/2+1); % Take 1st half of frequencies (i.e. positive freqs) and the N/2+1 index which is Nyquist
P1(:,2:end-1) = 2*P1(:,2:end-1); % Multiply all amplitudes by 2 except DC (index 1) and Nyquist(N/2+1)
figure (2)
subplot(3,1,1); plot(t,real(zz),'k') % Plot original time series
title('Original Time Series')
xlabel('Time (s)')
ylabel('Amplitude')
grid on
hold on
subplot(3,1,2); plot(fax_Hz,P2,'.-') % Plot output from FFT after normalisation
title('FFT Output after Normalisation')
xlabel('"Frequency bins" or FFT points')
ylabel('Amplitude')
hold on
grid on
subplot(3,1,3); plot(fax_Hz(1:N/2),P1(1:N/2),'.-') % Plot amplitude spectrum
title('1-sided Amplitude Spectrum')
xlabel('Frequency (Hz)')
ylabel('Amplitude')
grid on
Why in the case of a complex signal are there no negative frequencies and hence the amplitude is double what is expected? (if following same process as the real example)
My guess would be that it is due to Euler's formula:
hence only pos. frequency at amplitude A. Whereas for a real signal you have pos & neg at amplitude A/2:
In which case you only need to discard the latter half of the points from the fft output to achieve the correct spectrum?
Thanks in advance!
  1 Comment
Paul
Paul on 19 Nov 2023
This example illustrates why the concept of the "one-sided spectrum" doesn't apply for complex-valued signals.

Sign in to comment.

Answers (2)

Mehmet Selim Akay
Mehmet Selim Akay on 18 Nov 2023
Try on an example: A rotor tip has horizontal (x) and vertical (y) position signals each made up of a DC component, a forward whirl frequency and a backward whirl frequency (+ and - as written inside the exponents or xr). This motion is well expressed as xr, which gives the following x and y real signals shown. Then, let's combine the two signals together in a complex form as shown in xc (Note that xc is strictly different from xr's non-real() part because whatever you do with xr, it is still two distinct signals.)
% Define the signal
w1 = 5 ; w2 = 2*w1;%SIGNAL: 5 rad/s for FW, 10 rad/s for BW.
w_signal = gcd(w1,w2) ;%Greated common divider: Freq of the whole signal
T = 2*pi/w_signal * 10 ;%10 full periods is the total duration of the signal.
Nd = 1000 ;%Time series sample count (1000 samples for each period)
fs = Nd/T ;%Sampling rate, Nd samples in T sec.
Nh = Nd ; %=5 or whatever. Number of harmonics sought. Nh includes zero here.
t = (0:Nd-1)/Nd*T ;
wS = 2*pi*(-Nh/2:Nh/2-1)*fs/Nh ;
A = 1 ; B = 2 ; a = 3 ; b = 3; d = 5 ; e = 5; phi1 = 2*pi/30 ; phi2 = 0 ;
xr = real( [A;B] + [a ;-1j*b] * exp(+1j*(w1*t - phi1)) + [d ; -1j*e] * exp(-1j*(w2*t-phi2)) ) ;
x = A + a*cos(w1*t - phi1) + d*cos(w2*t - phi2) ;% xr(1,:)
y = B + b*sin(w1*t - phi1) - e*sin(w2*t - phi2) ;% xr(2,:)
xc = x - 1j*y ;
% DFT
for k = 1:Nh
Xc(k) = 1/Nd*dot(xc, exp(-1j*wS(k)*t) ) ;%Note the conjugation: Xc(k)=sum(conj(xc).*exp(...))
end
figure(100) ; plot(wS, abs(Xc)) ; title("abs(Xc)")
The whole goal of this complex signal form (xc=x-jy) is to get the FW freqs positive after the DFT, and vice versa (BW freqs as negative).
If you make the rotor orbit FW and BW components elliptic, because this actually includes other frequencies, you have a different plot:
A = 1 ; B = 2 ; a = 3 ; b = 4; d = 6 ; e = 8; phi1 = 2*pi/30 ; phi2 = 0 ;
x = A + a*cos(w1*t - phi1) + d*cos(w2*t - phi2) ;
y = B + b*sin(w1*t - phi1) - e*sin(w2*t - phi2) ;
xc = x - 1j*y ;

William Rose
William Rose on 8 Dec 2021
@Martin Doherty, the FFT of your complex sequence is zero for negative frequencies (i.e. for frequencies above the Nyquist frequency) because your signal is fortuitously simple. Consider another example:
N=100; dt=1/N; T=N*dt; df=1/(N*dt); t=dt*(0:N-1); f=df*(0:N-1);
a1=2+2*i; a2=1-i; %complex amplitudes
f1=9; phi1=pi/6; f2=17; phi2=4*pi/3;
x1=a1*cos(2*pi*f1*t+phi1)+a2*cos(2*pi*f2*t+phi2);
X1=fft(x1);
figure; subplot(2,1,1); plot(f,real(X1)); ylabel('Real(fft)'); grid on
subplot(2,1,2); plot(f,imag(X1)); ylabel('Imag(fft)'); xlabel('Frequency (Hz)'); grid on
As you see from the plots, the FFT of this complex signal has power at frequencies above and below the Nyquist frequency.
  1 Comment
William Rose
William Rose on 9 Dec 2021
@Martin Doherty, a general complex sinusoid
is represented in the discrete Fourier transform by
where a+, b+ are the real, imaginary parts of the positive spectral component, and a-, b- are the real, imginary parts of the negative spectral component. You solve for by applying Euler's formula and equating like terms. You get:
When x(t) is real, B=D=0, and in that case, the equations above simplify to
and
.
In other words, for a real signal, the positive and negative spectral components are complex conjugates.
For a complex signal, the positive and negative spectral components are NOT complex conjugates.
You happened to pick a special complex signal whose negative components are zero. Your signal is of the form
which equals , if you set
and . When you plug those values for A, B, C, D into the earlier equations, you get the following FFT components:
, , . In other words, the negative frequency part is zero.

Sign in to comment.

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!