Appendix C

Matlab codes

This appendix is a compilation of all the codes that were written in Matlab for the analysis
of raw data measured with the homemade Fourier transform spectrometer. The functions are
introduced in alphabetical order.

C.1 Function centerTheBurst.m

functi on [ i n t e r f e r o g r a m c e n t e r e d ] = centerTheBurst ( interferogram )

% DESCRIPTION : c e n t e r T h e B u r s t i s a f u n c t i o n t h a t f i n d s t h e
% c e n t e r b u r s t o f t h e i n t e r f e r o g r a m and p l a c e s
% i t a t t h e middle o f t h e data .
%
% INPUT : − i n t e r f e r o g r a m i s a v e c t o r c o n t a i n i n g t h e v a l u e s o f t h e
% measured i n t e r f e r o g r a m .
%
% OUTPUT : − i n t e r f e r o g r a m c e n t e r e d i s a v e c t o r c o n t a i n i n g t h e v a l u e s
% o f t h e c e n t e r e d i n t e r f e r o g r a m . The v e c t o r l e n g t h N i s
% t h e l a r g e s t power o f two t h a t i s s m a l l e r than o r e q u a l
% t o t h e l e n g t h o f i n t e r f e r o g r a m . The maximum i s p l a c e d a t
% i n d e x 2 ˆ (N−1) o r 2 ˆ (N−1)+1 depending on t h e w e i g h t o f
% t h e 64 f i r s t v a l u e s on both s i d e s o f t h e i n t e r f e r o g r a m
% maximum .

[ ˜ , i n d e x ] = max( i n t e r f e r o g r a m ) ;
h a l f l e n g t h n e w = pow2 ( nextpow2 ( min ( index , l e n g t h ( i n t e r f e r o g r a m ) −(
i n d e x +1) ) ) ) / 2 ;
i n t e r f e r o g r a m s h o r t = i n t e r f e r o g r a m ( index −64: i n d e x +64) ;

weight = ( −64:1:64) ’ ;
c e n t e r g r a v i t y = sum ( abs ( i n t e r f e r o g r a m s h o r t ) . ∗ w e i g h t ) ;

i f center gravity > 0

i n t e r f e r o g r a m c e n t e r e d = i n t e r f e r o g r a m ( index−h a l f l e n g t h n e w :
i n d e x+h a l f l e n g t h n e w −1) ;
e l s e i f center gravity < 0
i n t e r f e r o g r a m c e n t e r e d = i n t e r f e r o g r a m ( index−h a l f l e n g t h n e w +1:
i n d e x+h a l f l e n g t h n e w ) ;

73
end

end

C.2 Function filterTheReferenceSignal.m

function s i g n a l f i l t e r e d = filterTheReferenceSignal ( signal )

% D e s c r i p t i o n : This f u n c t i o n a l l o w s f o r d i g i t i l a l l y f i l t e r i n g
% t h e s i g n a l u s i n g a Butterworth f i l t e r ( o r d e r 2 0 ) .
%
% INPUT : − s i g n a l i s t h e He : Ne s i g n a l t o be f i l t e r e d
%
% OUTPUT : − s i g n a l f i l t e r e d i s t h e f i l t e r e d s i g n a l . I t has t h e same
% l e n g t h than s i g n a l .

[ b , a ] = b u t t e r ( 2 0 , [ 0 . 1 0 . 9 ] , ’ bandpass ’ ) ;
s i g n a l f i l t e r e d = f i l t e r (b , a , s ign al ) ;
s i g n a l f i l t e r e d = s i g n a l f i l t e r e d ( 2 1 : end ) ;
end

C.3 Function filterTheSignal.m

f u n c t i o n s i g n a l f i l t e r e d = f i l t e r T h e S i g n a l ( s i g n a l , wavelength max ,
wavelength min )
% DESCRIPTION : This f u n c t i o n d i g i t i l a l l y f i l t e r s t h e s i g n a l
% u s i n g a Butterworth f i l t e r ( o r d e r 5 ) .
%
% INPUTS : − s i g n a l i s t h e s i g n a l t o be f i l t e r e d
% − wavelength max i s t h e maximum wavelength i n nm o f t h e
% bandpass f i l t e r
% − wavelength min i s t h e minimum wavelength i n nm o f t h e
% bandpass f i l t e r
%
% OUTPUT : − s i g n a l f i l t e r e d i s t h e f i l t e r e d s i g n a l . I t has t h e same
% l e n g t h than s i g n a l .

% cutoff frequencies
w1 = 6 3 2 . 8 / wavelength max ; % low f r e q u e n c y
w2 = 6 3 2 . 8 / wavelength min ; % h i g h f r e q u e n c y

[ b , a ] = b u t t e r ( 5 , [ w1 w2 ] , ’ bandpass ’ ) ;
f v t o o l ( b , a ) ; % check t h e magnitude r e s p o n s e o f t h e f i l t e r

s i g n a l f i l t e r e d = f i l t e r (b , a , s ign al ) ;
s i g n a l f i l t e r e d = s i g n a l f i l t e r e d ( 2 1 : end ) ;
end

74
C.4 Function findTheZeros.m

function [ index z e r o c r o s s i n g ] = findTheZeros ( position , s i g n a l )

% DESCRIPTION : Compute an e s t i m a t i o n o f t h e z e r o −c r o s s i n g l o c a t i o n s
% o f t h e r e f e r e n c e l a s e r . They a r e o b t a i n e d by f i t t i n g
% t h e s i g n a l on an even number o f p o i n t s on both s i d e s
% o f t h e z e r o −c r o s s i n g and by f i n d i n g t h e r o o t s o f t h i s
% polynomial .
%
% INPUTS : − p o s i t i o n i s an o r d e r e d v e c t o r c o n t a i n i n g t h e i n d e x
% of the values in s i g n a l
% − s i g n a l corresponds to the i n t e r f e r o g r a m of the
% r e f e r e n c e l a s e r which i s sampled a t c o n s t a n t
% i n t e r v a l s o f t i m e s . ( The dc component has a l r e a d y
% been removed , i . e . t h e mean o f t h e s i g n a l i s
% equal to zero )
%
% OUTPUTS : − i n d e x z e r o c r o s s i n g i s a vector containing the
% i n d i c e s a t which s i g n a l i s v a n i s h i n g . These
% i n d i c e s a r e r e a l numbers a t which t h e i n t e r f e r o g r a m
% o f t h e o t h e r s o u r c e must be i n t e r p o l a t e d .

warning ( ’ o f f ’ , ’ a l l ’ ) ;

% I n t e r p o l a t i o n o f t h e s i g n a l . The amount o f data i s i n c r e a s e d

% by a f a c t o r 8 ( G r i f f i t h s ( 2 0 0 2 ) ) .
s i g n a l = i n t e r p ( s i g n a l , 8) ;

% The f i r s t 8 v a l u e s a r e not c o n s i d e r a t e d b e c a u s e they were checked

% i n t h e p r e v i o u s s e t o f data on which t h e f u n c t i o n has a l r e a d y been
% applied .
% The 7 l a s t v a l u e s a r e not c o n s i d e r a t e d s i n c e they o n l y c o n t a i n
% e x t r a p o l a t e d v a l u e s r e s u l t i n g from t h e a p p l i c a t i o n o f
% the f u n c t i o n i n t e r p .
l o c a t i o n = f i n d ( d i f f ( s i g n ( s i g n a l ( 9 : end −7) ) ) ) + 8 ;

% number o f p o i n t s f o r t h e c u r v e f i t t i n g .
% / ! \ data can o n l y t a k e even v a l u e s t h a t a r e s m a l l e r o r e q u a l t o 16
data = 1 0 ;

% I n i t i a l i z a t i o n o f v a r i a b l e s and v e c t o r s n e c e s s a r y f o r t h e l o o p :
k = 0;
index zero crossing = zeros ( length ( signal ) ,1) ;

for j = 1: length ( location )

% D e f i n i t i o n of the p o s i t i o n vector f o r the cubic polynomial f i t .
% An even number o f p o i n t s on both s i d e s o f t h e z e r o −c r o s s i n g i s
% considered .
% −The p o s i t i o n v e c t o r i s m u l t i p l i e d by 8 u s i n g t h e i n t e r p
% f u n c t i o n , t h e l a s t i n d e x b e f o r e t h e z e r o −c r o s s i n g i s g i v e n
% by 1+ f l o o r ( l o c a t i o n ( j ) / 8 ) .

75
 % −There a r e ( ( mod( l o c a t i o n ( j ) , 8 ) −1) i n t e r v a l s o f 0 . 1 2 5 between
% t h e l a s t p o s i t i o n p o i n t and t h e p o i n t b e f o r e t h e z e r o −c r o s s i n g
p o s i t i o n 2 = ( p o s i t i o n (1+ f l o o r ( l o c a t i o n ( j ) / 8 ) ) + 0 . 1 2 5 ∗ ( ( mod(
l o c a t i o n ( j ) , 8 ) −1) )+ 0.125∗(( − data /2+1) : ( data / 2 ) ) ) ’ ;

% D e f i n i t i o n of the s i g n a l vector f o r the cubic polynomial .

% An even number o f p o i n t s on both s i d e s o f t h e z e r o −c r o s s i n g
% are considered .
s i g n a l 2 = s i g n a l ( ( l o c a t i o n ( j )−data /2+1) : ( l o c a t i o n ( j )+data / 2 ) ) ;

% Determine t h e p o l y n o m i a l t h a t b e s t f i t s t h e data
polynomial = p o l y f i t ( position2 , signal2 , 3 ) ;

% Find t h e 3 r o o t s o f t h e c u b i c p o l y n o m i a l .
ra c i n e = roots ( polynomial ) ;
% The r o o t t h a t i s c o n s i d e r e d a s a z e r o −c r o s s i n g i s t h e one which
% i s t h e n e a r e s t t o t h e two v a l u e s l o c a t e d a t each s i d e o f t h e
% d e t e c t e d z e r o −c r o s s i n g .
[ ˜ , i n d e x ] = min ( ( r a c i n e −p o s i t i o n 2 ( data / 2 ) ) . ˆ 2 + ( r a c i n e −p o s i t i o n 2
( data /2+1) ) . ˆ 2 ) ;

% The z e r o −c r o s s i n g i s saved i n t h e output v e c t o r

k = k+1;
i n d e x z e r o c r o s s i n g ( k ) = r a c i n e ( index ) ;

end

% S u p p r e s s t h e z e r o s r e m a i n i n g a t t h e end o f t h e v e c t o r . A v a l u e was
% assigned to the k f i r s t po in ts .
index zero crossing = index zero crossing (1: k) ;

warning ( ’ on ’ , ’ a l l ’ ) ;
end

C.5 Function getTheSpectrum.m

f u n c t i o n [ wavenumber , spectrum ] = getTheSpectrum ( i n t e r f e r o g r a m )

% DESCRIPTION : This f u n c t i o n computes t h e n o r m a l i z e d i r r a d i a n c e from
% t h e i n t e r f e r o g r a m . The phase i s c o r r e c t e d by
% m u l t i p l y i n g t h e spectrum by exp(− phase ) .
%
% INPUT : −i n t e r f e r o g r a m i s a v e c t o r c o n t a i n i n g t h e v a l u e s o f t h e
% c e n t e r e d i n t e r f e r o g r a m , OPD i s c o n s t a n t between two
% c o n s e c u t i v e data p o i n t s .
%
% OUTPUT : −wavenumber i s a v e c t o r c o n t a i n i n g t h e wavenumber v a l u e s
% a s s o c i a t e d t o spectrum
% −spectrum i s t h e n o r m a l i z e d spectrum a s s o c i a t e d t o t h e
% interferogram .

76
% Fourier transform
fourier transform = f f t ( interferogram ) ;

% S u p p r e s s i o n o f t h e r e m a i n i n g dc term
fourier transform (1) = 0;

% phase c o r r e c t i o n
j j = s q r t ( −1) ;
f o u r i e r t r a n s f o r m = f o u r i e r t r a n s f o r m . ∗ exp(− j j ∗ a n g l e (
fourier transform ) ) ;

% d e f i n i t i o n o f t h e wavenumber a x i s
N = length ( fourier transform ) ;
c s t = 2 / 6 3 2 . 8 e −7;
wavenumber = ( c s t ∗ ( 0 : (N/2−1) ) /N) ’ ;

spectrum = r e a l ( f o u r i e r t r a n s f o r m ( 1 : (N/ 2 ) ) ) /max( r e a l (

f o u r i e r t r a n s f o r m ( 1 : (N/ 2 ) ) ) ) ;

end

77