How does a Hermitian FFT work in Numpy?

Question

Say, I create a Hermitian complex signal using,

import numpy as np
t = np.arange(-4, 4)
z = np.exp(1j * t)

Here z should be a complex signal with Hermitian Symmetry, as you can see below.

In [2]: t
Out[2]: array([-4, -3, -2, -1,  0,  1,  2,  3])

In [3]: z
Out[3]:
array([-0.65364362+0.7568025j , -0.98999250-0.14112001j,
       -0.41614684-0.90929743j,  0.54030231-0.84147098j,
        1.00000000+0.j        ,  0.54030231+0.84147098j,
       -0.41614684+0.90929743j, -0.98999250+0.14112001j])

But when I take the Fourier Transform of this signal, I don't get a real spectrum,

In [8]: np.fft.fft(z)
Out[8]:
array([-1.38531768+0.7568025j, -7.02599565+0.7568025j,
        2.57935201+0.7568025j,  0.93952139+0.7568025j,
        0.41344309+0.7568025j,  0.08151870+0.7568025j,
       -0.22205190+0.7568025j, -0.60961892+0.7568025j])

However, when I take the Hermitian FFT using the hfft function, I get

In [10]: np.fft.hfft(z[:5])
Out[10]:
array([-1.38531768, -7.02599565,  2.57935201,  0.93952139,  0.41344309,
        0.0815187 , -0.2220519 , -0.60961892])

This is the real component of the result I got using the regular fft function. I don't understand what I am doing wrong here. Shouldn't fft be giving a result with a zero or almost-zero imaginary component when the input is Hermitian?

I feel like I am doing something wrong here regarding how the signal is presented to fft.

Answer 1

You are facing the practical problem of representing a noncausal signal in a programs such as Matlab, Octave, or Python. These programs, in their default, assume that all represented sequences begin at time-index $n=0$.

Therefore, a Hermitian symmetric sequence ( $x[n] = x^*[-n]$ for all $n$ ) , which begins at $n=-4$, and ends at $n=4$, will lose this symmetry when interpreted in those programs. The effect is a shift of the signal to the right by $4$ for your particular example.

Hence, what you have compouted is not the FFT of the original symmetric sequence $x[n] = e^{j n}$, but a shifted, non-symmetric, sequence $w[n] = x[n-4]$. Using the DTFT / DFT time-shift property, we see what is actually computed as : $$ W(e^{j\omega}) = e^{-j4\omega} X(e^{j\omega}) \tag{1} $$ and $$ W[k] = e^{-j4 \frac{2\pi}{N} k} X[k] ~~~~~ , ~~\text{ for } ~~~k=0,1,...,N-1 \tag{2}$$

Where $N$ is the length of the sequences $x[n]$ and its DFT $X[k]$.

Therefore, eventhough the DTFT / DFT of the desired sequence $x[n]$ has zero imaginary part, as an Hermitian symmetric signal would have, the practically computed $W[k]$ has nonzero imaginary part.

If you wish, you can obtain the non-zero phase samples of the desired DFT $X[k]$ of $x[n]$, from the computed DFT $W[k]$ of $w[n]$, with the following:

$$ X[k] = e^{j 4 \frac{2\pi}{N} k} W[k] ~~~~ ,~~\text{ for } k=0,1,...,N-1 \tag{3}$$

The following MATLAB code demonstrates the issue:

% Demonstrates the Hermitian symmetry and its handling in Matlab/Octave/Python
% environments...

a = -4;
b = +4;
n = [a:b];
x = exp(1j*n);      % Samples of Hermitian symmetrix input: x[k] = x[-k]*

figure,stem( imag( exp(-1j*2*pi*a*[0:(b-a)]/(b-a+1)).*fft(x) ));

The output is: