function PulseFFT

% Eric Lundquist
% PulseFFT

clear all
close all
clc

global cf
global Zo1
global Zo2
global Zo3

% Zo (Characteristic Impedance) Changes
Zo1 = 50 ; % ohms
Zo2 = Zo1 ; % ohms
Zo3 = 125 ;

pulse = pulsegen;
N = length(pulse);
t_total = 1e-3;
dt = t_total/(N-1);
t = 0:dt:t_total ;

figure
plot(t,pulse)
title('Initial Pulse Shape')
xlabel('time (s)')
ylabel('(V)')

P = fftshift(fft(pulse))/N;
M=abs(P); % magnitude
A=angle(P); % phase (angle)
F = (0:N-1)/(2*dt); % frequency
cf = F(round(N/2));

figure
subplot(2,1,1)
plot(F,M)
ylabel('Magnitude')
subplot(2,1,2)
plot(F,A)
ylabel('Phase')

R=zeros(1,length(M));
for n=1:length(M)
    if M(n)>0.01
        R(n) = 2*PulseFDTD(F(n))*M(n)*exp(1i*A(n));
    end
end

Rmag = abs(R);
Rphase = angle(R);

analytic_refl_coeff = (Zo3-Zo1)/(Zo3+Zo1)
simulated_refl_coeff = max(Rmag)
simulation_error = (analytic_refl_coeff-simulated_refl_coeff)/analytic_refl_coeff*100

figure
plot(F,Rmag)
xlabel('Frequency (Hz)')
ylabel('Reflection Coefficient (| \Gamma|)')
title(['Expected | \Gamma_m_a_x| = ', num2str(analytic_refl_coeff)]) % applies for 2 regions only

IFFTR = abs(ifft(R))*2*N;

figure
plot(t,IFFTR)
title('IFFT Pulse Shape')
xlabel('time (s)')
ylabel('(V)')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [pulse] = pulsegen
% Eric Lundquist
% Pulse Generator

pulseplots = 0 ; % enter 1 if pulse-shape plots are desired

plength = 100;
p = zeros(1,plength);
p(25:75)=1; % square pulse

noise = randn(1,plength)/plength;
if pulseplots==1
    figure
    plot(noise)
    title('Noise')
end

pn=p+noise; % pulse + noise
if pulseplots==1
    figure
    plot(pn)
    title('Pulse + Noise')
end

n=1; % order of filter i.e., 1 realistic, 2 or 3 for more distortion
Wn=0.1;
[b,a] = butter(n,Wn,'low');
if pulseplots==1
    freqz(b,a)
end
pnf = filter(b,a,pn); % pulse + noise + filtered
if pulseplots==1
    figure
    plot(pnf)
    title('Lowpass Filtered Signal + Noise')
end

pnfn = pnf+noise/5; % pulse + noise + filtered + more noise
if pulseplots==1
    figure
    plot(pnfn)
    title('Lowpass Filtered Signal + Noise')
end

pulse=pnfn;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [R] = PulseFDTD(fp)
% Eric J. Lundquist
% FDTD Connector Characteristic Impedance / Reflection Simulation
% 3 Regions

global cf
global Zo1
global Zo2
global Zo3

medium = 'air' ; % enter air, teflon, or sea water
cable = 3 ; % 1 is Measured 3300 Cable (predetermined Zo), 2 is RG58 Coaxial Cable (predetermined Zo), 3 is Lossless Line, 4 is Lossy Line with user-controlled Zo
src = 'gaussian' ; % enter desired voltage source type, choose from pulse, sine, rsine, spulse, gaussian
pulsetype = 1 ; % if Gaussian pulse selected, 1 is raised only, 2 is monopulse, and 3 is Gaussian-modulated sinusoidal pulse

Zo = [Zo1,Zo2,Zo3];

% General Parameters
f = cf*100;
w = 2*pi*f;
eo = 8.854e-12;
u0 = 4e-7*pi;

% Copper Conductor
sigma_c =5.8e7;
u_c = u0;

% Parameters for Surrounding Medium
switch medium
    case 'air'
        sigma = 0 ;
        E_r = 1 ;
        u_r = u0 ;
    case 'teflon'
        sigma = 0 ;
        E_r = 2.1 ;
        u_r = u0 ;
    case 'sea water'
        sigma = 5 ;
        E_r = 81 ;
        u_r = u0 ;
    otherwise
        warning('Surrounding medium has not been specified! Air selected by default.')
        sigma = 0 ;
        E_r = 1 ;
        u_r = u0 ;
end

maxZ = 300 ; % cells in the transmission line
maxT = 10*maxZ ; % maximum number of time steps (time=2*length before wave hits the end of grid)

% Compute Constants
lossyon = 0 ; % lossyon=0 means the line is lossless, otherwise it is lossy  
switch cable
    case 1 % Measured 3300 Cable
        R0 = 3.6996 ;
        L0 = 2.7630e-007 ;
        G0 = 4.5602e-004 ;
        C0 = 9.0846e-011 ;
        lossyon=1;
    case 2 % RG58 Coaxial Cable
        a = 0.445e-3 ;
        b = 1.765e-3 ;
        rs = sqrt(pi*f*u_c/sigma_c) ;
        R0 = (rs/(2*pi))*(1/a + 1/b) ;
        L0 = (u_r/(2*pi))*log(b/a) ;
        G0 = (2*pi*sigma)/log(b/a) ;
        C0 = (2*pi*E_r*eo)/log(b/a) ;
        lossyon=1;
    case 3 % Lossless Line
        R0 = 0 ;
        L0 = 1e-6 ;
        G0 = 0 ;
        C0 = L0/(Zo1^2) ;
        C0_max = L0/(max(Zo)^2);
    case 4 % Lossy Line with user-determined Zo
        R0 = 4 ;
        L0 = 1e-6 ;
        G0 = 4e-4 ;
        C0 = L0/(Zo1^2) ;
        C0_max = L0/(max(Zo)^2);
    otherwise
        warning('No cable selected! Lossless line selected by default.')
        R0 = 0 ;
        L0 = 1e-6 ;
        G0 = 0 ;
        C0 = L0/(Zo1^2) ; % 2500 = (50 ohms)^2
        C0_max = L0/(max(Zo)^2);
end
R(1:maxZ) = R0 ;
L(1:maxZ) = L0 ;
G(1:maxZ) = G0 ;
C(1:maxZ) = C0 ;

% Constants computed in previous section
gamma = sqrt((R+1i.*w.*L).*(G+1i.*w.*C)) ;
alpha = real(gamma) ;
beta = imag(gamma) ;
vp = w./beta ;
lambda = vp./f ;

gamma_max = sqrt((R0+1i*w*L0)*(G0+1i*w*C0_max)) ;
beta_max = imag(gamma_max) ;
vp_max = w/beta_max ;
lambda_max = vp_max/f ;
dz = lambda_max/20 ;
dt = 0.5*dz/vp_max ;
length_total = dz*maxZ ;

% Additional Constants
A = -1./(dz.*(R./2 + L./dt));
B = -1.*((R./2 - L./dt)./(R./2 + L./dt));
D = -1./(dz.*(G./2 + C./dt));
E = -1.*(G./2 - C./dt)./(G./2+C./dt);

% Pulse Parameters
ramp = 1 ; % initialize parameter for ramped sine case
switch src
    case 'gaussian'
    t = -dt*maxT/2:dt:dt*maxT/2;
    switch pulsetype
        case 1 % raised only
            width = length(gmonopuls(t,fp));
            gp = exp((-([-width:width]).^2)/(2*width));
        case 2 % monopulse
            gp = gmonopuls(t,fp);
        case 3 % Gaussian-modulated sinusoidal pulse
            bw = .5 ; % bandwidth (i.e., 0.6=60%)
            gp = gauspuls(t,fp,bw);
        otherwise
            warning('Gaussian Pulse Type Unspecified! Monopulse selected by default.')
            gp = gmonopuls(t,fp);
    end
    x = length(gp);
    s = round(x*0.4); % shorten Guassian pulse by 40%
    g = gp(s:end);
    gp = zeros(1,x);
    gp(1:x-s+1)=g;
end

% Connector Parameters
cc = 2 ; % number of cells from center point of Region 2
connector_length = dz*(2*cc+1);
ccc = cc*5 ; % number of cells from center where reflection measurements take place

% Initialize voltage and current values of transmission line
V = zeros(1,maxZ) ; 
I = zeros(1,maxZ) ;
Vo = zeros(1,maxT); % original signal
Vr = zeros(1,maxT); % signal with reflections

% FDTD Loop 1 (without reflection)
for n=1:maxT
    switch src
        case 'pulse'
            if n<=25 % length of pulse
                V(1)=1;
            else
                V(1)=-1;
            end
        case 'sine'
            V(1) = sin(2*pi*fp*n*dt); % sine wave source
        case 'rsine'
            if n*dt < 1/f/4
                V(1) = sin(2*pi*fp*n*dt);
            else
                V(1) = sin(2*pi*fp*n*dt)*ramp;
                if ramp>0
                    ramp = ramp-.01;
                end
            end
        case 'spulse'
            if n*dt < 1/f/4
                V(1) = sin(2*pi*fp*n*dt);
            else
                V(1)=0;
            end
        case 'gaussian'
            V(1) = gp(n);
        otherwise
            error('Voltage Source Type Unspecified!')
    end
    for k=2:maxZ; % find voltage everywhere on the line
        V(k)= D(k)*(I(k)-I(k-1))+E(k)*V(k);
    end
    for k=1:maxZ-1; % find current everywhere on the line
        I(k)= A(k)*(V(k+1)-V(k))+B(k)*I(k);
    end
    Vo(n)=V(round(maxZ/2)-ccc);
    I(maxZ-1) = V(maxZ-1)/Zo1 ; % absorbing boundary (not needed if simulation stops before hitting end; maxT=2*maxZ)
    movieon = 0 ;
    if movieon==1
        x = 50 ; % movie rate, ie. x=2 causes movie to play every other frame 
        if logical(x*floor(n/x)==n) == 1 % test for xth values of n (shorter movie times, plays only every xth calculation)
            figure(1)
            plot(V)
            title('Original Signal')
            xlabel('Distance')
            ylabel('Voltage')
            axis([0 maxZ -2 2]) % control the axis for uniform pictures
            pause(.0001) % give the program time to plot to screen
        end
    end
end
if movieon==1 % option available for visual validation purposes in programming
    Vo(round(maxT/2):end)=0; % clear out slight reflections caused by imperfect absorbing boundaries
    figure(3)
    plot(Vo)
    title('Pure Signal')
end

% Initialize voltage and current values of transmission line
V = zeros(1,maxZ) ; 
I = zeros(1,maxZ) ;

% Compute Constants for Connector Impendances (Region 2)
loadC = round(maxZ/2)-cc:round(maxZ/2)+cc; % connector region
L(loadC)=1e-6;
C(loadC)=(1e-6)/(Zo2^2);
if lossyon==0
    R(loadC)=0;
    G(loadC)=0;
end

% Compute Constants for Cable 2 (Region 3)
r3=round(maxZ/2)+cc+1; % starting point after connector
L(r3:maxZ) = 1e-6 ;
C(r3:maxZ) = (1e-6)/(Zo3^2) ;
if lossyon==0
    R(r3:maxZ) = 0 ;
    G(r3:maxZ) = 0 ;
end

% Constants computed in previous section
gamma = sqrt((R+j.*w.*L).*(G+j.*w.*C)) ;
alpha = real(gamma) ;
beta = imag(gamma) ;
vp = w./beta ;
lambda = vp./f ;

% Compute Constants for End Terminal
loadZ = maxZ-1;
R(loadZ)=Zo3;
R(2)=Zo1;

% Additional Constants
A = -1./(dz.*(R./2 + L./dt));
B = -1.*((R./2 - L./dt)./(R./2 + L./dt));
D = -1./(dz.*(G./2 + C./dt));
E = -1.*(G./2 - C./dt)./(G./2+C./dt);

% FDTD Loop 2 (with reflections)
for n=1:maxT
    switch src
        case 'pulse'
            if n<=25 % length of pulse
                V(1)=1;
            else
                V(1)=-1;
            end
        case 'sine'
            V(1) = sin(2*pi*fp*n*dt); % sine wave source
        case 'rsine'
            if n*dt < 1/f/4
                V(1) = sin(2*pi*fp*n*dt);
            else
                V(1) = sin(2*pi*fp*n*dt)*ramp;
                if ramp>0
                    ramp = ramp-.01;
                end
            end
        case 'spulse'
            if n*dt < 1/fp/4
                V(1) = sin(2*pi*fp*n*dt);
            else
                V(1)=0;
            end
        case 'gaussian'
            V(1) = gp(n);
        otherwise
            error('Voltage Source Type Unspecified!')
    end
    for k=2:maxZ; % find voltage everywhere on the line
		V(k)= D(k)*(I(k)-I(k-1))+E(k)*V(k);
    end
    for k=1:maxZ-1; % find current everywhere on the line
		I(k)= A(k)*(V(k+1)-V(k))+B(k)*I(k);
    end
    I(loadZ) = V(loadZ)/Zo3 ; % right absorbing boundary
    Vr(n)=V(round(maxZ/2)-ccc);
    movieon = 0 ;
    if movieon==1
        x = 20 ; % movie rate, ie. x=2 causes movie to play every other frame 
        if logical(x*floor(n/x)==n) == 1 % test for xth values of n (shorter movie times, plays only every xth calculation)
            figure(2)
            subplot(2,1,1)
            plot(V)
            title('Signal through Varying Cables and Connectors')
            xlabel('Distance')
            ylabel('Voltage')
            axis([0 maxZ -2 2]) % control the axis for uniform pictures
            pause(.0001) % give the program time to plot to screen
            subplot(2,1,2)
            bar = zeros(1,maxT);
            bar(1,1:n)=1;
            status = sprintf(['\n\nSimulation is ',num2str(round(n/maxT*100)),'%% Complete']);
            time_status = sprintf(['\n\nTime Step ',num2str(n),' of ',num2str(maxT),'\n']);
            imagesc(bar)
            colormap(winter)
            daspect([8,100/maxT,1]) % x/y aspect ratio
            axis off
            title([status,time_status])
        end
    end
end

Vrefl = Vr-Vo;
ReflCo = max(Vrefl); % does not need to be divided by max(Vo) because max(Vo)=1 for this source

if movieon==1
    figure(4)
    plot(Vrefl)
    title('Signal Reflections')
    xlabel('time steps')
    ylabel('Voltage')
end

R = max(Vrefl);

if Zo2>=Zo1
    Simulated = ReflCo;
else
    Simulated = -ReflCo;
end
Analytic = (Zo2-Zo1)/(Zo2+Zo1); % applies for 2 regions