can someone explaing the computation for double sided and single sided spectrum in fft() example

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

[Shah Mahdi Hasan]

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

[RACHAEL AMOAH]

Wonderful! Thank you, very helpful website.

Answer 1

[dpb]

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

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.

[Syed Rumman]

Thanks for responding. I am still confused about few things. I still don't understand the plotting in frequency domain. My questions are:

f=Fs/2*linspace(-1,1,NFFT); 
plot(f,2*abs(Y))% why are you multiplying by 2?

- When I try running the previous code, Matlab can't find NFFT variable? Is this a predefined variable or am I suppose to declare it?

- Can you clarify the 1/L normalization again please? or direct to some resource? (didn't follow the role of NFFT there.)

[Syed Rumman]

So, I plotted with frequency range that you suggested.

f = Fs/2*linspace(-1,1,length(Y));
plot(f,2*abs(Y));

and I got this. But I was expecting the delta functions to be at -50and+50 and -120and+120. But they are at the very end of the spectrum. Can you explain what I did wrong here? Also for the y-axis, I don't think it shows right power value .7 and 1. Is there way to show those value? Thanks

[dpb]

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.

Answer 2

[Le Dung]

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!||

[dpb]

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.