Real and imaginary part of fourier transform using fft.

Hello.

Why am I getting a very different value from what is expected when we compute a fourier transform?

For example, let's consider a Gaussian distribution A*exp(-b^2*(y-y0).^2). Here, A is the amplitude, b is the parameter for controlling the width of the distribution and y0 is the point at which the distribution is centered.

I determined the fourier transform of this distribution analytically and obtained the transform as (A*sqrt(pi)/b)*exp(-(xaxis/(2*b)).^2).*exp(-1i*(xaxis*y0)).

I performed a simple test. I took y0 = 0. This would essentially give me the transform with only real part and no imaginary part. However, Matlab 'fft' gives a very different answer.

I am including the graphs obtained and the code so that I can paint a clear picture.

dy = 0.5;
y = dy*(-N/2+1:(N/2)).'; 
A = 52; % Amplitude 
b = 0.25; % Increase for decreasing the width of the distribution
y0 = 0; % Peak center
S = A*exp(-b^2*(y-y0).^2); % (Gaussian)
FT = fft(S)*dy; % Fourier transform
FT = ([FT(N/2+2:end,1);FT(1:floor(N/2+1),1)]);
df = (2*pi/(N*dy));
xaxis = df*(-N/2+1:(N/2)).'; 
FT_analytical = (A*sqrt(pi)/b)*exp(-(xaxis/(2*b)).^2).*exp(-1i*(xaxis*y0));

Can someone explain this to me, please?

Thanks.

Answer 1

[David Goodmanson]

Hi Venkata,

This is close, but you have to make a couple of adjustments. These are

y = dy*(-N/2:(N/2-1)).';
.
.
FT = fft(ifftshift(S))*dy; % Fourier transform

The first change puts y = 0 as the first point in the second half of the y array (assuming N is even). The second change swaps halves and moves the point corresponding to y = 0 down to the first point in the new y array. This is done because the fft always assumes that y=0 is the first point in the input array.

With those changes you get what you expect. Without the ifftshift, the fft is transforming an array that is shifted in the y coordinate, which produces an oscillatory phase factor in the result.

Also, if y corresponds to time and f is frequency, then the scaling is df = 1/(N*dy)

[Venkata Rajeshwar Majety]

Is it the same while computing ifft? Do we use ifft(fftshift())? Do we use it always? If not, when do we use it and when do we not use it?

I ask this because, I always prefer to write my own code so that I understand what's happening. For example, to calculate spectral density,

for ref = 1:17
  L = length(signal);
      t = (0:dt:(L-1)*dt).';
      T = t(end,1) + dt;
      tau = ((floor(-N/2+1:N/2).*dt)*10^4).';
      olap = 0.5; % 50 percent overlap
      trigger = olap*N;
      Nrecs = floor((L - N)/trigger + 1);
      G12 = zeros(floor(N/2+1),length(mic));S12 = zeros(N,length(mic));
      for count = 1:length(mic)
          w = zeros(N,Nrecs); p1 = w; p2 = w; Aw = zeros(Nrecs,1); X1 = p1;
          X2 = p2; Gxy = zeros(floor(N/2+1),Nrecs); Sxy = w; Sxy_pos = Gxy;
          Swxy = w; Gwxy = Gxy;

            for i = 1:Nrecs
                w(:,i) = 1 - cos(2*pi*t(((i-1)*trigger+1):((i-1)*trigger) + N,1)/(N*dt)); % Hann window
                p1(:,i) = (signal(((i-1)*trigger+1):((i-1)*trigger) + N,ref)).*(w(:,i));
                p2(:,i) = (signal(((i-1)*trigger+1):((i-1)*trigger) + N,count)).*(w(:,i));
                Aw(i,:) = 1/(mean(w(:,i).^2));
                X1(:,i) = fft(p1(:,i),N)*dt;
                X2(:,i) = fft(p2(:,i),N)*dt;
                Sxy(:,i) = (df).*conj(X1(:,i)).*X2(:,i);
                Sxy_pos(:,i) = Sxy(1:floor(N/2+1),i);
                Gxy(2:end-1,i) = 2*Sxy_pos(2:end-1,i);
                Gxy(1,i) = Sxy(1,i);
                Gxy(end,i) = Sxy(floor(N/2+1),i);
                Gwxy(:,i) = Gxy(:,i)*Aw(i,1);
                Swxy(:,i) = Sxy(:,i)*Aw(i,1);
            end
            G12(:,count) = mean(Gwxy,2);
            S12(:,count) = mean(Swxy,2);
    %       G12(:,count) = cpsd(signal(:,count),signal(:,ref),hann(N),[],N,Fs);
        end
  end

Is this correct? I checked it with last line of the code where it uses matlab's built in function to do that and it's almost similar except for the first few and last few points I suppose. when I use unwrap(angle()) to calculate the phase, it gives different answers. Now I wonder if I have been doing it all wrong?

Please clarify that to me. Thanks a lot.

[David Goodmanson]

I didn't check this in any detail, but yes, to do an fft the point corresponding to y = 0 should be the first point in input the array, and to do an ifft the point corresponding to f = 0 should be the first point in that input array. You just have to keep track of where things are. For example, if you do b = fft(a), b = fftshift(b) (say for plotting purposes), you need to do b = ifftshift(b) before using the ifft to go back. Either that or make a copy of b and fftshift that, leaving the original b alone to use in the ifft.