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.
% 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.
No products are associated with this question.
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?
Here is a way to clean up the code and make it easier to understand and debug:
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)';
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;