Suppose, 

I have a 1000 seconds normally distributed time domain samples sampled at 500KS/s. Take fourier transform using window of 1 sec. Now, I have 1000 sets of Fourier coefficients in which each set has frequency bins starting from 1Hz to 250,000Hz. If I will plot the magnitude distribution of certain frequency let's say 1000Hz (magnitude of 1000th bin from all 1000 sets = 1000 samples), then it should be Rayleigh distribution according to the theory. 

If I will separately plot the amplitude distribution of real and imaginary parts of the Fourier coefficients of 1000 Hz, then it comes out to be normally distributed. This observation is based on the simulation. 

Can anyone tell me if it theoretically correct and explain the intuition behind it? 

For ANY bin in the frequency domain the FFT coefficient will be linear combination of gaussian.So real and imaginary part will be gaussian.However the magnitude which is +ve square root of squared sum (All -ve values of gaussian will become +ve due to squaring).

The square root will be Rayleigh and phase will be uniform.

I see Rayleigh in my simulation:

%%%%%%%%%%%%%%%%%%%%%%%

clear all;
n = 10000;
f = zeros(n,1024);
for i = 1:n
   x=randn(1,1024);
   f(i,:) = abs(fft(x));
end

test = round(f(:,97)); %test this bin

for y = 0:max(test)
   h(y+1) = length(find(test == y));
end

plot(h);

%%%%%%%%%%%%%%%%%%%%%

Hi Amartansh,

Interesting observation.

I think you will find the answer here:

https://en.wikipedia.org/wiki/Rayleigh_distributio...

Mr C