781 views (last 30 days)

Show older comments

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
on 15 Apr 2020

dpb
on 12 Sep 2017

Edited: dpb
on 24 Jan 2018

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.

dpb
on 13 Sep 2017

Le Dung
on 24 Jan 2018

Edited: Walter Roberson
on 24 Jan 2018

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
on 24 Jan 2018

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.

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

Start Hunting!