Aeraging complex FFT to keep phase information

Hello, My question might seem naive but i can't find a clear answer to it so i'm asking for you opinion :

I am trying to find the correct method to average a complex Spectrum calculated overs sevreals Windows of a signal keeping the phase information. So the result i want is a complex averaged Spectrum, so that i can manipulate the phase of my signal in post processing.

Here is my startegy for now.

- I store the resulst of my FFTs in a matrix (4096*10 for a FFT size of 4096 points and 10 Windows)

- For each frequency, i avergae separatly the real and imaginary parts of my spectrums; Then i ake the angle of the result (sum real +i*imag)

- For each frequency, i average the absolute value to have the amplitude.

- Finally; for each frequency, i construc a complex assembling the two previous results.

The final result is a 4096 constructed complex Spectrum.

Is my startegy correct ? I can't help but think there is someting wrong ine the way i'm averaging complex numbers. Plese note that i have no phase reference, all my signal are real signals, recorded from an experimental set-up.

Here is the part of my code i'm talking about :

A = zeros(Nfft,1);
      for i=1:Nfft
                    %A(i)=mean(FFT(i,:));
          RES=mean(real(FFT(i,:)))+1i*mean(imag(FFT(i,:)));
          A(i)=mean(2.*2.*abs(FFT(i,:))./Nfft).*exp(1i*angle(RES));
                     %A(i)=mean(real(FFT(i,:)))+1i*mean(imag(FFT(i,:)));
      end