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

Question

2NOR_Kh 2022-6-22

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1746065-how-to-find-the-right-frequency-axis-when-we-take-the-fourier-transform-of-a-function

编辑： Paul 2022-6-26

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Paul 2022-6-22

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1746065-how-to-find-the-right-frequency-axis-when-we-take-the-fourier-transform-of-a-function#answer_991480

编辑：Paul 2022-6-23

在 MATLAB Online 中打开

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.

13 个评论
显示 11更早的评论隐藏 11更早的评论

2NOR_Kh 2022-6-24

在 MATLAB Online 中打开

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 2022-6-24

在 MATLAB Online 中打开

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])

Paul 2022-6-25

编辑：Paul 2022-6-25

在 MATLAB Online 中打开

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.

Paul 2022-6-25

在 MATLAB Online 中打开

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 2022-6-25

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 2022-6-26

编辑：Paul 2022-6-26

在 MATLAB Online 中打开

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.

请先登录，再进行评论。

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

13 个评论
显示 11更早的评论隐藏 11更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

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

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

13 个评论 显示 11更早的评论隐藏 11更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

13 个评论
显示 11更早的评论隐藏 11更早的评论