function spectrum = pulsefft(index,pulsemodel)
%  returns the FFT of a decaying exponential pulse waveform
%      index = 0 for input waveform, ivar for output waveform ivar
%      pulsemodel.t0 = time when pulse begins (ns)
%      pulsemodel.tau = decay time of signal (ns)
%      pulsemodel.trise = rise time of leading edge (ns)
%      pulsemodel.period = total length of signal (ns)
%      pulsemodel.N = total number of samplings during period
%  To get the default model, invoke with empty arglist, as pulsefft(0,[])

  if isempty(pulsemodel)
    pulsemodel.period = 1e3;
    pulsemodel.t0 = 30;
    pulsemodel.tau = 2.7;
    pulsemodel.trise = 1;
    pulsemodel.N = 10000;
  end
  T = pulsemodel.period;
  N = pulsemodel.N;
  t0 = pulsemodel.t0;
  tr = pulsemodel.trise;
  tau = pulsemodel.tau;

  varnames;
  response(N) = 0;
  spectrum(N) = 0;
  gaussian(N) = 0;

  i = sqrt(-1);
  for n=1:N,
    omega = 2*pi*(n-1)/T;
    alpha = -1/tau-i*omega;
    spectrum(n) = (1-exp(alpha*T))/(1-exp(alpha*T/N));
    t = T*(n-1)/N;
    gaussian(n) = exp(-0.5*((t-t0)/tr)^2);
    if (n-1 > N/2)
      response(n) = conj(response(N-n+2));
    elseif (index ~= 0)
      freq = 1e9*omega/(2*pi)+1;
      response(n) = equations(index,iIin,freq);
    end
  end
  gaussianfft = fft(gaussian);
  spectrum = spectrum.*gaussianfft;

  Qpixel = 7.5e5*1.6e-19;
  PEmean = 1;
  Qmean = PEmean*Qpixel;
  Imean = Qmean/(T*1e-9);
  fnorm = Imean*(N/spectrum(1));
  spectrum = spectrum.*fnorm;
  if (index ~= 0)
    spectrum = spectrum.*response;
  end
