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

2 次查看（过去 30 天）
Mehdi
Mehdi 2024-4-18
编辑： Star Strider 2024-4-26
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)|')

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

采纳的回答

Star Strider
Star Strider 2024-4-18
编辑：Star Strider 2024-4-18
Ran in:
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.
.
  2 个评论
Mehdi
Mehdi 2024-4-26
编辑：Mehdi 2024-4-26
Ran in:
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
Star Strider 2024-4-26
编辑：Star Strider 2024-4-26
Ran in:
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.
.

请先登录，再进行评论。

更多回答（1 个）

AH
AH 2024-4-18
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

类别

在 Help Center 和 File Exchange 中查找有关 Correlation and Convolution 的更多信息

标签

产品


版本

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by