Could anyone explain the example provided in thr help for understanding the use of fft?

I am trying to understand the using of fft in matlab, regarding the provided example in help I could not figure out what the f is and how it is defined as f = Fs*(0:(L/2))/L. Is there anybody who could explain why the f is defined as f = Fs*(0:(L/2))/L?

Fs = 1000;            % Sampling frequency                    
T = 1/Fs;             % Sampling period       
L = 1500;             % Length of signal
t = (0:L-1)*T;        % Time vector
S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
X = S + 2*randn(size(t));
plot(1000*t(1:50),X(1:50))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('t (milliseconds)')
ylabel('X(t)')
Y = fft(X);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;
plot(f,P1) 
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('|P1(f)|')

Answer 1

[Star Strider]

This calculation of the frequency vector for the displayed Fourier transform is not the easiest to understand:

f = Fs*(0:(L/2))/L;

The frequency vector generally extends from 0 Hz (the units are actually cycles/(independent variable unit)) to ½ the sampling frequency, called the Nyquist frequency. In this instance, the fft l;ength is the same as the signal length ‘L’ so that is used to calculate the frequency vector. It extends from 0 to ‘L/2’ appropriately.

Going further:

syms L Fs 
f = Fs*[0 1 2 3 (L/2)]/L

f = 

so the ‘f’ vector goes from 0 to the Nyquist frequency.in steps of ‘Fs/L’ (the default increment for the colon, : operator being 1).

The Nyquist frequency is the highest frequency at which a sampled signal can uniquely be described.

This then gets into sampling theory and details that are best left to texttbooks on digital signal processing. That quickly gets complicated, so I will stop here and encourage you to pursue that on your own.

Fs = 1000;            % Sampling frequency                    
T = 1/Fs;             % Sampling period       
L = 1500;             % Length of signal
t = (0:L-1)*T;        % Time vector
S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
X = S + 2*randn(size(t));
plot(1000*t(1:50),X(1:50))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('t (milliseconds)')
ylabel('X(t)')

Y = fft(X);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L

f = 1x751

0 0.6667 1.3333 2.0000 2.6667 3.3333 4.0000 4.6667 5.3333 6.0000 6.6667 7.3333 8.0000 8.6667 9.3333 10.0000 10.6667 11.3333 12.0000 12.6667 13.3333 14.0000 14.6667 15.3333 16.0000 16.6667 17.3333 18.0000 18.6667 19.3333

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

plot(f,P1) 
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('|P1(f)|')

EDIT — (18 Apr 2024 at 20:19)

Clarification added.

.

[Mehdi]

I am trying to analyse the fft of my response to get the natural frequencies as below:

20.9768932686948

31.6086353836362

32.3365811045944

38.4556687582899

38.8388607498704

45.1217285865391

But fail.

What do you think the problem cn be?

load('WXYN1')
S = 28451;
dt=  1.8171e-06;
plot((01:S)*dt,WXYN1(01:S))

Fs = 1/dt;            % Sampling frequency                    
L=size(WXYN1(1:S-2),2);
% L=S-1;
% plot([dt:dt:L*dt],WXYN1(1:S-2))
Yy = fft(WXYN1(1:S-2));
P2 = abs(Yy/L);
P1 = P2(1:L/2+1);

Warning: Integer operands are required for colon operator when used as index.

P1(2:end-1) = 2*P1(2:end-1);
f = (Fs)*(0:(L/2))/L;
plot(f,P1) 
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('|P1(f)|')

[Star Strider]

Calculating it differently, I do not believe your method is failing, there simply is not nuch information in the signal.

My code —

load('WXYN1')
whos  -file WXYN1

  Name       Size                Bytes  Class     Attributes

  WXYN1      1x28451            227608  double              

dt = 1.8171e-06;
Fs = 1/dt;
L = numel(WXYN1);
t = linspace(0, L-1, L)/Fs;

figure
plot(t, WXYN1)
grid
xlabel('Time (s)')
ylabel('Amplitude')
title('WXYN1')

Fn = Fs/2;
NFFT = 2^(nextpow2(L)+6)                                    % For Efficiency & Increased Frequency Resolution

NFFT = 2097152

wf = hann(L);                                               % Window Function
FTs = fft((WXYN1(:)-mean(WXYN1)).*wf, NFFT)/sum(wf);        % Fourier Transform (Subtract 'mean' So Other Peaks Are More Easily Visible)
% FTs = fft(WXYN1(:).*wf, NFFT)/sum(wf);                      % Fourier Transform
Fv = Fs*(0:NFFT/2)/NFFT;
Iv = 1:numel(Fv);

figure
plot(Fv, abs(FTs(Iv))*2)
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Fourier Transform of ‘WXYN1’')

[pks,locs] = findpeaks(abs(FTs(Iv))*2, 'MinPeakProminence',2.5E-3);
Major_Peaks = table(Fv(locs).', pks, 'VariableNames',{'Frequency (Hz)','Magnitude'})

Major_Peaks = 3x2 table
    Frequency (Hz)    Magnitude
    ______________    _________

        45.398        0.0089002
        66.654          0.00977
        489.93        0.0033595

DNF = [20.9768932686948; 31.6086353836362; 32.3365811045944; 38.4556687582899; 38.8388607498704; 45.1217285865391];
DNF_Mag = interp1(Fv, abs(FTs(Iv))*2, DNF);


figure
plot(Fv, abs(FTs(Iv))*2, 'DisplayName','Fourier Transform')
hold on
stem(DNF, DNF_Mag, 'v-r', 'filled',  'DisplayName','Desired Natural Frequencies')
hold off
grid
legend('Location','best')
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Fourier Transform of ‘WXYN1’ (Annotated)')
xlim([0 5E+2])
text(Fv(locs).', pks, compose('\\leftarrow Freq = %.2f  Hz,   Mag = %.6f',[Fv(locs).', pks]), 'Rotation',90)

figure
semilogx(Fv, mag2db(abs(FTs(Iv))*2))
grid
xlabel('Frequency (Hz)')
ylabel('Power (dB)')
title('Fourier Transform of ‘WXYN1’')

WXYN1_det = detrend(WXYN1,2);

figure
plot(t, WXYN1_det)
grid
xlabel('Time (s)')
ylabel('Amplitude')
title('WXYN1 Detrended')

FTs = fft((WXYN1_det(:)-mean(WXYN1_det)).*wf, NFFT)/sum(wf);        % Fourier Transform (Subtract 'mean' So Other Peaks Are More Easily Visible)
% FTs = fft(WXYN1(:).*wf, NFFT)/sum(wf);                      % Fourier Transform
Fv = Fs*(0:NFFT/2)/NFFT;
Iv = 1:numel(Fv);

figure
plot(Fv, abs(FTs(Iv))*2)
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Fourier Transform of ‘WXYN1’')

[pks,locs] = findpeaks(abs(FTs(Iv))*2, 'MinPeakProminence',2.5E-3);
Major_Peaks = table(Fv(locs).', pks, 'VariableNames',{'Frequency (Hz)','Magnitude'})

Major_Peaks = 3x2 table
    Frequency (Hz)    Magnitude
    ______________    _________

        21.518           0.1729
        60.093         0.020701
        489.93        0.0033691

DNF_Mag = interp1(Fv, abs(FTs(Iv))*2, DNF);

figure
plot(Fv, abs(FTs(Iv))*2, 'DisplayName','Fourier Transform')
hold on
stem(DNF, DNF_Mag, 'v-r', 'filled',  'DisplayName','Desired Natural Frequencies')
hold off
grid
legend('Location','best')
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Fourier Transform of ‘WXYN1’ (Annotated)')
xlim([0 5E+2])
ylim([0 0.8])
text(Fv(locs).', pks, compose('\\leftarrow Freq = %.2f  Hz,   Mag = %.6f',[Fv(locs).', pks]), 'Rotation',90)

There is a considerable D-C offset that cannot completely be elimiinated by subtracting the mean of the signal. Some peaks are still evident, although the desired natural frequency peaks are not evident, even with significantly increased frequency resolution.

I thought that detrending it might help (and it does, to an extent), however only one of your desired resonant frequencies is now evident. If you do this experiment again, it might help to increase the sampling frequency, and to provide the actual sampling times as a vector, however that may not be a practical possibility. (I gave some thought to filtering the signal to detrend it, since that is cometimes preferable as a detrending method, however there is not enough frequency information to design an appropriate filter.)

EDIT — (26 Apr 2023 at 10::53)

Minor update to last figure.

.

Answer 2

[AH]

When you take a L-point DFT of a signal sampled at rate , then the frequency interval is uniformly sampled. The distance between two frequency points is . Hnece, the discretized frequency is given by . For a real-valued signal, it's sufficient to look at the positive side of the spectrum (a.k.a one-sided spectrum) as its Fourier transform symmetric with respect to frequency 0. Hence, the one sided frequency vector is .

Hope this clarifies the main question.

Two caveats below are worth mentioning:

1. How to deal with the right edge.
2. In the code above, the spectrum of teh lowpass-equivalent (analytic signal) is shown.

For further detailed discussion on thistopic you may want to take a look at this example: Practical Introduction to Frequency-Domain Analysis