- 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
4 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
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 comentarios
Rick Rosson
el 1 de Sept. de 2012
I do not have a copy of the book you are using, so I cannot see the equations you mentioned or Appendix C.
Respuesta aceptada
Rick Rosson
el 31 de Ag. de 2012
Editada: Rick Rosson
el 31 de Ag. de 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
2 comentarios
Más respuestas (1)
Rick Rosson
el 1 de Sept. de 2012
Editada: Rick Rosson
el 2 de Sept. de 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));
Ver también
Categorías
Más información sobre Dates and Time en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!