can someone explaing the computation for double sided and single sided spectrum in fft() example
86 次查看(过去 30 天)
显示 更早的评论
Fs = 1000;
T = 1/Fs;
L = 1500;
t = (0:L-1)*T;
S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
X = S + 2*randn(size(t));
Y = fft(X);
P2 = abs(Y/L); %why are they dividing by L
P1 = P2(1:L/2+1); %Not sure what is being done here
P1(2:end-1) = 2*P1(2:end-1); %Not sure what is being done here
f = Fs*(0:(L/2))/L; % not sure what is going on here as well
%I want to have the frequency span from - to + so that I can see both delta functions of cosine and sine
plot(f,P1)
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('|P1(f)|')
%I want to plot cos(2*pi*60*t) in frequency domain with two delta functions at -60 freq and +60 freq
2 个评论
Shah Mahdi Hasan
2020-4-15
I actually did write a Medium article which essentially covers all of your questions regarding this type of example. Have a look: https://medium.com/@tanweer.mahdi/breaking-down-confusions-over-fast-fourier-transform-fft-1561a029b1ab
回答(2 个)
dpb
2017-9-12
编辑:dpb
2018-1-24
The example is single-sided from F=0:Nyquist (Fs/2). Since the FFT returned by Matlab is double-sided, from -Fmax:+Fmax,
f = Fs/2*linspace(0,NFFT/2+1);
as the example uses is the positive frequency vector. NFFT/2+1 is the location of the Nyquist frequency component in the returned FFT; (+)ive frequencies start from the origin of the array df apart positive with higher index to Nyquist, then negative with -Fs at the next element past NFFT/2+1 --> NFFT/2+2, the midpoint plus one being -(Fs/2-df).
The 1/L normalization is the normalization of the output for the length of the input signal; in the example there are L actual data points with power associated with them but NFFT points in the FFT--the normalization to return the amplitude of the PSD to match the input is only for the number of elements in the vector with data; not including the additional NFFT-L points that are just augmented zeros.
To do the two-sided plot, just change the range for f and the subscripting for the FFT data--
f=Fs/2*linspace(-1,1,NFFT);
plot(f,2*abs([Y(NFFT/2+2:end) Y(1:NFFT/2+1)]))
I get the following figure for this case--
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/178235/image.jpeg)
Again, read the code for the example at FFT very carefully to see the difference between L and NFFT...the former is the length of the actual time series (1000 in the example, 1500 in your case) whereas the latter ( NFFT ) is the length of the transform which was determined using nextpow2 to get the next higher power-of-two length series for the (slight) efficiency in computation that makes internally. Thus NFFT would be 1024 in the example but 2048 for your case as 1500>1024 so must augment to next.
But, for proper normalization, only the length of the original signal L is used as all the additional points are simply introduced zeros and so to divide by NFFT where NFFT > L will reduce the apparent power by that ratio.
If you use the default argument for size in the call to FFT, then the routine will return the same length transform as the input time series and you should set
NFFT=L;
to get proper normalization.
ADDENDUM
Oh, the factor of two in the plot is since Matlab FFT returns double-sided, half the input energy is in each the positive and negative frequencies...strictly speaking the factor there should be included only for the one-sided plot but I kept it just for presentation purposes to make the actual peak match the input time series amplitude -- in reality, it should be unity for the two-sided so the sum of the two peaks is the total.
3 个评论
dpb
2017-9-13
My bad...the FFT return is two-sided and the midpoint is at NFFT/2+1 but it's the Fs point in the middle, not the DC component. I knew what I meant but didn't write it correctly, nor plot it right. As you can see, from 1:NFFT gives the mirror image from left to right. See updated answer for correction.
Le Dung
2018-1-24
编辑:Walter Roberson
2018-1-24
Hi dpb.
I have a problem about DC value.
At code above,
Fs = 1000;
T = 1/Fs;
L = 1500;
t = (0:L-1)*T;
S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
X = S + 2*randn(size(t));
Y = fft(X);
P2 = abs(Y/L);
i create two axis-freq:
f1 = Fs/N*(0:N-1); % f1=0:N1
f2 = Fs/N*(-N/2:N/2-1); % f2=-Fs/2:F2/2 or (-Fmax : + Fmax as you said above)
and i plot:
subplot(2,1,1);
plot(f1,P2)
subplot(2,1,2);
plot(f2,P2)
We will receive two the same double-sided spectrum. (you an see below) So, in first case, f=0:N-1, so does DC value lie in f=N/2 ? and -Fmax lie in f=0, +Fmax lie in f=N-1?
and second case, does DC value lie in f=0 ? and -Fmax lie in f=-N/2, +Fmax lie in f=N/2?
thank you so much!||
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/181307/image.png)
1 个评论
dpb
2018-1-24
编辑:dpb
2018-1-24
Be better to ask question as comment rather than an Answer but I'm too harried to move it at the moment...
The two-sided FFT is as discussed above from DC-->Fmax-->DC with negative frequency in the first location in the returned vector.
To plot on absolute frequency scale -Fmax:+Fmax with DC in the middle of the axis, you have to split and flip the PSD data in the middle and plot against a continuous frequency vector. This is illustrated in my previous answer altho I see there was a typo in that the plot command had a variable "ff" instead of just "f" as it should have been; guess my ol' fat fingers double-hit the keyboard in copying/pasting somehow.
If you'll run that example as given, you'll generate the same plot that shows the two peaks where on expects on the frequency axis.
Note carefully the discussion about the use of L and NFFT if one uses a different length for the transform as that of the input signal and the question of normalization of the peak for double-sided vis a vis one-side spectra as to whether or not to include the 2X factor for total energy.
Also note for the two-sided example, the example you want to use from above is the last where the linspace() term is over Fx*[-1:1] range, not the initial discussion that is one-sided and shows how the data are arranged.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Spectral Measurements 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!