MATLAB: Scaling the FFT and the IFFT

fftifftMATLABscalescaling

What is the correct way to scale results when taking the Fast Fourier Transform (FFT) and/or the Inverse Fast Fourier Transform (IFFT)?

Best Answer

Does scaling matter?

In most situations, scaling is really not all that important. The overall shape of the spectrum matters much more than the absolute scale.

What are the conventions?

But if you really are worried about it, there are several different conventions from which you can choose (see definitions below):

Scale by dt for the FFT, and by Fs for the IFFT
Scale by 1/M for the FFT, and by M for the IFFT
Scale by 1 for the FFT, and by 1 for the IFFT
Scale by 1/sqrt(M) for the FFT, and by sqrt(M) for the IFFT.
and so on.

I generally use either option #1 or option #2 depending on my mood and whether it's raining outside.

Definitions

Here I am assuming that I have a discrete-time signal x represented as an M x N matrix, where M is the number of samples and N is the number of channels.

   [M,N] = size(x);

Furthermore, I am assuming that the sampling rate is Fs and that I have defined the time increment as

   dt = 1/Fs;

and the frequency increment as:

   dF = Fs/M;

What do all these conventions have in common?

All of these conventions have one thing in common: The product of the two scaling factors is always 1. Please note that the ifft function in MATLAB includes a scaling factor of 1/M as part of the computation, so that the overall round-trip scaling is 1/M (as it should be).

HTH.

Rick

Related Solutions

MATLAB: Physically correct Normalization of fft + Implementing Parseval’s Theorem

Hello fgb

Here is a revised answer that attempts to better address the issue. I believe there are a couple of things going on here. In your code,

   df = 2*Fs/L

should not have the factor of 2. Also, if the length of the x array is not a power of 2, then the procedure

NFFT = nextpow2(L)
Y ~ fft(x,NFFT)

changes the length of the x array that gets fft'd. So both uses of L afterwards have to be changed, and df has to become

df = Fs/NFFT.

For a transient response the transform and the energies in the time and frequency domains are

Y = fft(x)*dt                % dt is the same as T
ex = sum(conj(x).*x)*dt      % dt is the same as T 
eY = sum(conj(Y).*Y)*df

The first two you have. The frequency domain calculation eY agrees with ex, and it uses the entire frequency array. All frequency components have identical energy content. Each energy is the abs value squared of the amplitude. Life is good. The issue is with using positive frequencies only.

When you fft a real quantity and double the positive positive frequency components Y(f), you are going to a representation (just at one frequency for simplicity)

x = 2 |Y(f)| cos(2pi f t + phase angle)

But the energy content is different, because the average value of cos(x)^2 is 1/2. So if somebody has handed you a positive-frequency-absolute-value spectrum like the one that you have after this step

YA(1)=YA(1)/2;  %Correction for the fact that DC doesn't occur twice

then you also have to make sure that point N/2+1 has been included (as you did), and that that point is not doubled. So there needs to have been a step

YA(end)=YA(end)/2;  %Correction because nyquist point doesn't occur twice

Then the expression

eY_new = (YA(1)^2 + (1/2)*sum(YA(2:end-1).^2) + YA(end)^2)*df;

works. But it's not straightforward.

Stuff like this is why I try to stay away from positive-frequencies-only as much as possible for real calculations. It's not bad for illustrative plotting purposes.

MATLAB: I get the peaks at the wrong frequency and with the wrong amplitude.

hello Anna

look at my fft function at the end of my code

clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% data
data = importdata('SFC5S_nov25_ST_1_1.txt');
dt = 1e-6; % 1 micro seconds
signal = data(:,1);
samples = length(signal);
Fs = 1/dt; % sampling frequency (Hz)
%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 1;
if decim>1
    signal = decimate(signal,decim);
    Fs = Fs/decim;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 512;    % 
OVERLAP = 0.95;
% spectrogram dB scale
spectrogram_dB_scale = 80;  % dB range scale (means , the lowest displayed level is XX dB below the max level)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% options 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if you are dealing with acoustics, you may wish to have A weighted
% spectrums 
% option_w = 0 : linear spectrum (no weighting dB (L) )
% option_w = 1 : A weighted spectrum (dB (A) )
option_w = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% display 1 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[freq, sensor_spectrum] = myfft_peak(signal,Fs,NFFT,OVERLAP);
% convert to dB scale (ref = 1)
sensor_spectrum_dB = 20*log10(sensor_spectrum);
% apply A weigthing if needed

if option_w == 1
    pondA_dB = pondA_function(freq);
    sensor_spectrum_dB = sensor_spectrum_dB+pondA_dB;
    my_ylabel = ('Amplitude (dB (A))');
else
    my_ylabel = ('Amplitude (dB (L))');
end
figure(1),plot(freq,sensor_spectrum_dB,'b');grid
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
xlabel('Frequency (Hz)');ylabel(my_ylabel);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% display 2 : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[sg,fsg,tsg] = specgram(signal,NFFT,Fs,hanning(NFFT),floor(NFFT*OVERLAP));  
% FFT normalisation and conversion amplitude from linear to dB (peak)
sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg));     % NB : X=fft(x.*hanning(N))*4/N; % hanning only
% apply A weigthing if needed
if option_w == 1
    pondA_dB = pondA_function(fsg);
    sg_dBpeak = sg_dBpeak+(pondA_dB*ones(1,size(sg_dBpeak,2)));
    my_title = ('Spectrogram (dB (A))');
else
    my_title = ('Spectrogram (dB (L))');
end
% saturation of the dB range : 
% saturation_dB = 60;  % dB range scale (means , the lowest displayed level is XX dB below the max level)
min_disp_dB = round(max(max(sg_dBpeak))) - spectrogram_dB_scale;
sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;
% plots spectrogram
figure(2);
imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
axis('xy');colorbar('vert');grid
title([my_title ' / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(fsg(2)-fsg(1)) ' Hz ']);
xlabel('Time (s)');ylabel('Frequency (Hz)');
function pondA_dB = pondA_function(f)
	% dB (A) weighting curve
	n = ((12200^2*f.^4)./((f.^2+20.6^2).*(f.^2+12200^2).*sqrt(f.^2+107.7^2).*sqrt(f.^2+737.9^2)));
	r = ((12200^2*1000.^4)./((1000.^2+20.6^2).*(1000.^2+12200^2).*sqrt(1000.^2+107.7^2).*sqrt(1000.^2+737.9^2))) * ones(size(f));
	pondA = n./r;
	pondA_dB = 20*log10(pondA(:));
end
function  [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
% FFT peak spectrum of signal  (example sinus amplitude 1   = 0 dB after fft).
% Linear averaging
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer overlap % (between 0 and 0.95)
signal = signal(:);
samples = length(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,1);
    s_tmp((1:samples)) = signal;
    signal = s_tmp;
    samples = nfft;
end
% window : hanning
window = hanning(nfft);
window = window(:);
%    compute fft with overlap 
 offset = fix((1-Overlap)*nfft);
 spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows
    % main loop
    fft_spectrum = 0;
    for i=1:spectnum
        start = (i-1)*offset;
        sw = signal((1+start):(start+nfft)).*window;
        fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft);     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end
    fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling
% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_spectrum = fft_spectrum(select);
freq_vector = (select - 1)*Fs/nfft;
end

Best Answer

Related Solutions

MATLAB: Physically correct Normalization of fft + Implementing Parseval’s Theorem

MATLAB: I get the peaks at the wrong frequency and with the wrong amplitude.

Related Question