FFT example on MATLAB help

Hi everybody, I am trying to learn FFT in MATLAB by understanding the example available in help file(<http://www.mathworks.co.uk/help/techdoc/ref/fft.html>). I don't know why in the line : Y = fft(y,NFFT)/L; the fft result is divided by L.

your help is apreciated

Answer 1

[Rick Rosson]

You don't need to divide it by L, it is purely a matter of scaling the result by a constant, which does not affect the shape of the spectrum, but really only affects the units of measure. The choice of scaling is largely a matter of convention rather than anything of significance.

Some people will disagree with this assessment and argue that the scaling does matter. The reality is that it depends to a large degree on your objectives and what you need to get from taking the FFT.

[Alzapoa]

Thanks Rick for the nice response. Actually the magnitude does matter in my case. I need to see how different are my magnitudes from time domain to frequency domain. when I remove divide by L, and for simplisity let me take the noise out of the game, the amplitude of the harmonics are 700 and 1000 for 50Hz and 120Hz respectively, but I know my time domain amplitude were 0.7 and 1. so to me seems all the existing harmonics are added L times at the output of the FFT. consequently, if i need the ampliyudes like the real wold the divide by L is inevitable. is n't it?

[Juan]

importand of time sampled

I guessed that the longer one samples, the most exact resault with fft. And so it seems to do with frequency, but not at all with amplitude. The amplitud decrees as long as longer the time sampled.

Did I do something wrong ? Is this reasonable?

Thanks in advance

% T=[.01 .02 .1 .2 1];
% T=[1:10]
T=[1 5 10];
n=10; % maximal number of harmonic

freq = 5;
Fs = 150e4; % Sampling frequency (frecuencia de muestreo)
nfft = 1e7; % Length of FFT % number of fft bins

tic

for m=1:length(T) 

      t = 0:1/Fs:T(m); % Time vector of 1 second
      x = cos(2*pi*t*freq); % sine wave of f Hz.
      x = x-mean(x);% restamos de la fcn 'x' su DC offset

      figure(1);
      subplot(length(T),1,m);
      plot(t,x,'.');

      % Take fft, padding with zeros so that length(X) is equal to nfft
      X = fft(x,nfft);
      % FFT is symmetric, throw away second half
      X = X(1:nfft/2);
      % Scaling is done here using the number of samples: length(x)/2
      X=X/(length(x)/2);
      % Take the magnitude of fft of x
      mx = abs(X);
      % Frequency vector
      f_fft = (0:nfft/2-1)*Fs/nfft;

      f_Axis=f_fft(1:nfft/2);

      % picos de amplitud
      pks= findpeaks(mx(1:nfft/2));
      n=min(n,length(pks));
      pks_sort=sort(pks,'descend');% vector de picos

      for k =  1:n
          locs=find(mx(1:nfft/2)==pks_sort(k));
          f(k)=f_Axis(locs);
      end

      % bar plot freq y magnitude
      figure
      uds='Amp';
      pks_sort=pks_sort(:); % to make it column array for text of bar plot
      f=f(:); % to make it column array for text of bar plot
      subplot(1,2,1)
      bar(pks_sort(1:n),'r');
      ylabel({sprintf('Amplitude (%s)'...
          ,uds)},'fontweight','bold','fontsize',16);
      text(1:n,pks_sort(1:n),num2str(pks_sort(1:n),'%.2f'),... 
          'HorizontalAlignment','center',...
          'VerticalAlignment','bottom')

      subplot(1,2,2)
      bar(f(1:n),'b');
      ylabel('Frequency (Hz)','fontweight','bold','fontsize',16);
      text(1:n,f(1:n),num2str(f(1:n),'%.2f'),...
          'HorizontalAlignment','center',...
          'VerticalAlignment','bottom')

      bar3d_pks_sort(m,:)=pks_sort(1:n);
      bar3d_f(m,:)=f(:);
  end

figure;
col(1,:)=[0 0 1]; col(2,:)=[0 .5 0];
plot(bar3d_pks_sort(:,1),'-+','Color',col(1,:),'Linewidth',2);
h1 = gca;
h2 = axes('Position',get(h1,'Position'));
plot(bar3d_f(:,1),'-o','Color',col(2,:),'Linewidth',2);
set(h1,'YColor',col(1,:))
set(h2,'YAxisLocation','right','Color','none','XTickLabel',[],'YColor',col(2,:))
set(h2,'YAxisLocation','right','Color','none','XTickLabel',[])

Answer 2

[E K]

when you write the command Y=fft(y,NFFT) you calculate the fft of y on NFFT and when you divide it by L you just divide the FFT matrix.

lets say a=fft(y,NFFT) what you are doing basicly

a/L.

[Alzapoa]

Ok! but the question is what is the meaning of this dividation? I mean why we need to divide the fft matrix by L?

Answer 3

[Honglei Chen]

This is basically done to preserve the power at each frequency sample point. The original series has L samples in it. At each frequency sample point, L copies of signal at corresponding frequency are coherently added together via FFT. So to preserve the power, you need to divide by L.

This is best seen when there is no noise involved

Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 1000;                     % Length of signal
t = (0:L-1)*T;                % Time vector
x = 0.7*sin(2*pi*Fs/8*t) + sin(2*pi*Fs/4*t); 
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(x,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);

% Plot single-sided amplitude spectrum.
plot(f,2*abs(Y(1:NFFT/2+1))) 
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')

[Alzapoa]

Thak you Honglei for the nice response. so let me understand your statement "At each frequency sample point, L copies of signal at corresponding frequency are coherently added together via FFT" let consider dc frequency (k=0Hz) and recall fft definition from MATLAB help: http://www.mathworks.co.uk/help/techdoc/ref/fft.html. the exponential function stays as 1 and the value of the fft function is summation of N sample of the function x[j] in time domain. right?

[Honglei Chen]

Yes and No. There are N samples added together. But because your L is less than N, the signal is zero-padded. Therefore, in terms of power, you only have L effective samples. That's why you need to divide by L, not N to preserve the power.

Since you bring up the DC point, I have to mention that the way the DC is treated is not entirely correct in the example. To get the one-sided spectrum, you don't need to scale both DC and Nyquist frequency as these two points are unique.

Answer 4

[Wayne King]

Both Honglei and Rick have given you good responses. If you want the magnitudes recovered from the DFT to equal the time domain amplitudes: yes, you have to scale by the length of the input vector and multiply by 2 if you have a real-valued signal, because the real-valued signal results in complex exponentials scaled by 1/2.

   Fs = 1000;
   t = 0:1/Fs:1-1/Fs;
   x = 0.7*cos(2*pi*50*t)+ cos(2*pi*100*t);
   xdft = fft(x);
   % the DFT bin for 50 Hz is 51
   % the DFT bin for 100 Hz is 101
   amp50 = 2/length(x)*xdft(51);
   amp100 = 2/length(x)*xdft(101);
   abs(amp50)
   abs(amp100)

[Alzapoa]

Thanks Wayne, your complementary note is actually what I needed to be sure about the magnitude of the fft in f-domain. so what if I have a range of discrete data instead on function x (x=[a,b,c,d,...,z]). I mean do I still need to divide the fft output by length of x to get the amplitudes the same in the time domain?

Answer 5

[ajay munaga]

Compute the 8 point DFT of the sequence x(n)={ 1, 0, 1, 0, 0.5, 0, 0.5, 0} using radix-2 DIF FFT algorithm. Implement using MATLAB