# Need help re: FFT output scaling

Asked by Domi on 31 Aug 2012
Latest activity Commented on by Domi on 4 Sep 2012

Hello,

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 numbers [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 http://www.mathworks.com/matlabcentral/answers/15770-scaling-the-fft-and-the-ifft, the peak is as 0.025.

Following http://www.mathworks.com/matlabcentral/answers/15770-scaling-the-fft-and-the-ifft'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, Domi

## 5 Comments

Rick Rosson on 31 Aug 2012

In the following line of code:

`   y = y.*tri; %dot product of interferogram with triangular apodization function`

the comment claims that the right-hand expression computes the "dot product" of y and tri. I am not sure if your intent was to compute a dot product in the formal sense of the term (as in the dot product of one vector with another), or if your intent was to compute the element-wise multiplication of the two vectors, which is what the code actually does. I assume it is the latter, but I just want to make sure.

Also, why are you zero-padding the data to 2^18 elements before computing the FFT? Why not just compute the FFT without zero-padding?

Domi on 31 Aug 2012

Rick, you are correct, it actually refers to element-wise multiplication. I assume the author was referring to the .* operation. This code is based on the code in Appendix C provided by the author.

In the book I've mentioned, and in others I've studied, they recommend to do zero-padding to smooth out the spectrum, as the effect will be similar to interpolation. I'm not sure why pad till 2^18 though. Seems too much. I get the same data with 2^13.

Computing the FFT without zero-padding has a negative effect though, as the magnitude is decreasing the bigger the wavenumber. If with zero padding all the peaks are at 0.25, without padding they are 0.244, 0.229 and 0.205 for wavenumbers 5000, 10000 and 15000 [cm-1], respectively

Rick Rosson on 1 Sep 2012

I do not have a copy of the book you are using, so I cannot see the equations you mentioned or Appendix C.

## Products

No products are associated with this question.

## 2 Answers

Answer by Rick Rosson on 31 Aug 2012
Edited by Rick Rosson on 31 Aug 2012
Accepted answer

I think I know why you are getting 0.25 instead of 1.00 for the magnitudes. There are two separate reasons, each one of which is causing an error by a factor of 2x in the scaling, for a total error of 4x.

The first error is that you are looking at only half of the spectrum (the positive frequencies). You are "throwing away" the negative frequencies, which is fine, but if you want to have the "correct" magnitude, you have to adjust for that fact by multiplying the spectrum by a factor of 2.

The second is more subtle, and perhaps speculative on my part. The vector tri is a linear function that represents a triangular window that goes from a value of 1 on the left-hand-side of the axis to a value of 0 on the right-hand side. Since it is a triangle, it is easy to see that the average value of this triangular window is exactly 0.5. So, by multiplying the source signal by this trangular window, you have effectively cut the average amplitude by 0.5. As a result, the magnitude of the spectral lines will also be 1/2 of what they would otherwise be. So here again, you have to multiply the spectrum by another factor of 2 in order to adjust for this effect.

So I think the correct scaling would be:

```    % FFT operation
fft_y = 4 * fft(y) / n_y;```

Does that make sense?

Rick

## 2 Comments

Domi on 31 Aug 2012

You're right about scaling with 2x, due to using only half of the spectrum, but I'm not sure about the reason for the second scaling with 2x -- the apodization.

I have another input with 300+ wavenumbers, but if I use the scaling you mention, i don't get the same result. Let me check this again though.

Thanks!

Domi on 31 Aug 2012

Ok, so modified the code from above, and this time the polychromatic source contains 350 wavenumbers. Again, I assume all of the wavenumbers have an intensity of 1.

In this case I expect the output to be a horizontal line at y=1. Instead I get a bathtub curve, with most of it around 0.58.

```   wavenumbers = linspace(600, 2000, 350);
samp_period = 1/max(wavenumbers)/4/4 ;% Nyquist criteria. An additional factor of 4.
no_elements = 4096;
x_samp = 0:samp_period:samp_period*no_elements; %Vector of clock signals in cm```
```   y = 0;
for index = 1:length(wavenumbers)
y = y + cos(2 * pi * 2 * wavenumbers(index) * x_samp); % interferogram with phase error
end```
```   total_data = 2^13;
no_zero_fills=[total_data - length(y)];
zero_fills=zeros(1, no_zero_fills);```
```   figure();
plot(x_samp,y);
title('Interferogram');```
```   %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) * 4 / length(y);
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```
```   for i=1:1:length(freq)
if freq(i)>wavenumbers(1)
startIndex = i;
break;
end
end```
```   for i=1:1:length(freq)
if freq(i)>wavenumbers(end)
endIndex = i;
break;
end
end```
```   if endIndex < startIndex
spec_y = spec_y(endIndex:startIndex);
freq = freq(endIndex:startIndex);
else
spec_y = spec_y(startIndex:endIndex);
freq = freq(startIndex:endIndex);
end```
```   figure();
plot(freq, spec_y);   % plot of spectrum vs wave number
xlabel('Wavenumber [cm-1]');
ylabel('Intesity [V]');
legend('Received spectrum');```
Answer by Rick Rosson on 1 Sep 2012
Edited by Rick Rosson on 2 Sep 2012

Here is a way to clean up the code and make it easier to understand and debug:

`   %% Initialize`
```   close all;
clear all;
clc;```
`   %% Spatial domain`
```   % HeNe laser wavelength (in cm):
HeNe = 632.8e-7;```
```   % Spatial increment (in cm per sample):
dx = HeNe/4;```
```   % Number of samples:
N = 4096; ```
```   % Spatial domain (in cm):
x = dx*(0:N-1)';```
`   %% Signal`
```   Ac = [    1 ;     1 ;     1 ];
Vc = [ 5000 ; 10000 ; 15000 ];```
`   y = cos(2*pi*x*Vc')*Ac;`
```   % Apodization:
win = 1 - 1/N * (0:N-1)';
y = win .* y;```
```   % Zero-padding:
M = 2^18;
y = [ y ; zeros(M-N,1) ];```
`   %% Wavenumber domain:`
```   % Sampling rate (in samples per cm):
Vs = 1/dx;```
```   % Wavenumber increment (in waves per cm):
dV = Vs/M;```
```   % Wavenumber domain (in waves per cm):
v = (-Vs/2:dV:Vs/2-dV)';```
`   %% Fourier transform`
`   Y = 4*fftshift(fft(y))/N;`
```   figure;
plot(v,abs(Y));```

## 1 Comment

Domi on 4 Sep 2012

Thanks!

Contact us