1
$\begingroup$

I have been looking around on the internet for a really clear explanation for why fft and fftshift give different results.

I have a very long time series for which I have applied a Hanning window. The sample rate is $8\textrm{ Hz}$ and the recording is $\approx 9$ hours long ($262144$ point time series). The time series is assumed to be stationary.

When I do an fft in MATLAB, I get the following result:

enter image description here

When I do an fftshift in MATLAB, I get the following result:

enter image description here

Four questions:

  1. How do I decide which one to use?
  2. Which is "correct"
  3. Why are they giving such different answers?
  4. Why is the fftshift only real-valued while the fft is complex?

Pseudo-code:

t=load('t.mat'); %load time series
N = length(t);
dt = 0.125; Fs = 1/dt; df = Fs/N;
f = -Fs/2:df:Fs/2-df;  

t_hann = hanning(t);
F_fft = fft(t);
F_fftshift = fftshift(t);
figure(1); plot(f,F_fft);
figure(2); plot(f,F_fftshift);

Any help is greatly appreciated.

Cheers

$\endgroup$
  • $\begingroup$ could you show a snippet of your code? like where hanning() is applied and where fft() and fftshift() are applied. in what order? and how are you getting the limits of -4 to 4 on the horizontal axis of your plot? $\endgroup$ – robert bristow-johnson Apr 20 '17 at 18:02
  • 2
    $\begingroup$ fftshift is not performing FFT, it's merely shifting the output of the fft command so that it becomes symmetric around 0 Hz. All of that is explained in the MATLAB docs. Unless I don't understand your question. $\endgroup$ – jojek Apr 20 '17 at 18:02
  • $\begingroup$ I'm voting to close this question as off-topic because this is the RTFM question. If you have any doubts about it, please edit the question. $\endgroup$ – jojek Apr 20 '17 at 18:03
  • $\begingroup$ dunno "RTFM"... $\endgroup$ – robert bristow-johnson Apr 20 '17 at 18:04
  • $\begingroup$ @robertbristow-johnson: ;) $\endgroup$ – jojek Apr 20 '17 at 18:06
5
$\begingroup$

I believe you simply misunderstand what these functions are supposed to do. You get different results because they are not meant to do the same thing. I think it would be relevant that you read the documentation of these 2 functions.

Answers below:

  1. How do I decide which one to use?

    • fft : Use it when you need to compute the Discrete Fourier Transform of your time series
    • fftshift : Use it when you want to to swap the position of the negative and positive frequencies of a spectrum gotten with fft()
  2. Which is "correct"?

    What do you consider correct ? fft will be correct if you need a frequency spectrum of your time domain signal. It doesn't really make sense to apply fftshift on a time domain series.

  3. Why are they giving such different answers?

    Because they don't do the same thing.

  4. Why is the fftshift only real-valued while the fft is complex?

    fft computes the discrete Fourier transform and by definition the output is complex. fftshift doesn't compute anything except swaping the position of the samples, so if your input is real, you get real output.

Normally we use these two functions this way : spectrum = fftshift(fft(x))

$\endgroup$
4
$\begingroup$

all fftshift() does is swap the first half and the second half of the array supplied to it. the reason for doing so have several, mostly historical and conventional, roots.

  1. for some dumb reason, the creator of MATLAB (a nice guy named Cleve) was never able to grok that the fixed array indices of each array is a problem to be solved by allowing the MATLAB user to change the array limits. in MATLAB, all arrays have an origin of 1. so if you have an N-length array, X, it always begins with X(1) and ends with X(N). so if those limits, 1 to N, are not the natural limits of time or frequency of your array, the MATLAB solution to deal with that is for the user to have to create and maintain a separate array, of the same size N, that has for the array values the indices that have the natural range the user wants.

it appears to me that the user must have created such an array for the plot() function to yield a horizontal axis from -4 to +4. now if (to keep thinks simpler) we leave out the df factor (call it 1), the f array you have goes from -N/2 to N/2-1 rather than from 1 to N. either way there are N elements.

now, suppose MATLAB allowed you to index your arrays just like we do in the DSP textbooks. nearly always we include 0 as an index (as in $x[0]$) and often we include "negative" time or negative frequency, for mathematical convenience. so in DSP we can have a discrete-time function $$x[n] \qquad -128 \le n < 128$$ and likewise a discrete-frequency function $$X[k] \qquad -128 \le k < 128$$ where the center of our windowed time function is at $n=0$ and DC is at $k=0$.

  1. normally, we define the Discrete Fourier Transform (which is what the FFT is a "fast" implementation of) and its inverse as

$$\begin{align} X[k] = & \sum\limits_{n=0}^{N-1} x[n] \, e^{-j 2 \pi \frac{nk}{N}} \qquad & 0 \le k < N \\ \\ x[n] = \frac{1}{N} & \sum\limits_{k=0}^{N-1} X[k] \, e^{+j 2 \pi \frac{nk}{N}} \qquad & 0 \le n < N \\ \end{align}$$

if you examine the expressions for both $X[k]$ and $x[n]$ you will see that both are periodic with period of $N$:

$$\begin{align} x[n+N] = & x[n] \qquad & -\infty < n < \infty \\ \\ X[k+N] = & X[k] \qquad & -\infty < k < \infty \\ \end{align}$$

this means, for any integers $n_0$ and $k_0$ the DFT and inverse can be expressed as:

$$\begin{align} X[k] = & \sum\limits_{n=n_0}^{n_0+N-1} x[n] \, e^{-j 2 \pi \frac{nk}{N}} \qquad & k_0 \le k < k_0+N \\ \\ x[n] = \frac{1}{N} & \sum\limits_{k=k_0}^{k_0+N-1} X[k] \, e^{+j 2 \pi \frac{nk}{N}} \qquad & n_0 \le n < n_0+N \\ \end{align}$$

Now, if Cleve and company had wished to make MATLAB truly compatible with the mathematical expressions already existing in the literature, and to make MATLAB flexible so that the user can define what is the minimum and maximum index for each array, they could have defined the MATLAB array variables with user-definable origins where the default value of the origin is 1 for each dimension of a newly created array, just to be backward compatible.

and being locked into this "1-origin" has forced MATLAB to redefine the DFT as:

$$\begin{align} X[k] = & \sum\limits_{n=1}^{N} x[n] \, e^{-j 2 \pi \frac{(n-1)(k-1)}{N}} \qquad & 1 \le k \le N \\ \\ x[n] = \frac{1}{N} & \sum\limits_{k=1}^{N} X[k] \, e^{+j 2 \pi \frac{(n-1)(k-1)}{N}} \qquad & 1 \le n \le N \\ \end{align}$$

and this should have been a big red flag to Cleve and TMW that there was something wrong in their fixed indexing scheme, because it forces users to remember that indices will be off by one when using the fft() and they should have known that signal processing engineers often benefit by modeling problems and signals and systems with negative time or negative frequency. like we sometimes want to display negative frequency along with positive frequency (just like we do in the textbooks and the lit) and, since the Fourier transform treats time and frequency identically (except for the sign of $j$, but that is purely a conventional difference), then we want to display negative time along with positive time. as such, this is an extremely useful and elegant manner to define and make use of the DFT:

$$\begin{align} X[k] = & \sum\limits_{n=-N/2}^{N/2-1} x[n] \, e^{-j 2 \pi \frac{nk}{N}} \qquad & -N/2 \le k < N/2 \\ \\ x[n] = \frac{1}{N} & \sum\limits_{k=-N/2}^{N/2-1} X[k] \, e^{+j 2 \pi \frac{nk}{N}} \qquad & -N/2 \le n < N/2 \\ \end{align}$$

  1. if MATLAB had adjustable index origins embedded into their array variables, then this would be painless. because MATLAB's fft() and ifft() (or i would call it "dft()" and "idft()") would be fully aware of the $n_0$ and $k_0$ inherent to the arrays supplied and can do exactly the right mathematical operation on those arrays. and if you were to use max() or min() or find() on the results of fft(), the index returned would be exactly the value you need for your mathematics. but, because MATLAB is the way it is, you must maintain a separate array of indices and use one level of indirection to get what you want.

so, suppose you had a simple chore of windowing a snippet of audio, FFTing it, finding the frequency of the peak amplitude, and displaying both time and the DTFT. your MATLAB code might look like:

x_input = load('x.mat');          % load time-domain input
N = 2*floor(length(x_input)/2);   % make sure N is even
x = x_input(1:N);

t = linspace(-N/2, N/2-1, N);              % values of time in units of samples
omega = linspace(-pi, pi*(1-2/N), N);      % values of (normalized) angular frequency

X = fftshift( fft( fftshift( x.*hamming(length(x)) ) ) );

[X_max k_max] = max( abs(X) );

figure(1);
plot(t, x, 'g');

figure(2);
plot(omega, abs(X), 'b');
hold on;
plot(omega(k_max), X_max, 'or');
hold off;

even though i haven't (yet) run this code, i think it will work. the fftshift() is necessary to put the first half of the time series x, into the "negative time" position (which is the last half going into the fft()) and to put the latter half of the output of the fft() into the "negative frequency" position in X.

now if MATLAB allowed adjustable index origin with a simple call to a function i'll call reorigin() (similar to reshape()), we wouldn't need the t or omega arrays nor the fftshift() operations and the indices of our x and X arrays would mean what the literature would expect them to mean.

(comp.dsp people know where i am on this issue, as it is about my principal pet peeve regarding MATLAB. and it's a shame that Octave and R appear to have followed along.)

$\endgroup$
  • 2
    $\begingroup$ How is this relevant ? Op asked "how to switch the light on" and you answer with "Why do Maxwell Equations are what they are". $\endgroup$ – Pier-Yves Lessard Apr 21 '17 at 21:29
  • $\begingroup$ the reason fftshift() exists and the operation behind it is necessary is really a technically and historically nuanced reason (which is why, outside of MATLAB, there is no intrinsic "FFT shift" or, better, "FFT swap halves" operation in the DSP textbooks). $\endgroup$ – robert bristow-johnson Apr 22 '17 at 5:05
1
$\begingroup$

fftshift only changes the order of the input "samples" or "coefficients". If applied on real data, the output is real.

It is meant to be applied on some fft result, not on the analysed time signal. Standard fft outputs do not emphasize on the natural symmetry around the $0$ index frequency, which fftshift corrects. Applying fftshift on the time signal does some wrapping which can be inconsistent if the signal does not exhibit FFT symmetries.

So a more natural use is:

plot(abs(fftshift(fft(t))));
$\endgroup$
0
$\begingroup$

There are some wonderful answers here. I think the difficulty here is you dont know what a fourier transform is. I also have reason to believe that this is one of the first times you are handling an audio signal in matlab. I'll try to keep my answer low key

In short the fourier transform is a transform to the frequency domain. The frequency domain contains information as to which frequencies are predominant in the signal, in speech and audio processing you call them modes (and even more names), in image and video processing you look at a whole frequency domain image (kind of looks like a shining bright star at the end of a tunnel). Point is how your signal changes in the spatial domain is reflected in the frequency domain.

Thing is, this frequency domain is itself repeating, hence when you take the fft, you see two copies (since fft here is the DFT, I'm also omitting a lot when I say that the freq domain is repeating, your DSP textbook has more). You dont want to see "two" copies unnecessarily do you? So you would want to circularly shift the fft and delete half the portion so that you can look at the nice long frequency denoting peaks.

So you do that shift by fftshift, its made to be used with the fft itself.

Does that help?

$\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.