Using fft in matlab to retrieve input signal magnitude

Question

I am trying to understand how you use an FFT in matlab to measure/retrieve the magnitude of the input signal.

For my understanding I have put this script together, creating a 15Hz sine wave with a sampling frequency of 1kHz.

fftsize    = 2048;            % How many samples for FFT
fs         = 1000;            % Sampling frequency (Hz)
sim_length = fftsize/fs ;     % run sim for this time
t          = 0:1/fs:sim_length-1/fs;                
x          = (0.5)*sin(2*pi*15*t); % 15 Hz component, 0.5 magnitude

Window and FFT:

%% Windowed FFT
data                  = x'; % Assign input data
window                = blackman(fftsize);
windowed_data         = window.*data;

% Perform FFT and take magnitude (of complex result)
fftd_data             = abs( fft(windowed_data, fftsize) );

The input frequency was not coherent, add in the windowing and there is going to be some spectral leakage, therefore I need to sum the power over a few bins around where I expect the fundamental to be.

% Each bin represents
bin_freq = fs/fftsize ;

% 15 Hz fundamental is centred in bin 
fundamental_bin = ceil( 15/bin_freq );

% Allow spectral spreading of +-4 bins
lower_limit = fundamental_bin - 4;
upper_limit = fundamental_bin + 4; 

Sum of values from fft:

sum( fftd_data(lower_limit:upper_limit) )

ans =

  512.3738

Which is not the 0.5 I was looking for. What is the best way to get the (0.5) input magnitude from the fft result.

After some digging about I have come up with the following but not quite sure why it works or if I have named the calculations appropriately.

incoherent_power_gain = sum(window.^2);
power_spectrum        = fftd_data.^2 ;
power_fftd_data       = (4 * power_spectrum) / (incoherent_power_gain * fftsize);

sum_fundamental_power = sum( power_fftd_data(lower_limit:upper_limit) )

sum_fundamental_power =

    0.2500

We converted to Power so we can take the square root to get back to a signal level.

sqrt( sum_fundamental_power )

ans =

    0.5000

Great, but is there a better way or can some one explain:

If power spectrum is the correct name for fftd_data.^2.

I am also unsure what is happening in (4 * power_spectrum) / (incoherent_power_gain * fftsize).

Answer 1

I just fixed 3 things in your code and added comments where I did. The reason is basicly just the way the FFT is defined/implemented in Matlab.

fs         = 1024;            
fftsize    = 2048;           % Choose as a power of 2 of the sample rate       

sim_length = fftsize/fs ;     
t          = 0:1/fs:sim_length-1/fs;                
x          = 0.5*sin(2*pi*15*t);

data                  = x';
window                = blackman(fftsize);
windowed_data         = window.*data;


fftd_data   = abs( fft(windowed_data, fftsize)/length(x) ); % Normalize the fft

bin_freq = fs/(fftsize) ;

fundamental_bin = ceil( 15/bin_freq );

lower_limit = fundamental_bin-4;
upper_limit = fundamental_bin+4; 

2*sum( fftd_data(lower_limit:upper_limit) ) % Multiply by 2, because of the DFT-symmetry

ans =

    0.5000

For a more detailed explanation, check out this thread at Mathworks forum.