From: arun dixit on
HI FRIENDS,
i am using ifft to obtain inverse fourier transform of of a set of data obtained from a given expression.but i am unable to reconstruct the orignal result.my codes are following

function example_recons_givenprofile
clc;
clear all;
close all;
global f;
x_cord_original=[0:.04:4];


frequency=0:5e7:5e9; %frequency range
n=length(frequency);

%For computing gamma
gamaf=0;
xinv=[4:-.04:0];
for j=1:n
f=frequency(j);
[xinv_returned,gama]=ode45(@differential_recons_example,xinv,gamaf); % inverted x axis to make gamma final value as initial value input for ode45
d=length(gama);
gama_at_x_eq_0(j)=gama(d); %from inverse calculation from epsilon in ode 45
calculatedgamma(j)=1/(1+i*2*2*pi*f/3e8);
end

epsilon(1)=1; %continuus boundary with air
R01=(1-sqrt(epsilon(1)))/(1+sqrt(epsilon(1)));

gamanot=((gama_at_x_eq_0-R01)./(1-(R01.*gama_at_x_eq_0)));
gamahat=(gama_at_x_eq_0./(1-gama_at_x_eq_0.^2))-((1/3)*(gama_at_x_eq_0.^3));
k0=2*pi*frequency/3e8;
expected_gamahat=((4i*k0.*(k0*i+1)).^(-1))+((2*(i*k0+1)).^(-1))-((24*((i*k0+0.5).^3)).^(-1)); % value of expected gama obtain by solving manually
figure(2)
plot(frequency,gamahat,'r',frequency,expected_gamahat,'--')
figure(3)
plot(frequency,expected_gamahat,'b')

r=ifft(gamahat)
expected_inv_fourier_transform=(5e7)*r;
figure(4)
plot(frequency,(r),'b')
r_hat=inv_fourier_transform(expected_gamahat,frequency);
figure(5)
plot(frequency,r_hat,'r')

t(1)=0; % t0 =0
x(1)=t(1); %eq 14
g(1)=0; % g0=0
%no of points processed in t and n are same so tn can be replaced
%by sample location so g(t(n+1)) is eqvt to g(n+1)
expectedepsilon(1)=1;

for t_position_in_array=2:n
del_t=1/(n*5e7); % for ifft with N sample points and deltaf frequency spacing delta time =1/(N*deltaf)
t(t_position_in_array)=t(t_position_in_array-1)+del_t;
g(t_position_in_array)=g(t_position_in_array-1)-(2*del_t*((expected_inv_fourier_transform((t_position_in_array-1)))+(expected_inv_fourier_transform(t_position_in_array))));
epsilon(t_position_in_array)=epsilon(1)*exp((g(t_position_in_array)));
x(t_position_in_array)=(x(t_position_in_array-1)+(del_t/(sqrt(epsilon(t_position_in_array-1))+sqrt(epsilon(t_position_in_array)))));
end
expectedepsilon=(1+3*x_cord_original).^(-4/3);
figure(6)
plot(x,abs(epsilon),'r')%,x_cord_original/max(x_cord_original),abs(expectedepsilon),'b')
figure(7)
plot(x_cord_original/max(x_cord_original),abs(expectedepsilon),'b')
% % -------------------------------------------------------------------------

function gamadot = differential_recons_example(x,gama) %differential(x values, final value)
global f
k0=2*pi*f/3e8;
%epsilon_r=x+1;
epsilon_r = (1+3.*x).^(-4/3);
gamadot=((2*i)*sqrt(epsilon_r).*gama*k0)+((1-gama^2)./(4*epsilon_r)).*(-4*(1+3*x).^(-7/3));