How does a Hermitian FFT work in Numpy?

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 : t
Out: array([-4, -3, -2, -1,  0,  1,  2,  3])

In : z
Out:
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 : np.fft.fft(z)
Out:
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 : np.fft.hfft(z[:5])
Out:
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.

You are facing the practical problem of representing a noncausal signal in a computing environment such as the Matlab/Octave or Python. These systems, in their default, assume that all represented signals begin from their time index $n=0$.

Therefore a signal which begins at the index $n=-4$ and ends at $n=4$ that has the hermitian (complex-conjugate) symmetry such as $z[n] = z[-n]^*$ will loose this symmetry when represented in those environments. The effect is equivalent to a shift of the signal to the right by $4$ for this particular example.

Hence what you take the FFT of is not the original symmetric signal $z[n] = e^{j n}$ bu the following shifted non-symmetric version $w[n] = z[n-4]$. From the DTFT/DFT properties it can be shown that what is actually computed is: $$W(e^{j\omega}) = e^{-j4\omega} X(e^{j\omega})$$ and $$W[k] = e^{-j4 \frac{2\pi}{N} k} X[k] ,~~\text{ for } k=0,1,...,N-1$$

Where $N$ is the length of the signal $z[n]$ and the size of the DFT $X[k]$.

Therefore, eventhough the theoretical signal $X[k]$ has zero imaginary part as its input was a hermitian symmetric signal, the practically computed DFT $W[k]$ has nonzero imaginary part.

Finally you can obtain the theoretical DFT $X[k]$ with the inverse shift operation, applied in frequency domain:

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

The following MATLAB code demonstrated the issue:

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

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

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

The output is: 