Correct choice of one sided frequency axis after fft

There are a lot of queries on fft frequency. I guess the following point not discussed anywhere explicitly (at least at one place). Hope someone can provide an insight here.

If we have and even number of data points, N=10, the fft output arranges the data as

fft = [c0, c1, c2, c3, c4, c-5, c-4, c-3, c-2, c-1], where the numbers are subscripts corresponding to positive and negative frequencies. I read somewhere that MATLAB calculates the negative coefficient first, hence we have c-5 but not c5. The author did not explain the reason.

Point no. 1, that the c values are not symmetric.

When we wish to make two-sided frequency spectrum, the frequency axis ranges from [-(N/2): (N/2)-1]*Fs/ N. Fs is the sampling rate, N is the number of even data points.

If we wish to make a one-sided positive frequency spectrum, should we choose

A) [0:(N/2)]*Fs/N and ignore the fact the we are using the values corresponding to the negative frequency axis, given that the data is a real number and it is just a mirror image.

B) [0: (N/2)-1]*Fs/N represents the true positive frequency axis?

If Fs= 250 Hz, the true positive frequency axis will end at 124.9980 Hz

If we happen to choose the negative frequency axis values and ignore the frequency sign, the frequency axis ends at 125 Hz exactly.

The same data when plotted in Origin ends the frequency axis at 125 Hz when plotted single sidedly.

Which approach is rigorously correct? Thanks.

[dpb]

Why would you choose anything but the actual DC to max positive frequency that is actually calculated? That there is one sample short of Fs/2 is just a fact of "the way things are".

The time history of N points beginning at T=0 is one element short of N*dt as well...maybe that realization will help to understand the "why" of it all...

[FW]

Why would you choose anything but the actual DC to max positive frequency that is actually calculated? That there is one sample short of Fs/2 is just a fact of "the way things are".

Yes it is. The key question is a conceptual problem rather. I am not worrying about the positive frequency axis ending at 124.9980 Hz rather than 125 Hz. but why does OriginPro have positive the frequency axis ending at 125 Hz exactly and MATLAB doesn't?

Are they including fft = [c0, c-1, c-2, c-3, c-4, c-5]?

[dpb]

I have no idea; you'll have to ask them or see what the values are that they have used are empirically and determine.

This is the wrong forum for Origin-specific queries; this is Matlab.

[FW]

Is it officially documented somewhere by Matlab that the fft out is in this order as shown below? This problem is very much relevant to Matlab's way of indexing, Origin was just used for comparison.

fft = [c0, c1, c2, c3, c4, c-5, c-4, c-3, c-2, c-1]

Thanks.

[Bruno Luong]

The fft formula is documented in black and white (and MATLAB is no different than many other SW).:

Discrete Fourier Transform of Vector

Y = fft(X) and X = ifft(Y) implement the Fourier transform and inverse Fourier transform, respectively. For X and Y of length n, these transforms are defined as follows:

Y(k)=SUM j=1...n X(j) W_n(j−1)​(k−1)

X(j)=1/n SUM k=1...n Y(k) W{n}-(j−1)​(k−1),

where

W_n=e(−2πi)/n

is one of n roots of unity.

The problem often arises when people try to interpret the formula as frequency and they often mess up or confuses themselve.

[FW]

Bruno, I am talking about the indexing and ordering of the fft output by Matlab. I have seen this Matlab page you linked. Also see a similar discussion here, and the problem is perhaps not as trivial.

https://dsp.stackexchange.com/questions/59316/one-sided-frequency-spectrum-matlab-vs-origin?noredirect=1#comment118369_59316

[Bruno Luong]

The index location / ordering is fully defined in this formula. Don't you see it ?

What is not trivial for you is when you talking about c5 and c-5 as if they are two different things, they actually represents the same function if you look at for formula and it located at y(6) (and it's Nyquis frequency). Period.

[dpb]

As the doc says, ML uses the fftw implementation so there's nothing unique about it other than ML does store all arrays with one-based indexing. That's simply bookkeeping.

ML returns the DC component in the first array position, then the fix(N/2) positive followed by the same number of negative frequencies in obverse order. If N odd, then the returned spectrum is symmetric about DC.

When N even, it is simply not possible to have a DC component and N/2 positive and negative frequency components...just the way it is. fftw will have the same "problem" and Origin will too unless they do something like either augment or truncate the input vector to be odd internally or simply add another element or somesuch on the output.

[Bruno Luong]

The DFT formula never tells which Y elements corresponds to "positive" or "negative" frequency. Often user interprets them as such but DFT never state signed frequency, since one can add an (k'*n), where k' is an arbitrary integer in "frequency" to give the same result.

In the example of n=10, the c-5, c5 are the same coefficient of the fast oscillated sequence (1,-1,1,-1, ....-1) and I can arbitrary view it as c25 or c1005, nothing is wrong with it since it gives the same result.

[dpb]

True, but if one is looking at it from the viewpoint of baseband frequency analysis as is the OP, the ML storage order begins with the DC component up, then the obverse.

>> x=randn(1,11);
>> y=fft(x.')
y =
   1.9397 + 0.0000i
   2.7291 + 0.7234i
  -0.6829 + 0.2768i
   4.2973 + 1.0150i
   3.0071 + 2.0347i
  -0.0981 + 2.2543i
  -0.0981 - 2.2543i
   3.0071 - 2.0347i
   4.2973 - 1.0150i
  -0.6829 - 0.2768i
   2.7291 - 0.7234i
>> mean(x)
ans =
    0.1763
>> y/11
ans =
   0.1763 + 0.0000i
   0.2481 + 0.0658i
  -0.0621 + 0.0252i
  ...

The first element returned by FFT after normalization is the mean of the input--the DC component if doing typical signal processing of a measured input.

[FW]

Hello dbp, I did a quick check with OriginPro double sided spectrum. The frequency scale is automatic there when sampling rate is provided. At 250 Hz, with OriginPro, the negative freq. ends at -124.9980 Hz, and positive frequency ends at +125 Hz. However, if we follow the Matlab output order, the positive frequency ends at +124.9980 Hz and negative frequency end is -125 Hz. It seems like a matter of convention and both approaches by Matlab or OriginPro are "mathematically" correct.

This was my Matlab code

t=[0:1/250:10.004]'; % ensuring even number output

Fs= 250; % Hz units

N= length(t);

g=normpdf(t, 5, 0.2); % Gaussian peak of unit area

f_Hz=[(-N/2):(N/2)-1]'*Fs/N;

G_double=fft(g);

fftshift_G = fftshift(G_double); % fftshift rearranges the vectors so that zero frequency is in the center

Y=abs(fftshift_G);

plot(f_Hz, Y, '.')

[Bruno Luong]

"if we follow the Matlab output order, the positive frequency ends at +124.9980 Hz and negative frequency end is -125 Hz."

Again no : MATLAB FFT doc/help doesn't never ever tell such thing. This is an interpretation by you (or whoever you read and believe). We can select either convention of frequency intervals

[-124.9980,125] Hz, or

[-125,124.998] Hz.

Both are equally correct, or even an odd choice of frequency interval:

[-125,124.998] + 250 Hz.

Odd but still correct.

[FW]

Again no : MATLAB FFT doc/help doesn't never ever tell such thing. This is an interpretation by you (or whoever you read and believe). We can select either convention of frequency intervals

[-124.9980,125] Hz, or

[-125,124.998] Hz.

This makes much more sense. My source of information was a primer on FFT for Physicists. The author does not cite anything. Basically you are saying that there is nothing in MATLAB which supports form of the table. This is the table I was following (full pdf attached), see page 10.

[Bruno Luong]

The paper cites about the choice of FFTSHIFT applied on even-data spectrum.

This choice is arbitrary. User can take it or leave it, in the later case placing the nyquise as positive if he/she prefers (that makes compliance with OriginPro).

[dpb]

I think your author goes out of his way to confound what is really quite simple... :)

I'd basically ignore the table as being of no benefit for any purpose under the sun.

As noted earlier, ML returns DC component in first array location; whether there are even or odd number in the series determines whether the other components are or are not symmetric in pairs or whether there's an "odd man out" when even.

As Bruno says, depending on your purpose you can either fftshift or not depending upon objective at hand. I rarely bother, simply take the first half as I don't recall in 40 years a need for anything but the one-sided spectrum. Others' needs clearly may differ... :)

[FW]

Thanks dpb, I guess it all boils down to convention. This FFT primer was pretty good and I thought this is the way Matlab arranges the output (as shown in the table). However I could not find any other page which confirmed it. It is now clear that this the author's own interpretation, which is not wrong either. I am not a mathematican nor a signal processing engineer. As a chemist, I was trying to implement a fft step in resolution enhancement protocols for a spectroscopic signals. The take home message is that fft output, say N=6, it is

fft = [c0, c1, c2, c3, c4, c5], where the number k refers to N-1 values.

for which we can construct the frequency axis as

-Fs/2: (Fs/2-Fs/N)

or -(Fs/2-Fs/N): Fs/2 is completely up to us, and both frequency axes are correct.

[dpb]

>> n=0:9;
>> [n;fftshift(n)]
ans =
     0     1     2     3     4     5     6     7     8     9
     5     6     7     8     9     0     1     2     3     4
>> n=n+1;
>> [n;fftshift(n);fftshift(n)-1]
ans =
     1     2     3     4     5     6     7     8     9    10
     6     7     8     9    10     1     2     3     4     5
     5     6     7     8     9     0     1     2     3     4
>> n=1:9;
>> [n;fftshift(n);fftshift(n)-1]
ans =
     1     2     3     4     5     6     7     8     9
     6     7     8     9     1     2     3     4     5
     5     6     7     8     0     1     2     3     4
>> 

[Bruno Luong]

@dpb: Don't understand what you are at by doing

fftshift(n)-1

???

Don't you want to show something along this line:

ffts_l=@(x) circshift(x,floor((size(x)-1)/2)) % OriginPro compatible
ffts_r=@(x) circshift(x,floor((size(x))/2));  % equilalent to MATLAB fftshift

>> x=0:8;
>> [x; ffts_l(x); ffts_r(x)]

ans =

     0     1     2     3     4     5     6     7     8
     5     6     7     8     0     1     2     3     4
     5     6     7     8     0     1     2     3     4

>> x=0:9;
>> [x; ffts_l(x); ffts_r(x)]

ans =

     0     1     2     3     4     5     6     7     8     9
     6     7     8     9     0     1     2     3     4     5
     5     6     7     8     9     0     1     2     3     4

[dpb]

Bruno, the n-1 is just the offset difference between 1- and 0-based arrays--a trivial difference but seems to throw many into complete tizzy.

Your latter illustration is also pertinent, I was just too lazy to write the extra code figuring simply showing OP the indexing is the root answer of how Matlab returns the order and one can build the quoted author's table from it by inspection.

[Bruno Luong]

"difference between 1- and 0-based arrays"

It escapes me about who was discussing of 1/0 base indexing in this thread.

As I understood the discussion is around the FFT output coefficients c_i and its subindex i which makes the correspondance to the frequency:

f(i) := Fs * i / n. 

and whereas i should vary beween

[-floor((n-1)/2) : ceil((n-1)/2)] % OriginPro

or

[-ceil((n-1)/2) : floor((n-1)/2)] % MATLAB FFTSHIFT

(As we discuss, both are mathematically correct and it's a matter of interpretation or arbitrary choice)

[Bjorn Gustavsson]

@dpb (on 7 Jul 2019 at 14:50): One "common" use of "not choosing" base-band is that one has undersampled a high-frequency signal with signal at frequencies inside the N-th Nyqvist zone. The spectra of that signal can then be accurately determined with sampling-rates much lower than the centre-frequency of the band.