1,438 views (last 30 days)

Show older comments

Jan
on 21 Jan 2019

Rick Rosson
on 13 Sep 2011

Raunaq Nayar
on 20 Feb 2020

What about the case where I perform scaling by factor(x) in time domain perform FFT followed by IFFT and then scale back by factor (1/x) to recover signal. Do I need to take care of the value of x then? I am asking this because in this case we will be scaling up and scaling down both in time domain, so I think for this specific case the Parseval's Theorum should hold true irrespective of the value of x.

In other words, I am following the following sequence of steps-

Scaling by x --> FFT --> IFFT -->Scaling by 1/x

I would also like to confirm the same for scaling by x in frequency domain followed by IFFT then followed by FFT and then scaling by 1/x in frequency domain itself.

In other words, I am following the following sequence of steps-

Scaling by x --> IFFT --> FFT -->Scaling by 1/x

Can x assume any value to recover the signal for the mentioned cases?

Dr. Seis
on 28 Jan 2012

You are right about scaling being unimportant if only the shape of the spectrum is desired. However, if it is necessary that the amplitudes in the frequency spectrum be correct, then there is only one option for scaling - your option #1. In order for Parseval's theorem to hold, the energy in the time domain must equal the energy in the frequency domain. The example below demonstrates this:

> N = 8;

> dt = 0.1;

> df = 1/(dt*N)

df =

1.25

> a = randn(N,1)

a =

0.70154

-2.0518

-0.35385

-0.82359

-1.5771

0.50797

0.28198

0.03348

> b = fft(a)*dt

b =

-0.32813

0.10746 + 0.30519i

-0.080365 + 0.075374i

0.34826 + 0.17802i

0.13866

0.34826 - 0.17802i

-0.080365 - 0.075374i

0.10746 - 0.30519i

> energy_a = sum(a.*conj(a) * dt) % Not necessary to use conj here

energy_a =

0.83314

> energy_b = sum(b.*conj(b) * df) % Necessary to use conj here

energy_b =

0.83314

Raunaq Nayar
on 20 Feb 2020

What about the case where I perform scaling by factor(x) in time domain perform FFT followed by IFFT and then scale back by factor (1/x) to recover signal. Do I need to take care of the value of x then? I am asking this because in this case we will be scaling up and scaling down both in time domain, so I think for this specific case the Parseval's Theorum should hold true irrespective of the value of x.

In other words, I am following the following sequence of steps-

Scaling by x --> FFT --> IFFT -->Scaling by 1/x

I would also like to confirm the same for scaling by x in frequency domain followed by IFFT then followed by FFT and then scaling by 1/x in frequency domain itself.

In other words, I am following the following sequence of steps-

Scaling by x --> IFFT --> FFT -->Scaling by 1/x

Can x assume any value to recover the signal for the mentioned cases?

Jérome
on 1 Aug 2017

Hello

I have seen so many subjects related to this issue and so don't know on which to put my answer...I hope this can help some people. I can see there are mainly 3 groups of people, those saying that is not important, those dividing by the number of points of the signal, and those dividing by the sampling frequency. first, for some applications,yes the correct amplitude is important. for example, the magnitude of an earthquake is computed from the amplitude of the spectrum.

then I agree with Dr Seis, the correct way of scaling spectrum is multiplying by dt.

People saying fft has to be divided by the number of points often take the example of the sin wave with amplitude A and want to see 2 peaks with amplitude A/2 on the spectrum. however, this is not the Fourier transform of a continuous sin wave. The Fourier transform has two Diracs. The value of each peak is infinite and the integration over the frequency domain is A/2. So the value of a correctly-scaled discrete spectrum we should have on both peaks is

A/df/2 = Npoints*A/Fs/2

Fs being the sampling frequency, df the step of the frequency vector.

the matlab fft outputs 2 pics of amplitude A*Npoints/2 and so the correct way of scaling the spectrum is multiplying the fft by dt = 1/Fs. Dividing by Npoints highlights A but is not the correct factor to approximate the spectrum of the continuous signal.

The second point is the parseval equation. I have seen many people saying the fft can not respect this relation or is not applicable in discrete mode. first, in discrete mode, if should tends to the continuous value. if can not be Npoints lager or smaller. And, if the fft is multiplied by dt, the energy of the input signal equals the energy of the spectrum.

I have seen quite often people using the Parseval equation for discrete signal like this, which is incorrect

sum(abs(xi).^2) = sum(abs(Xi).^2) with X = fft(x)

The correct discrete form of the Parseval relation is:

sum(abs(xi).^2)*dt = sum(abs(Xi).^2)*df

the relation is satisfied if the fft is multiplied dt and df is correctly defined.

Moreover, there are many simple typical Fourier transforms such as exponential decay, triangle function.. you can model the temporal signals and the known continuous transforms and check that fft*dt is the correct way of approaching the continuous transform.

Alex Hanes
on 14 May 2020

In MATLAB, your statement is basically correct (e.g. the correct discrete form of Parseval's theorem is:

P_avg_TD = sum(abs(x).^2)*dt;

P_avg_FD = sum(abs(X).^2)*df;

provided that the Fourier transform result is appropriately scaled by a factor of dt (your samling period). You can see this by:

where and so the above equation becomes:

Canceling terms on the left and right, you get the "incorrect" expression you wrote above:

where is just defined as the fft result (). For reference, Numerical Recipes uses the "incorrect" definition as the discrete form of Parseval's theorem, so maybe there is confusion for varying definitions and conventions.

Here's an example code that shows the correct Fourier transform scaling with comparison to the analytical result for the Fourier transform of a Gaussian. The analytic solution to the integral:

clear; clc;

%% Analytic Function: Gaussian

sigma = 1; % "width" of Gaussian

% Define the time-domain function, |f(t)|^2:

fun_t = @(t) exp(-t.^2/(2*sigma^2)).*conj(exp(-t.^2/(2*sigma^2)));

% Define the analytic Fourier transform of f(t), e.g. |f(w)|^2

fun_w = @(w) (sigma*sqrt(2*pi)).*exp(-(sigma*w).^2/2).*conj(...

(sigma*sqrt(2*pi)).*exp(-(sigma*w).^2/2));

int_t = integral(fun_t,-inf,inf); % P_avg of f(t)

int_w = integral(fun_w,-inf,inf)/(2*pi); % P_avg of f(w)/2pi

%% Numerical Example: Gaussian

% Define time domain for numerical signal

N = 1e7 ; % no. points in domain

domain = [-5000 5000]; % domain to evaluate f(t)

L = domain(1,2) - domain(1,1); % length of domain (e.g. [-5e3,+5e3])

dt = L/N; % step size, dt

t = min(domain):dt:max(domain)-dt; % t_n = t_min + n*dt (n = 0, ..., N-1)

% Define the function in the time-domain

sigma = 1; % "width" of Gaussian

ft_num = exp(-t.^2/(2*sigma^2)); % f(t), numerical

fw_num = fft(ft_num); % f(w), numerical (unscaled)

% Consider Parseval's Theorem (Numerical Recipes, no dt or df factors)

E1_timedomain = sum(abs(ft_num.^2))

E1_freqdomain = sum(abs(fw_num.^2))/N

% Consider Parseval's Theorem (based on "Riemann" integration)

df = 1/(N*dt); % calculate frequency step

fx = fw_num.*dt; % scale fw_num according to dt (redefine as fx)

E2_timedomain = dt*sum(abs(ft_num.^2))

E2_freqdomain = df*sum(abs(fx.^2)) % factor of 2*pi is dropped by

% differences in FT definitions*

I should also point out that the "incorrect" method satisfies Parseval's theorem (e.g. P_avg_TD = P_avg_FD in the first line of code), however, you get the incorrect average power content, as shown by the example code for E1_timedomain versus E2_timedomain.

Chani
on 16 Sep 2015

Edited: Chani
on 16 Sep 2015

Hi everyone,

I tried both options mentioned above, namely option #1 and option #2 with a simple sine curve. The results of both options confuse me and I am hoping you can help clear things up.

The code I am referring to is the following:

dt = 0.05; %Time step

t = [0:dt:10]; %Time vector

s = sin(2*pi*t); %Signal vector

N = length(s); %Number of data points

f_s = 1/dt; %Sampling frequency

df = f_s/N; %Frequency increment

NFFT = N; %FFT length used as second argument of fft()

y_option1 = fft(s,NFFT); %Compute FFT using sampling interval as scaling factor

y_option1 = y_option1*dt;

y_shiftOption1 = fftshift(y_option1);

y_option2 = fft(s); %Compute FFT using signal length as scaling factor

y_option2 = y_option2/N;

y_shiftOption2 = fftshift(y_option2);

if mod(N,2) == 0 %N is even

f = -f_s/2 : df : f_s/2-df;

else %N is odd

f = [sort(-1*(df:df:f_s/2)) (0:df:f_s/2)];

end

Plotting the sine signal and both result vectors y_shiftOption1 and y_shiftOption2 leads to the following Figure:

From my point of view the result of option #2 exhibits the correct amplitude of 0.5V, since the time signal has an amplitude of 1V and this value is split into the two impulses on the negative and the positive frequency scale.

However, regarding Parseval's Theorem, the energy is only preserved using option #1:

sum(abs(s).^2)*dt = 5

sum(abs(y_shiftOption1).^2).*df = 5

sum(abs(y_shiftOption2).^2).*df = 0.0495

Does this mean, if I am interested in the real amplitude of the signal, I have to use option #2 and if I want to preserve the energy of the signal option #1 has to applied?

Thanks in advance for your help!

Best Regards Chani

Kenny
on 14 Feb 2018

Doc Seis.....have you got an example to show this? I'm thinking of a portion of a sinusoidal signal that digitally spans 1 period..... such as N points, where point #0 is the start point, and the point #(N-1) is the LAST point (such that the following data point would theoretically be the same as point #0 again). The data set will be 'N' points.

For this case, the scaling should still be 'divide by N', right? If instead we multiplied the result by the sampling period, then the FFT magnitude amplitudes wouldn't look right..... right?

Michele Marconi
on 24 Oct 2017

How is it possible that the fft computed with matlab (no particular scaling!) and the FFT obtained by utilizing Xilinx module for SystemGenerator, have a 2000x factor scaling difference? i.e. Xilinx one is 2000 times lower amplitude?

(I thought this topic was spot-on regarding amplitude issues, sorry for gravedigging!)

Hossam Fadeel
on 23 May 2019

I have a situation where the values of the spectrum from the MATLAB FFT is double the values of the Xilinx FPGA FFT for the same input data.

I used below MATLAB Code,

clc;

clear all;

close all;

NFFT=2048; %NFFT-point DFT

fileID = fopen('Sine2048_test.txt','r');

A_Data = fscanf(fileID, '%d');% 2048 Sample points

%% Generate Sine with different freucies frpom the samples

mem_len = length(A_Data)

Sine0=[];% Base Frequency

Sine1=[];% Required Frequency

index=1;

N=32;% represend the frequency multiplier

for i=1:mem_len

Sine0 = [Sine0 A_Data(i)]; % get sin value, in range 0.0-1.0

if (index >= mem_len) %|| (index > mem_len)

index =1;

else

Sine1 = [Sine1 A_Data(index)];

index = index +N;

end

end

figure;

subplot(2,2,1);

plot(Sine0)

title('Original Sine Wave 1 Hz')

subplot(2,2,2);

plot(Sine1)

title('Sine Wave with Higher Frequency (5 Hz)')

%%

X_Sine0=fftshift(fft(Sine0,NFFT)); %compute DFT using FFT

fVals=(-NFFT/2:NFFT/2-1)/NFFT; %DFT Sample points

subplot(2,2,3);

% plot(fVals,abs(X_Sine0));

plot(fVals,abs(X_Sine0));

title('Double Sided FFT - with FFTShift');

xlabel('Normalized Frequency')

ylabel('DFT Values');

X_Sine1=fftshift(fft(Sine1,NFFT)); %compute DFT using FFT

fVals=(-NFFT/2:NFFT/2-1)/NFFT; %DFT Sample points

subplot(2,2,4);

plot(fVals,abs(X_Sine1));

title('Double Sided FFT - with FFTShift');

xlabel('Normalized Frequency')

ylabel('DFT Values');

%%

figure;

subplot(2,1,1);

plot(Sine1)

title('Sine Wave with Higher Frequency (5 Hz)')

X_Sine1=fftshift(fft(Sine1,NFFT)); %compute DFT using FFT

fVals=(-NFFT/2:NFFT/2-1)/NFFT; %DFT Sample points

subplot(2,1,2);

plot(fVals,abs(X_Sine1));

title('Double Sided FFT - with FFTShift');

xlabel('Normalized Frequency')

ylabel('DFT Values');

%%

X_Sine1=fft(Sine1,NFFT)/NFFT; %compute DFT using FFT --(/NFFT) scale factor for FFT

%fVals=(-NFFT/2:NFFT/2-1)/(NFFT); %DFT Sample points

FFT_imag = imag(X_Sine1)';

FFT_real = real(X_Sine1)';

figure;

plot(abs(FFT_imag));

figure;

plot(abs(FFT_real));

figure;

plot(abs(X_Sine1));

title('Double Sided FFT - with FFTShift');

xlabel('Normalized Frequency')

ylabel('DFT Values');

===========

The output of MATLAB - Last figure

The amplitude is 8000

The output of the Xilinx FPGA is:

The amplitude is 4000

Any explaination for this will be appreciated.

RSK
on 19 Apr 2018

Dr. Seis,

I think here and at other places where you have mentioned the word "amplitude", you meant the amplitude of the resulting transform after multiplying by dt, correct? (and that is not going to be the amplitude of the sine wave at that frequency and having that transform). However, if people are interested in the amplitude of sine wave that generated the transform, should option#2 in Chani's example be used(although that is not the real amplitude of the transform)? If I have something, say, machine bed(for the sake of example) vibrating at 1hz with unit magnitude and I am presented with the transform generated by "option#1" as the measured vibration, I should certainly not conclude that the machine bed is vibrating with the amplitude shown by "option#1". So I guess dividing by N (option#2) would tell me the magnitude with which the machine bed is vibrating. Is this understanding correct? Thanks

Bruno Luong
on 21 Sep 2018

Edited: Bruno Luong
on 21 Sep 2018

IMHO, if the input signal is x(t) is real, only about half of X(f) = fft(x)(f) is relevant (index 1:ceil((N+1)/2))) the other half is just conjugate of the flip of the first part, so the Parseval would be matching the integration of X(f)^2 on [0,1/(2*dt)] (up to Nyquist frequency) and one must multiply by an extra factor of abot sqrt(2) on top of dt on fft(x(t_i)) to get the Parseval equality up to Nyquist frequency (and not beyond that).

Some subtle consideration: if N is even, the frequency ceil((N-1)/2)*df count only once so it should not be multiplied by sqrt(2). The same for static term (f=0) for all the cases.

A factor of sqrt(2): it can make a bridge resists or fallen down.

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

Start Hunting!