-Techniques and tricks‎ > ‎Matlab‎ > ‎

### f1/f0 demo

 One can imagine several ways that a neuron might respond to a sinusoidal input, which is our normal visual input.First, let's express the sinusoidal input signal S as a function of time. We will write out 20 seconds of time in 1ms steps. Before time 0, the signal will be 0, and it will start being sinusoidal at time 0.`start_time = -10;``stop_time = 20;``t = start_time:0.001:stop_time;``f = 4; % 4 Hz frequency, or 4 cycles per second``S = (t>=0).*sin(2*pi*f*t);``plot(t,S)``box off; % turn off the box around the plot`` ``% zoom in a bit`` ``axis([-1 1 -1.5 1.5])`` ``ylabel('Signal')`` ``xlabel('Time(s)')` There are many ways that our neuron could respond to this signal. Let's first think about its voltage response. The voltage could simply fluctuate sinusoidally as a function of time.`delay = 0.1; % the response is a bit delayed in a real system like a neuron``gain = 0.5; % 0.5 response units / input signal unit``Y_1 = (t>=0+delay).*gain.*sin(2*pi*f*(t-delay));``  ``% .* means element-by-element multiplication``plot(t,S)``hold on;``plot(t,Y_1,'g');``box off; % turn off the box around the plot`` ``% zoom in a bit`` ``axis([-1 1 -1.5 1.5])`` ``ylabel('Voltage (V)')`` ``xlabel('Time(s)')`![png](output_3_0.png) If we had the output signal Y_1, how can we measure the magnitude of the response? We'd ideally like to measure the amplitude of the sinusoidal response. The technique for doing this is called Fourier analysis (Lab 3.2). We have a Matlab function that performs Fourier analysis: **fouriercoeffs_tf2**. fouriercoeffs_tf2 takes the signal, the  frequency to be examined, and the sample rate. `sample_rate = 1/0.001; % we had 1 ms steps, so the sample rate is 1000 Hz, or samples per second``A = fouriercoeffs_tf2( (Y_1((numel(Y_1)-1)/2:numel(Y_1)))', f, sample_rate)`        A =          -0.2732 + 0.4187i    The fourier coefficient at 4 Hz is a complex number that has a magnitude and phase. Let's check out the magnitude by using the **abs** (absolute value) command:`abs(A)`        ans =            0.5000    We got a value of 0.4940, which is very close to the magnitude of the acutal response (0.5). If there were no delay, then this would be exactly right. Another cell might give a response that looks like this:`Y_0 = (t>=0+delay).*gain;``  ``% .* means element-by-element multiplication``plot(t,S)``hold on;``plot(t,Y_0,'g');``box off; % turn off the box around the plot`` ``% zoom in a bit`` ``axis([-1 1 -1.5 1.5])`` ``ylabel('Voltage (V)')`` ``xlabel('Time(s)')` Here, the voltage response to the input with a fixed step. Taking the mean of the response does a good job of capturing the response.`B = mean( (Y_0((numel(Y_0)-1)/2:numel(Y_0)))')`        B =            0.5000    So, here again we have a value very close to the true function that was used to generate, with the only small difference arising from the delay in the response.We can also imagine a mixed response, like the following:`Y = Y_0 + Y_1;``plot(t,S)``hold on;``plot(t,Y,'g');``box off; % turn off the box around the plot`` ``% zoom in a bit`` ``axis([-1 1 -1.5 1.5])`` ``ylabel('Voltage (V)')`` ``xlabel('Time(s)')` We can calculate the mean of this response, and the F1 response (the amplitude of the response at the fundamental stimulus):`B_ = mean( (Y((numel(Y)-1)/2:numel(Y)))')`        B_ =            0.5000    Even though this captures the mean response very well, it does not describe the full response. It does not describe the sinusoidal modulation of the response. To capture this, we must also examine the F1 of the response:`A_ = fouriercoeffs_tf2( (Y_1((numel(Y_1)-1)/2:numel(Y_1)))', f, sample_rate)``abs(A_)`        A_ =          -0.2732 + 0.4187i            ans =            0.5000    Now, together, these two numbers (the mean, also called F0, and F1, the amplitude of the signal at the stimulus frequency) capture the response of the cell pretty well.We could also imagine that we do not have a voltage response, but, instead, we have a firing rate that is based on this voltage response. For the moment, we will let the **firing rate** of a neuron be equal to 10 spikes/(sec/volt) * Y.`FR = 10 * Y;``plot(t,FR);``box off; % turn off the box around the plot`` ``% zoom in a bit`` ``axis([-1 1 0 10])`` ``ylabel('Firing rate (spikes/sec)')`` ``xlabel('Time(s)')` We can calculate the probability that the cell fires in any given millisecond by calculating:`P = FR * 0.001;``plot(t,P);``box off; % turn off the box around the plot`` ``% zoom in a bit`` ``axis([-1 1 0 0.02])`` ``ylabel('Probability of firing')`` ``xlabel('Time(s)')` Now we can generate some spikes by drawing random numbers:`F = rand(size(P)) < P;``plot(t,S)``hold on;``plot(t,F,'r')``box off; % turn off the box around the plot`` ``% zoom in a bit`` ``axis([-1 3 0 1])`` ``ylabel('Signal (blue) and spikes (red)')`` ``xlabel('Time(s)')` We can calculate the mean response:`F_mean = sum(F)/(stop_time) % count the spikes, divide by the time`        F_mean =            5.1500    The mean response has units of spikes/sec. But knowing the mean response doesn't fully describe the whole response. We need to consider the rhythmicity of the signal as well by calculating the F1. For calculating this for spikes, we will find the spike times and compute the fourier transform: `spiketimes = t(find(F));``F1_resp_manual = (2/stop_time)*sum(exp(-sqrt(-1)*2*pi*f*(spiketimes)))``abs(F1_resp_manual)`` ``% or``F1_resp = fouriercoeffs_tf_spikes(spiketimes, f, stop_time)``abs(F1_resp)`        F1_resp_manual =          -3.2691 + 3.7591i            ans =            4.9817            F1_resp =          -3.2691 + 3.7591i            ans =            4.9817    This indicates that the response has a high degree of rhythmicity; the magnitude is the same size as the mean response. It has the units of spikes/sec.
ċ
f1_f0_demo.ipynb
(169k)
Steve Van Hooser,
Jun 1, 2020, 2:33 PM