- What is the unit of measure for v1, v2, and v3?
- Why are you using 4*pi to compute arg?
- Why is your sampling interval related to the wavelength of HeNe?
- What do you expect for the value of the magnitudes?
Need help re: FFT output scaling
3 views (last 30 days)
Show older comments
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 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.
Accepted Answer
Rick Rosson
on 31 Aug 2012
Edited: Rick Rosson
on 31 Aug 2012
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
More Answers (1)
Rick Rosson
on 1 Sep 2012
Edited: 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));
See Also
Categories
Find more on Dates and Time in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!