# Problems Using FFT to Compute Impedance in a Model Neuron

I'm a neuroscientist currently investigating the resonance properties of a single neuron model that a colleague and I have constructed. The language we code in is Julia, which I hope is similar enough to Python for most people to understand. I've run in to a very particular problem, so bear with me please.

In order to investigate the resonance properties of the model neuron, I am applying an exponentially swept sinusoid current to it and analyzing its voltage response. In my field, such a current is called a ZAP current, though a more familiar term is probably Chirp signal. Just so there is no confusion, the current is generated from the following function \begin{gather*} ZAP(t) = I_{app}sin(2\pi f(t)t)\\f(t) = (0.001)\frac{f_{min}}{Lt}(e^{Lt}-1)\\ L = \frac{\ln{\frac{f_{max}}{f_{min}}}}{t_{max}} \end{gather*} where $t_{max}$ is the time duration of the simulation, $f_{min}$ and $f_{max}$ are the starting and ending frequencies of the sweep, and $I_{app}$ is the magnitude of the sinusoid. The voltage response of the model neuron to this input looks like this:

Clearly, there is a resonance in the model at whatever the input frequency is at around 50,000 ms. This means there should be a peak in the impedance of the model at this frequency. The impedance is calculated as the ratio of the discrete Fourier transform of the output voltage to the input current $$Z = \frac{FFT(V(t))}{FFT(ZAP(t))}$$ I then plot the magnitude of the impedance against the natural logarithm of the frequency. I expect to get something like the one on the left, but instead I obtain the one on the right:

As I vary the parameter $E_K$, the resonance frequency changes, and I expect to see this reflected in the impedance-frequency plot, but every one of the impedance-frequency plots looks almost identical to the one shown above despite the voltage responses looking quite different.

I have spent hours going through my code line by line and experimenting with different modifications to figure out where I'm going wrong, but have not been able to produce the desired output. Here is the stripped down code:

## Time Parameters
dt      = 0.01         # ms
tmax    = 100000        # ms
tspan   = convert(Int64,tmax/dt)
t       = linspace(dt,tmax,tspan)

## Simulation and Analysis
# magnitude of sinusoid
Iapp = 1.0
# minimum frequency
f0 = 0.1
# maximum frequency
ff = 15.0
# create an array of the frequency sweep values for plotting
Zfreq = ZAPfreq(f0, ff)
# optimize Fast Fourier Transform for an array of size and type ZAP()
P = plan_fft(ZAP(Iapp, f0, ff))
# calculate FFT for input ZAP()
inFFT = P*ZAP(Iapp, f0, ff)
output    = simulate(ZAP(Iapp, f0, ff))
# calculate the FFT of the membrane voltage response
outFFT = P*output
# calculate the impedance of the cell at each frequency ZAPfreq[i] from Z = V/I
Z = outFFT ./ inFFT
# calculate the magnitude of the complex impedance
MagZ = abs(Z)
# center the DC response (f = 0 at center of frequency axis)
shiftZ = fftshift(MagZ)
# plot the impedance at each frequency
plot(Zfreq, shiftZ, grid = false,
legend = false,
linecolor= :black,
xscale= :ln,
xlabel = "log(\$f\$)",
ylabel = "Impedance (Ω)",