Why does an fftshift in the time domain give a different result compared to the frequency domain?

Question

I am trying to speed up some of my code and for physical reasons I need to do an FFT shift to get the FFT data in the standard DFT fourier representation. In Matlab this is performed in the frequency domain, with a data swap.

I had believed that I found a numerical trick, by performing the shift in the time domain, before transformation I was able to save a reasonable amount of computation time.

The problem is that for some reason the results are slightly different. I am curious as to why.

The Matlab code is as follows for a 1K signal

%1khz
f=1000;
t=0:.0001:5;
x=sin(2*pi*f*t);
%Wave
subplot(4,1,1);
plot(x)
%FFT
subplot(4,1,2)
ffted=fft(x);
plot(real(ffted));
%FFT Shifted
subplot(4,1,3)
shifted=fftshift(ffted);
plot(real(shifted));
%FFT Shifted by preshifting
subplot(4,1,4)
s=zeros(1,length(x));%In languages like C++ this is computationally inexpesnive
for i=1:length(x)
    if (bitget(i,1))%odd evens, works for real signal
        s(i)=-1;
    else
        s(i)=1;
    end
end
chirpshifted=fft(s.*x);
plot(real(chirpshifted));
%Difference
A=sum(real(chirpshifted));
B=sum(real(shifted));
A
B

Answer 1

There are a couple issues with your code. First, let's talk about why you can do what you're trying to do.

What you're trying to take advantage of is the duality between a multiplication in the time domain and convolution in the frequency domain. The fftshift() function simply swaps the two halves of a vector (which is assumed to be the result of a call to fft()). This operation is equivalent to a circular rotation of half of the vector length in the frequency domain.

One can accomplish the same operation by multiplying by a complex exponential function in the time domain before the FFT, as you noted. To get an $M$-sample shift, you need to multiply by a complex exponential that looks like:

$$ s[n] = e^{\frac{j2\pi nM}{N} } $$

For the case where $M = \frac{N}{2}$, this simplifies to:

$$ s[n] = e^{j\pi n} = [1, -1, 1, -1, \ldots] $$

This illustrates problem #1 with your script: it multiplies by the inverse of the above: $-s[n]$. You need to invert the sense of the if statement in the for loop to ensure that the vector s has the appropriate value.

Problem #2 stems from another simple MATLAB mistake: all of the vectors in your simulation are length 50001. So, the time-domain rotation that you're doing by multiplying by the $s[n]$ illustrated above actually results in an effective frequency domain shift of 25000.5 bins (instead of the 25000-bin shift implemented by fftshift()). The mismatch in results is due to this extra half-bin shift.

You can fix this problem in two ways:

1. Change the time-domain multiplication to multiply by an actual complex exponential vector at the appropriate frequency (i.e. set $M=25000$ and $N=50001$ in the first equation). The resulting signal at the input to the FFT will be complex.

2. Change your simulation to use even-length vectors (i.e. chop off the last element of the t vector that you initialize at the beginning).

I tweaked your simulation to fix the generation of s and truncate the vector length to 50000. It is as follows:

%1khz
f=1000;
t=0:.0001:5;
t=t(1:length(t)-1);
x=sin(2*pi*f*t);
%Wave
subplot(4,1,1);
plot(x)
%FFT
subplot(4,1,2)
ffted=fft(x);
plot(real(ffted));
%FFT Shifted
subplot(4,1,3)
shifted=fftshift(ffted);
plot(real(shifted));
%FFT Shifted by preshifting
subplot(4,1,4)
s=zeros(1,length(x));%In languages like C++ this is computationally inexpesnive
for i=1:length(x)
    if (~bitget(i,1))%odd evens, works for real signal
        s(i)=-1;
    else
        s(i)=1;
    end
end
chirpshifted=fft(s.*x);
plot(real(chirpshifted));

I then compare the results of the two methods:

>> diff = abs(chirpshifted-shifted);
>> sum(diff)
ans =

   3.3406e-21

As you can see, the results agree as perfectly as you would expect from floating-point math.