1
$\begingroup$

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

$\endgroup$
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Thanks. In my final application I do not have control over the sample rate, that is determined by the system (8kHz) I am measuring which introduces a small error in to the final magnitude. $\endgroup$ – Morgan Oct 23 '13 at 7:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.