Is this a mistake in documentation or my misunderstanding of FFT?

조회 수: 1 (최근 30일)
Fahad
Fahad 2023년 2월 15일
댓글: Fahad 2023년 2월 16일
In MATLAB Help Center for FFT here; there is one example whose code is given below.
My question is why at line 6 from below (in the following code) it is P1(2:end-1) = 2*P1(2:end-1); instead of P1(2:end) = 2*P1(2:end)?
My rational is that while calculating the Single sided Amplitude spectrum from Double sided spectrum it is only the DC component to be excluded from multiplying by 2. (All AC component is supposed to be multiplied by 2.) If so then why did the 5001 element in vector P1 (which is an AC component) is keep off from getting multiplied by 2 ?
regards
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;
figure
plot(f,P1)
title("Single-Sided Amplitude Spectrum of X(t)")
xlabel("f (Hz)")
ylabel("|P1(f)|")
Y = fft(S);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1); % Why is it NOT P1(2:end) = 2*P1(2:end);
figure
plot(f,P1)
title("Single-Sided Amplitude Spectrum of S(t)")
xlabel("f (Hz)")
ylabel("|P1(f)|")

채택된 답변

Bruno Luong
Bruno Luong 2023년 2월 15일
편집: Bruno Luong 2023년 2월 15일
"Single sided Amplitude spectrum from Double sided spectrum it is only the DC component to be excluded from multiplying by 2"
No if L is even, the nyquist frequency is Y(L/2+1) or P1(end) is shared between one-sided spectrum and the "negative" frequency so it must not devided by 2. So TMW code is correct.
On the otherhand if L is odd the Niquist frequency is in between 2 bins and skip. The correct code is
P1 = P2(1:(L+1)/2)
P1(2:end) = 2*P1(2:end);
Of course one can do the same code for both, a little bit more complicated
P1 = P2(1:ceil((L+1)/2));
P1(2:end-mod(L+1,2)) = 2*P1(2:end-mod(L+1,2));
  댓글 수: 2
Paul
Paul 2023년 2월 16일
@Bruno Luong is correct, of course, wrt to amplitude at the freuqency samples. But there's another subtlety to consider.
If you don't multiply the first and last points by two, then you have to be careful about how to interpret what's happening between the first two points and between the last two points if you're going to "connect the dots" to estimate the Discrete Time Fourier Transform. That's because every point in between the first two and between the last two should be multiplied by two as well, but if you connect the dots between the first and second points, the connecting line won't give correct interpolation. Same for the last two points.
Fahad
Fahad 2023년 2월 16일
Oh, so spectrum is already adding up from positive and negative frequencies at P1(end) (or Nyquist frequency: fs/2).
Cool ! How come I was never aware of this subtle and key concept in all my 20 years with DSP.

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Get Started with MATLAB에 대해 자세히 알아보기

태그

제품


릴리스

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by