how to find the right frequency axis when we take the Fourier transform of a function?

I want to know how to find the correct frequency axis when I am taking the fourier transfor. To make it clear, we know that for sinc function, we have the below relationship between the lenght of the sinc function and its fourier (the rect function):

based on this fact, each time I tried different frequency axix and then finally I find the correct one. But, actually I still do not know how to find the correct axis. This was a simple function that I already knew about the lenght of the function and its transforemed but for other functions I do not have any idea about the correct frequency axis.

Any help will be appreciated

Answer 1

[Paul]

Hi Donya,

By default, the frequency variable for fourier in the Symbolic Math Toolbox is in rad/sec.

Also, Matlab definition of sinc in the SMT is sinc(x) = sin(pi*x)/(pi*x).

So with the example in the question we have:

syms t omega f T real
s(t,T) = sinc(2*t/T);
S_omega(omega,T) = simplify(fourier(s(t,T),t,omega));

Now, if we want the CTFT expressed in Hz

S_f(f,T) = simplify(S_omega(2*sym(pi)*f,T))

S_f(f, T) = 

Now make some plots for T = 2

figure
fplot(s(t,2),[-3 3]),grid
xlabel('t (sec)')

figure
fplot(S_f(f,2),[-3 3]),grid
xlabel('f (Hz)')

As expected and consistent with the plot in the Question, the mainlobe of s(t) has a span of T the CTFT spans -1/T to 1/T.

See fourier and sympref for more discusson the how the SMT parameterizes the CTFT and its inverse, and how to control the parameterization.

[2NOR_Kh]

my question is more theory based. I think for having the correct frequency axis in plotting the fourier function, we should have nyquist frequency as a parameter.

After trying several times and changing the parameters, I finally could write a code to give me the lengths I was looking for, I'll share it in this comment.

It's been three months that I started studying signal processing, so I need to know what is the exact theory behind the relationship between the x axis in the fourier transformed function with its counterpart in the time domain.

close all;
clear all;
Fs=80;%same as the Nyq rate
Ts=1/Fs;
t=-1:Ts/2:1;
N=length(t);
f=5;
cx=sinc(t*f);
figure;
plot(t,cx);
xlabel('time');
ylabel(' magnitude');

%ff=fourier(cx(t,f));
fy=(fft(cx));
dx = (t(end)-t(1))/(N-1); % Time increment
Nyq = 1/(2*Ts); % Nyquist frequency
df = 1/(N*Ts);  % Frequency increment
if mod(N,2) == 0 % N is even
  k = -Nyq : df : Nyq-df;
else            % N is odd
  k = [sort(-1*(df:df:Nyq)) (0:df:Nyq)];
end
%fr=(-(N-1)/4:0.5:(N-1)/4); use this one or k
figure;
 plot(k,fftshift(abs(fy)));
 xlabel(' frequency t^{-1}');
 ylabel(' magnitude');

[Paul]

At least three different Fourier transforms are of interest: continuous time Fourier transform (CTFT), discrete-time Fourier transform (DTFT), and discrete Fourier transform (DFT). As the question had no mention of the Nyquist frequency, I assumed it was about the CTFT and answered as such. But now it's apparent the question is really about the DFT.

It looks like there is an error in the posted code

clear all;
Fs=80;%same as the Nyq rate
Ts=1/Fs;
t=-1:Ts/2:1;
N=length(t);
f=5;
cx=sinc(t*f);
figure;
plot(t,cx);
xlabel('time');
ylabel(' magnitude');

%ff=fourier(cx(t,f));
fy=(fft(cx));
dx = (t(end)-t(1))/(N-1); % Time increment

The Nyquist freuqency is 1/2 the sampling frequency. In the code above, the sampling frequency is 1/(Ts/2), so Nyq really should be

%Nyq = 1/(2*Ts); % Nyquist frequency
Nyq = 1/(2*Ts/2); % Nyquist frequency

Similiarly, df should be computed from the sampling frequency, which is 1/(Ts/2)

%df = 1/(N*Ts);  % Frequency increment
df = 1/(N*Ts/2);  % Frequency increment

if mod(N,2) == 0 % N is even
  k = -Nyq : df : Nyq-df;
else            % N is odd
  k = [sort(-1*(df:df:Nyq)) (0:df:Nyq)];
end
%fr=(-(N-1)/4:0.5:(N-1)/4); use this one or k
figure;
plot(k,fftshift(abs(fy)));
xlabel(' frequency t^{-1}');
ylabel(' magnitude');
xlim([-10 10])

I strongly recommend that the symbots Ts and Fs be used for the sampling periond and sampling frequency respectively, which I believe is standard notation and, IMO, makes the analysis much more clear. With this update, and some other modifications to the code we have

clear all;
Fs=160;  % sampling frequency
Ts=1/Fs; % sampling period
t=-1:Ts:1;
N=numel(t); % numel recommended instead of length
f=5;
cx=sinc(t*f);
figure;
plot(t,cx);
xlabel('time');
ylabel(' magnitude');

fy=(fft(cx));
Nyq = Fs/2;  % Nyquist frequency is 1/2 of the sampling frequency
% dx = (t(end)-t(1))/(N-1); % Time increment, should rename to dt, but not used
% suggest not using colon with df stride, stick colon with unit stride
if mod(N,2) == 0 % N is even
  k = ( (-N/2) : ((N-2)/2) )/N*Fs;
else            % N is odd
  k = ( (-(N-1)/2) : ((N-1)/2) )/N*Fs;
end
figure;
plot(k,fftshift(abs(fy)));
xlabel(' frequency t^{-1}');
ylabel(' magnitude');
xlim([-10 10])

[2NOR_Kh]

Hi Paul,

Thank you for your comment. I just didn't get why you're saying the sampling frequency is 1/(Ts/2).

The sampling frequency is Fs which is equal to 1/Ts. So the Nyq frequency should be half of this value, and in my code I wrote as:

Fs=80;
Ts=1/Fs;
Nyq = 1/(2*Ts) % Nyquist frequency

Nyq = 40

Also, the length of the rect function should be 5 but with this change, the length is 7 which is not correct.

[Paul]

The original code in the question is:

Fs=80;%same as the Nyq rate
Ts=1/Fs;
t=-1:Ts/2:1;

The time interval from one sample to the next is Ts/2 on the third line. Based on that, and the comment on the first line, I assumed that goal was to have the Nyquist frequency (not rate) be 80, which implies the sampling frequency is 160, which makes the sampling period 1/160 = 1/80/2.

Perhaps you really meant to have:

Fs=80;    % sampling frequency
Ts=1/Fs;  % sampling period
t=-1:Ts:1;
Nyq = 1/Ts/2

Nyq = 40

% or
Nyq = Fs/2

Nyq = 40

I'm not sure what what you mean by the length of the rect function being 7 and not 5, but I think you need to get Fs, Ts, and t sorted out first.

[2NOR_Kh]

if we assume the length of the sinc is T then the length of the rect must be 2/T.

[Paul]

Please go back and look at my original answer.

If we want the width, as defined by the first zero crossings, of the sinc function to be T, and the width of its Fourier transform to be be 2/T, then we must define the signal as

y(t) = sinc(2*t/T)

if using Matlab's definition of the sinc function.

In your case, we have

f=5;
cx=sinc(t*f);

Therefore f = 2/T, or T = 2/f = 2/5 = 0.4.

The width of the sinc function shown on the plot in this comment is 0.4 as exepcted, and the width of the Fourier transform in that same comment is 2/0.4 = 5, excactly as expected.

[2NOR_Kh]

there is a reason behind why I am insisting to solve it by using the signal processing parameters, such as sampling rate, or nyquist frequency, because I shoulld apply some filters on RF data in frequency domain and I have sampling rate, and time interval. So, this sinc function is just a simple example for me to know how to work with fourier on my data. But I don't know how to correct my code by your first solution.

[Paul]

I'm quite confused now. i already showed how to correct your code, at least the code posted above here. If you still have questions about this specific problem for a sinc function, post your current code and explain why you think it's not working. If you're talking about a different problem, perhaps you should post the code for that problem and explain why you think it's not working.

[Paul]

Maybe a different example will be clearer on how to calculate the frequency vector

Let the signal of interest be a sine wave with period P

P = 0.1;

Let the sampling period be

Ts = P/10;

Therefore, the sampling frequency (Hz) is

Fs = 1/Ts;

Define a vector of sample times, spaced by Ts

N = 50;
t = (0:(N-1))*Ts;

Compute the samples of the sine wave at the sample times

y = sin(2*pi/P*t);

Compute the FFT and apply fftshift()

Y = fft(y);
Y = fftshift(Y);

Compute the associated frequency vector in Hz for the fftshifted samples of Y.

Note that the frequency vector depends only on the number of samples and the sampling frequency, not on the signal itself, so it doesn't matter if we're working with a sine, or a sinc, etc.

if mod(N,2) == 0 % N is even
  k = ( (-N/2) : ((N-2)/2) )/N*Fs;
else            % N is odd
  k = ( (-(N-1)/2) : ((N-1)/2) )/N*Fs;
end

Now plot

figure
plot(k,abs(Y),'-o'),grid
xline(1/P,'r')
xline(-1/P,'r')
xlim([-15 15])
xlabel('Frequency (Hz)')

We see that the FFT spikes at the frequency of the sine wave (1/P), as expected.

[2NOR_Kh]

Sorry Paul, in my previous comment, I didn't read your comment correctly.

I wanted to check something else with you. I read that for normalizing the DFT coefficients, we should divide them to the length of the vector, but I don't know the scientific reason behind this. If you could help me to know why I'll appreciate it.

[Paul]

Before talking about normalization, I'd like to point out another issue to consider. The code we have now is correct for plotting the amplitude of the DFT. However, the phase plot would be incorrect. The reason for this is that the DFT (as implemented by fft()) assumes that the time domain samples start at n = 0, but in our case the samples of cx start at some value of n < 0. So we need to use the shift property to correct the output of fft() (we ccould also use ifftshfit(),which would actually be better for this problem).

Here is the code to illustrate

clear all;
Fs=160;  % sampling frequency
Ts=1/Fs; % sampling period
t=-1:Ts:1;
N=numel(t); % numel recommended instead of length
f=5;
cx=sinc(t*f);
% figure;
% plot(t,cx);
% xlabel('time');
% ylabel(' magnitude');

fy=(fft(cx));
Nyq = Fs/2;  % Nyquist frequency is 1/2 of the sampling frequency
if mod(N,2) == 0 % N is even
  k = ( (-N/2) : ((N-2)/2) )/N*Fs;
else            % N is odd
  k = ( (-(N-1)/2) : ((N-1)/2) )/N*Fs;
end

figure;
plot(k,abs(fftshift(fy)));
xlabel(' frequency t^{-1}');
ylabel(' magnitude');
xlim([-10 10])

figure
plot(k,imag(fftshift(fy)))

We see that the imaginary part of fy is nonzero, which should not be the case as we are expecting a rect function.

But if we correct fy by the appropriate phase shift, we get the expected result with imag(fy) = 0 (essentially)

N0 is the shift in cx that makes the input to fft() 0-based

N0 = -round(t(1)/Ts)

N0 = 160

fy = fy.*exp(1j*(0:(N-1))/N*2*pi*N0);
figure
plot(k,imag(fftshift(fy)))

I will come back later and respond to the question about normalization.

[2NOR_Kh]

However, in my work I totally ignore the phase, but it was a great point to have in mind in case I wanted to work on the phase one day.

Thank you again Paul, you helped me a lot.

[Paul]

for normalizing the DFT coefficients, we should divide them to the length of the vector,

Instead of "normalizing," I thnk a better word is "scaling." Whether or not we "should" scale the DFT coefficients depends on a) what we want the DFT coefficients to represent, and b) how the DFT is defined for the particular implementation being used.

One motivation for dividing the DFT coefficients by the length of the input vector, N, can be illustrated as follows. Consider again the case of the sine wave with period P and we take the DFT of time domain samples from 5 periods.

P = 0.1;
Ts = P/10;
Fs = 1/Ts;
N = 50;
t = (0:(N-1))*Ts;
y = sin(2*pi/P*t);
figure
stem(t,y)

Compute the FFT and apply fftshift()

Y = fft(y);
Y = fftshift(Y);

Compute the associated frequency vector in Hz for the fftshifted samples of Y.

if mod(N,2) == 0 % N is even
    k = ( (-N/2) : ((N-2)/2) )/N*Fs;
else            % N is odd
    k = ( (-(N-1)/2) : ((N-1)/2) )/N*Fs;
end

Now plot

figure
stem(k,abs(Y)),grid
xlim([-15 15])
xlabel('Frequency (Hz)')

We see that, using Matlab's fft() implementation, the DFT amplitudes at +-10 Hz are 25. Suppose now we want to take those DFT coefficients and reconstruct the original 5 periods of the sine wave. We know that the DFT coefficients multiply complex exponentials of the form exp(1j*2*pi*k*n*Ts). Note that n*Ts = t and, in this case, we only need to worry about k = -10 and k = 10. After the fftshift, Y(21) corresponds to k = -10 and Y(31) to k = 10

yrecon = Y(21)*exp(1j*2*pi*(-10)*t) + Y(31)*exp(1j*2*pi*10*t);
figure
stem(t,yrecon)

This plot is the same as that above, except the amplitude is too high by a factor of 50, which not coincidentally is the value of N

N

N = 50

So, if we want to reconstruct the sine wave, we can divide Y by N

Y = Y/N;
figure
stem(k,abs(Y))

yrecon = Y(21)*exp(1j*2*pi*(-10)*t) + Y(31)*exp(1j*2*pi*10*t);
figure
stem(t,yrecon)

After scaling by N, the DFT coeffcients have amplitide of 1/2 that of the sine wave.

An alternative approach to reconstructing the sine wave is to divide by N in the reconstruction equation instead of scaling the DFT coeffeicients

Y = Y*N; % get back the original DFT
yrecon = (1/N) * ( Y(21)*exp(1j*2*pi*(-10)*t) + Y(31)*exp(1j*2*pi*10*t) );

which we won't plot because we know what it will look like.

Now look the documenation for fft, particularly the equation for the ifft. We see that equation includes the 1/N factor in front of the summation of the complex exponentials, exactly as shown in the yrecon equation immediately above.

In summary, using Matlab's implementation of fft/ifft, you can scale the DFT coefficients output from fft by 1/N if you want them to be "consistent" in some sense with the amplitude of the exponentials that would be used to reconstruct the signal. I think some people call this "amplitude scaling." IMO, I don't think you should ever explicitly scale Y by N, i.e., do not do

Y = fft(y)/N

because then you run the risk of forgetting later in the analysis thay Y is scaled and something else won't work, e.g., ifft(Y) yields the wrong answer if you forget to multiply the result by N. Instead, just divide by N on the fly as needed, e.g.,

plot(k.abs(Y)/N)

if you decide you want to see the "amplitude scaled" coefficients.

The CTFT, DTFT, and DFT are all related to each other, so if you change the definition of the DFT by scaling it by 1/N, you may need to reconsider how to maintain its relationshihp to the CTFT and DTFT as well if want all three to be consistent with each other.

Finally, before you do any scaling you really need to know how the DFT is defined for the particular software you're using. As we've seen, Matlab scales the DFT sum by 1 and the iDFT sum by 1/N (I think this convention is very common). However, the only requirement is that the product be 1/N, so other s/w might have the 1/N scale factor on the DFT sum and unity scale factor on the iDFT sum (or 1/sqrt(N) on both!). In this case, obviously, you wouldn't want to divide Y by N.