CODEWITHSUNDEEP
matlabfft
Domi
Domi

Reputation: 198

Need help re: FFT output scaling

I'm currently studying the following book: "Fourier Transform Spectroscopy Instrumentation Engineering", by Vidi Saptari. My question is related to the code below, based on the code from the book, Appendix C. The code below computes the interferogram of 3 waves with wavenumbers [cm-1] 5000, 10000 and 15000, respectively, and than performs an FFT to retrieve the information. The unscaled output has a magnitude of 1600, instead of 1.

clear;
% Sampling clock signal generation
samp_period_nm = 632.8 / 4; % sampling period in nm. 632.8 is the HeNe laser's wavelength
samp_period = 1 * samp_period_nm * 10^-7; % sampling period in cm. 
scan_dist   = 0.1; % mirror scan distance in cm. 
no_elements = floor(scan_dist/samp_period);
x_samp = 0:samp_period:samp_period*no_elements; %Vector of clock signals in cm
xn_samp = x_samp .* (1 + rand(1, length(x_samp)));
v1 = 5000;
v2 = 10000;
v3 = 15000;
arg  = 4 * pi * x_samp;
y   = cos(arg*v1) + cos(arg*v2) + cos(arg*v3) ;
total_data = 2^18;
no_zero_fills=[total_data - length(y)];
zero_fills=zeros(1, no_zero_fills);
%triangular apodization
n_y = length(y);
points = 1:1:n_y;
tri = 1 - 1/(n_y) * points(1:n_y);
y = y.*tri; %dot product of interferogram with triangular apodization function
y = [y zero_fills];   %zero filling
% FFT operation
fft_y = fft(y);
% fft_y = fft_y / n_y;
% fft_y = fft_y * samp_period;
fft_y(1) = [];
n_fft=length(fft_y);
spec_y = abs(fft_y(1:n_fft/2)); %spectrum generation
nyquist = 1 / (samp_period * 4);
freq = (1:n_fft/2)/(n_fft/2)*nyquist; %frequency scale generation
figure();
plot(freq, spec_y);   % plot of spectrum vs wave number
xlabel('Wavenumber [cm-1]');
ylabel('Intesity [V]');

By multiplying the result of the fft (fft_y) with dt = samp_period, as suggested here, the peak is as 0.025.

Following the same link's second solution, by dividing fft_y by n_y (the length of y), the magnitude is 0.25.

Clearly, I'm doing something wrong. Any help is appreciated.

Thanks,

Upvotes: 1

Views: 6960

Answers (2)

Felix Bors
Felix Bors

Reputation: 21

if you want the same energy in the IFFT input and ouput you have to multiply the IFFT output (time signal) by sqrt(N) where N is the size of the transform.

here's the code :

same thing for the FFT (divide output by sqrt(N)); hope it helps Felix

N = 4096;
Freq = randn(N,1)+1j*randn(N,1);
Time = sqrt(N)*ifft(Freq,N);
FreqEn = sum(real(Freq).^2 + imag(Freq).^2); 
TimeEn = sum(real(Time).^2 + imag(Time).^2);

TimeEn/Freq

Upvotes: 0

Deve
Deve

Reputation: 4648

The only thing you're doing wrong here is expecting the peaks in the spectrum to be 1. According to Parseval's theorem of DFT the energy of the time domain signal is equal to the energy of the frequency domain signal divided by the lenght of the sequence N. You can check this in your example:

td_energy = sum( abs(y).^2 )
fd_energy = sum( abs(fft_y).^2 )
td_energy - fd_energy / length(y) % won't be exactly zero because you deleted the zero frequency bin.

So the peaks in your spectrum don't represent the amplitudes of cosine waves in the time domain but their energy. Also note at this point that the energy is lower than you might expect as you padded a lot of zeros.

In practice, the average power of a certain frequency is often of greater interest. Consider the following code example

t = linspace(-4*pi, 4*pi, 2^16);
N = length(t); % DFT length
y = cos(t); % single cosine wave
y_pow = sum( abs(y).^2 ) / N; % is 0.5
fft_y = fft(y);
fft_y_pow = (sum( abs(fft_y).^2 ) / N) /N; % is 0.5
figure; plot(abs(fft_y)./N);

The power is obtained by averaging the energy by the length of the sequence N. If you divide the spectrum by N you obtain the average power per frequency. In the above example you recognize a single peak with height 0.5 that represents the single cos wave of amplitude 1 (and hence power 0.5).

Personally, I prefer scaling MATLAB's FFT output by 1/sqrt(N) and its IFFT output by sqrt(N). In this way, the energy of the time and frequency domain sequence are always equal.

Upvotes: 2

Related Questions