3
$\begingroup$

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
$\endgroup$
  • 1
    $\begingroup$ So length(x) is 50001? Does it work with an even length array, but not an odd length? Probably similar to dsp.stackexchange.com/q/10250/29 Remember that odd-length ffts are like [DC, 1, 2, 3, 3, 2, 1] while even-length are like [DC, 1, 2, 3, 4 (Nyquist), 3, 2, 1]. The DC component means it's not symmetrical in the obvious way, and odd-length FFTs don't include the Nyquist frequency. I think inverting every other sample only works for even-length arrays. $\endgroup$ – endolith Aug 20 '13 at 13:35
5
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.