Need help re: FFT output scaling

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

I have a few questions:

1. What is the unit of measure for v1, v2, and v3?
2. Why are you using 4*pi to compute arg?
3. Why is your sampling interval related to the wavelength of HeNe?
4. What do you expect for the value of the magnitudes?

Hey Rick,

Thanks for the reply. To answer your questions:

1. v1, v2, v3 are wavenumbers measured in cm-1

2. I'm using equation 3.22 provided in the book at page 25. According to it, the interference of a polychromatic source is

   I_ac (x) = sum(C(v_i)*cos(4*pi*v_i*x)), 

wavenumber v_i = v_min .. v_max C(v_i) is the instensity of each particular wavenumber v_i. In this example, to keep it simple, C(v_i) for each of the wavenumber is 1.

3. it's related to the wavelength of the HeNe laser as it's simulating the interferogram provided by a HeNe Fourier Transform Infra-Red (FTIR) spectrometer. Since this FTIR is using a HeNe laser, as far as I understand, the HeNe laser's wavelength is the minimum wavelength the radiation will contain. Anyone, please correct me if I misunderstood this. Nyquist criterion states that the sampling distance should be at least l_min / 4, where l_min is the minimum wavelength in the radiation.

4. By applying an FFT to the interferogram, we obtain the values of C(v_i), that is the intensity of each wavenumber v_i (according to the theory). Thus, I'm expecting to obtain peaks at 5000, 10000 and 15000 of magnitude 1. Not 0.25.

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?

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

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

Answer 1

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

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!

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 2

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

Thanks!