What is this fft code doing?

3 views (last 30 days)
kcee N
kcee N on 6 Jan 2022
Edited: David Goodmanson on 10 Jan 2022
Hi I am new to Matlab and performing signal processing. I am trying to understand what this code is doing? how and why are we are determining the indexNyquist, why are we scaling, trancating and compensating for trancating? Thank you
spectrum = fft(Signal,k); %generate spetrum of signal with FFT to k points
indexNyquist = round(k/2+1); %vicinity of nyquist frequency
spectrum= spectrum(1:indexNyquist); %truncate spectrum to Nyquist frequency
spectrum = spectrum/(length(Signal)); %scale spectrum by number of points
spectrum(2:end) = 2 * spectrum3(2:end); %compensate for truncating negative frequencies, but not DC component

Answers (2)

dpb
dpb on 6 Jan 2022
See the documentation of FFT and read the example for a one-sided PSD from the FFT carefully...
Basically, the reason for the above is because the FFT returns a two-sided DFT from -Fmax : +Fmax.
So to return only the one-sided positive frequency, take only the first half of the output. That happens to correspond to Fmax which is the Nyquist frequency by definition so the writer chose to use that as a variable name.
Since the base FFT output is two-sided, all the energy in the input signal is spread over the full range of 2*Fmax (-Fmax -> +Fmax), so if we only keep the positivie half of the output, the magnitude of the remainder must be doubled to conserve the input signal total energy, hence the factor of two.
The nuance there is that there is only one DC component in the output FFT so it is not doubled -- hence the starting index of 2 in the (2:end) indexing expression. There is, however, also only one Fmax component and so that expression should be (2:end-1) in general. It depends on whether k is odd or even as to whether the round() expression picks up an extra element or not; generally I will truncate instead of round.
  3 Comments
dpb
dpb on 9 Jan 2022
OK, I was loose on power/energy, I'll agree. The above was referenced to the FFT example from which the code snippet appears to have been extracted -- it illustrates amplitude normalization.
For illustration, let's revert back to the example there --
>> 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);
>> S=S+pi/10; % let's introduce a DC component for to see where it goes...
>> P2=abs(fft(S))/L;
>> P1=P2(1:L/2+1);
>> P1(2:end-1)=2*P1(2:end-1);
>> P1(1:2) % and, as expected, the DC component is full magnitude
ans =
0.31 0.00
>> L=L-1 % let's go to odd length now
L =
1499.00
>> t = (0:L-1)*T; % Time vector
>> S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
S=S+pi/10;
>> P2=abs(fft(S))/L;
>> P1=P2(1:L/2+1); % now the odd number generates warning...we'll just ignore for now
>> P1(2:end-1)=2*P1(2:end-1);
>> P1(1:2) % but the DC component is still full amplitude
ans =
0.31 0.00
>> mean(S) % to confirm (I had "format bank" on, is why two decimal places)
ans =
0.31
>> find(P2>0.1) % locate the nonzero peaks other tha rounding noise...
ans =
1.00 76.00 181.00 1320.00 1425.00
>> P2(ans)
ans =
0.31 0.35 0.49 0.49 0.35
>>
There's still only the one DC component...

Sign in to comment.


David Goodmanson
David Goodmanson on 9 Jan 2022
Edited: David Goodmanson on 10 Jan 2022
Hi kcee,
The code (I'll refer to it by line below) is not doing what it is doing very well. Let's start with the good part, the fourth line, normalization. Take a cosine signal with a certain amplitude A.
A cos(2pif0t) = A ( e^(2pif0t) + e^(-2pif0t) ) / 2
which in the frequency domain has peaks of height A/2 at -f0 and +f0. If the signal has length N, then to obtain the peaks A/2 in the frequency domain**, you have to divide the fft output by N. That's due to Mathwork's choice for initial normalization of the fft, nothing either good or bad about it. Note that N has to be the length of the original signal, not some later truncated or padded signal.
**If there are an exact integral number of oscillations in the time domain window, you get A/2 exactly, and if not you get something close to A/2
Now for the not so good parts.
ooo For a signal with sampling frequency fs, the output frequency spacing of an fft is fs/N. That’s a hard and fast rule. N can be either even or odd, and taking 8 or 7 points as an example, the fft output corresponds to the following frequency array:
[0 1 2 3 ny -3 -2 -1]*(fs/8) even
[0 1 2 3 -3 -2 -1]*(fs/7) odd
If the fftshift function is used to put 0 frequency in the center of the array, then
[ny -3 -2 -1 0 1 2 3]*(fs/8) even
[-3 -2 -1 0 1 2 3]*(fs/7) odd
The nyquist frequency is the largest frequency that can be expressed with a given sampling frequency, that is fs/2. In the arrays above, ‘ny’ was shown as a label in place of 4. For even N, If the 4 is plugged in, that frequency is 4/8 = fs/2, nyquist. But for odd N there is no way to get fs/2, so there is no nyquist.
To be clear, a time array with an odd number of points clearly has a nyquist frequency. It’s just not available from the fft. All this means that for odd N, the second and third lines are incorrrect.
ooo If you are transforming a real signal, then the amplitudes of the negative frequencies are the complex conjugates of the amplitudes of the positive frequencies. So if you know one you know the other. What people will do is toss out the negative frequencies and double the amplitude at the positive frequencies, so A/2 --> A and you have the time domain amplitudes. This is all right for plotting purposes but I think it’s done way more often than necessary.
At f=0 you don’t have a positive and negative frequency, and since A cos(0) = A there is no doubling involved. It’s also true that for even N, there is no doubling needed for the nyquist frequency either. So you have to be careful to not-double the point corresponding to f=0 in all cases, and when N is even not- double the nyquist frequency either. But the fifth line doesn’t make that distinction and is incorrect.

Community Treasure Hunt

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

Start Hunting!

Translated by